Crop Type and Land Cover Mapping in Northern Malawi Using the Integration of Sentinel-1, Sentinel-2, and PlanetScope Satellite Data

: Mapping crop types and land cover in smallholder farming systems in sub-Saharan Africa remains a challenge due to data costs, high cloud cover, and poor temporal resolution of satellite data. With improvement in satellite technology and image processing techniques, there is a potential for integrating data from sensors with different spectral characteristics and temporal resolutions to effectively map crop types and land cover. In our Malawi study area, it is common that there are no cloud-free images available for the entire crop growth season. The goal of this experiment is to produce detailed crop type and land cover maps in agricultural landscapes using the Sentinel-1 (S-1) radar data, Sentinel-2 (S-2) optical data, S-2 and PlanetScope data fusion, and S-1 C 2 matrix and S-1 H/ α polarimetric decomposition. We evaluated the ability to combine these data to map crop types and land cover in two smallholder farming locations. The random forest algorithm, trained with crop and land cover type data collected in the ﬁeld, complemented with samples digitized from Google Earth Pro and DigitalGlobe, was used for the classiﬁcation experiments. The results show that the S-2 and PlanetScope fused image + S-1 covariance (C 2 ) matrix + H/ α polarimetric decomposition (an entropy-based decomposition method) fusion outperformed all other image combinations, producing higher overall accuracies (OAs) (>85%) and Kappa coefﬁcients (>0.80). These OAs represent a 13.53% and 11.7% improvement on the Sentinel-2-only (OAs < 80%) experiment for Thimalala and Edundu, respectively. The experiment also provided accurate insights into the distribution of crop and land cover types in the area. The ﬁndings suggest that in cloud-dense and resource-poor locations, fusing high temporal resolution radar data with available optical data presents an opportunity for operational mapping of crop types and land cover to support food security and environmental management decision-making. mapping over Turkey and found that the dual polarimetric Sentinel-1 SAR data can be used to produce accurate land cover maps. SAR products are mostly used in combination with optical images crop types and land cover for use as decision support tools to facilitate food security decision-making. Our study makes a valuable contribution to the literature on image fusion for crop type mapping. Intercropping is predominant in most smallholder contexts in SSA, but our experiments could not adequately identify the individual crops on such intercropped ﬁelds due to data constraints. Since the H/ α polarimetric decomposition image contributed to the improved accuracy of the experiments, future studies should apply H/ α polarimetric decomposition from quad-polarized data, if available, with early-season optical data and more ground truth samples to better identify the mix of crops on intercropped farms. Given the current fast pace of sensor and image processing algorithm development, it is possible to apply the method explored in this study on an operational basis to develop crop inventory data at larger scales to guide decision-making for improved food security and environmental management in cloud-dense, resource-poor locations. The approach used in this study can be used by local agriculture decision-makers to map and monitor cropping patterns and land cover change dynamics over time to build land cover inventory datasets.


Introduction
Mapping agricultural landscapes to identify crop types, analyze the spatial distribution of crops and cropping systems, and document land cover types in countries in the Global South is critical for guiding agricultural and environmental planning decisionmaking, especially in areas experiencing rapid climate change and chronic food insecurity. The type of crops and land cover on the landscape contribute to preventing soil degradation and maintaining soil health [1][2][3]. The diversity of crops on the landscape also contributes to weed control on farmlands and yield improvements by reducing the ability and Sentinel-2 data to map wetlands in Turkey while Slagter et al. [28] fused Sentinel-1 and Sentinel-2 data to map wetlands in South Africa. Kannaujiya et al. [29] integrated electrical resistivity tomography and ground-penetrating radar to map landslides in Kunjethi, while Yan et al. [30] integrated Landsat-8 optical data and Sentinel-1A to detect underground coal fires in China. Venter et al. [31] assessed the efficacy of mapping hyperlocal Tair over Oslo in Norway by integrating Sentinel, Landsat, and light detection and ranging (LiDAR) data with crowdsourced Tair measurements.
In this study, the overall goal is to examine how the integration of multitemporal dualpolarized Sentinel-1 SAR data and multispectral Sentinel-2 and PlanetScope optical data can be used for crop type and land cover mapping in complex heterogeneous agricultural landscapes using a machine learning classification algorithm. The objective is to explore the potential of several possible integrations of Sentinel-1 SAR, Sentinel-2 optical, and highresolution PlanetScope optical data with diverse data processing techniques for accurate mapping of crop types and land cover in a smallholder agricultural system. Specifically, we compare the results of classifying (i) Sentinel-1-only images, (ii) Sentinel-2-only images, (iii) the fusion of Sentinel-2 and PlanetScope images, (iv) Sentinel-2 and PlanetScope fused image + Sentinel-1 C 2 matrix, and (v) Sentinel-2 and PlanetScope fused image + Sentinel-1 C 2 matrix + H/α polarimetric decomposition image. The experiments did not consider PlanetScope-only analysis because our goal was to use the PlanetScope image to sharpen the Sentinel-2 images which have more bands, including the red edge band that is more useful for identifying the biophysical variables in vegetation [32]. The experiments are focused on two smallholder agricultural landscapes in rural northern Malawi. The hilly topography in the study locations means that cropping types and land cover characteristics tend to vary over short distances. Our final map, therefore, shows crop types and land cover categories and their distribution. Understanding the distribution of crop and land cover types in these two locations will provide a good overview of the agricultural landscape in the entire region since cropping patterns are generally similar.

Study Area Description
The study was conducted in the Edundu village area (land area =~29 km 2 ) with center location of latitude 11 • 22.545 S, longitude 33 • 46.982 E, and Thimalala village area (land area =~22 km 2 ) with center location of latitude 11 • 16.736 S, longitude 33 • 50.995 E, in northern Malawi ( Figure 1). A village area comprises several smaller communities predominantly engaged in smallholder agriculture. The study locations are in the Mzimba district in northern Malawi. Soils in the district are moderately fertile, generally medium-to lighttextured, mostly sandy-loam and loamy, with moderate to good drainage [33]. The climatic type in the area is semi-humid, with average monthly maximum temperatures ranging from 27 to 33 • C in the summer (November to April), and from 0 to 10 • C during the winter months (May to August). Annual rainfall amounts range from 800 to 1000 mm [34] which, in a good rainy season, makes the area suitable for the cultivation of a wide variety of crops, including cereals, legumes, and tobacco. The main economic activities in the area are subsistence agriculture and commercial tobacco cultivation [33]. As these smallholders predominantly rely upon rainfed agriculture, the growing season coincides with a rainy season that begins in November or December and ends in April. Maize is the most important subsistence crop cultivated in much of Malawi and is often intercropped with soybeans, groundnuts, beans, and pumpkins in the two study locations. In the last few decades, the district has often been impacted by extreme climate events such as floods and droughts, with predictions that these extremes will worsen as global climate change intensifies [33]. Understanding the dominant crop types and land cover and their distribution in the area will facilitate decision-making to build resilience to these current and projected changes in the climate. Dry season farming (dimba gardening) in valleys contributes significantly to household food security and income in the area for households that have access to the wetlands [35]. This research was part of a broader participatory interdisciplinary research project aimed at understanding the rela-tionship between farm management practices, wild biodiversity, and ecosystem services, to develop scenarios for community action plans.
intercropped with soybeans, groundnuts, beans, and pumpkins in the two study locations. In the last few decades, the district has often been impacted by extreme climate events such as floods and droughts, with predictions that these extremes will worsen as global climate change intensifies [33]. Understanding the dominant crop types and land cover and their distribution in the area will facilitate decision-making to build resilience to these current and projected changes in the climate. Dry season farming (dimba gardening) in valleys contributes significantly to household food security and income in the area for households that have access to the wetlands [35]. This research was part of a broader participatory interdisciplinary research project aimed at understanding the relationship between farm management practices, wild biodiversity, and ecosystem services, to develop scenarios for community action plans.

Satellite Data
A combination of radar (Sentinel-1) and optical (Sentinel-2 and PlanetScope) satellite images were used in this study. The sensor specifications of the data used are presented in Table 1. The Sentinel-1 constellation has two satellites-Sentinel-1A and Sentinel-1B, which were launched in April 2014 and April 2016, respectively. Sentinel-1 Single Look Complex (SLC) products, provided by the European Space Agency (ESA), consist of SAR data in the C-band and capture 5-20 m spatial resolution imagery. The two satellites have a combined revisit period of 6 days [36]. SAR images contain coherent (interferometric phases) and incoherent (amplitude features) information. The interferometric wide (IW) swatch mode, which acquires images with dual-polarization (vertical transmit, vertical receive (VV) and vertical transmit, horizontal receive (VH)), was used [36]. Abdikan et al. [37] investigated the efficiency of Sentinel-1 SAR images in land cover mapping over Turkey and found that the dual polarimetric Sentinel-1 SAR data can be used to produce accurate land cover maps. SAR products are mostly used in combination with optical images

Satellite Data
A combination of radar (Sentinel-1) and optical (Sentinel-2 and PlanetScope) satellite images were used in this study. The sensor specifications of the data used are presented in Table 1. The Sentinel-1 constellation has two satellites-Sentinel-1A and Sentinel-1B, which were launched in April 2014 and April 2016, respectively. Sentinel-1 Single Look Complex (SLC) products, provided by the European Space Agency (ESA), consist of SAR data in the C-band and capture 5-20 m spatial resolution imagery. The two satellites have a combined revisit period of 6 days [36]. SAR images contain coherent (interferometric phases) and incoherent (amplitude features) information. The interferometric wide (IW) swatch mode, which acquires images with dual-polarization (vertical transmit, vertical receive (VV) and vertical transmit, horizontal receive (VH)), was used [36]. Abdikan et al. [37] investigated the efficiency of Sentinel-1 SAR images in land cover mapping over Turkey and found that the dual polarimetric Sentinel-1 SAR data can be used to produce accurate land cover maps. SAR products are mostly used in combination with optical images to improve crop classification accuracy [38] and class discrimination [39] since they are not affected by clouds, haze, and smoke. The Sentinel-2 data used in this study were also acquired from the ESA. The Sentinel-2A (launched in June 2015) and 2B (launched in March 2017) satellites have a combined revisit period of 5 days, making them suitable for monitoring crop growth compared to other satellites, such as the Landsat, that have a lower temporal resolution. Both Sentinel-2A and Sentinel-2B satellites carry a single multispectral instrument with 13 spectral bands. We used a top of the atmosphere (TOA) reflectance product (Level-1C) provided freely by the ESA for this study. A higher-level surface reflectance product (Level-2A) can be obtained from the Level-1C using the Sentinel Application Platform (SNAP) toolbox version 7.0 [40]. Spectral bands of the Level-1C products used in this study are bands 2, 3, 4, and 8, at 10-m spatial resolution, and bands 5, 6, 7, and 8A, at 20-m spatial resolution. All the 20-m spatial resolution bands were resampled to 10-m resolution to make them comparable with the 10 m bands. The relatively higher 10-m resolution produces improved crop and land cover classification accuracy. Several previous studies have shown that Sentinel-2 data can be used to identify crop types [21,23,41].
The PlanetScope constellation of satellites presently has about 130+ CubeSats (4-kg satellites) operated by Planet Labs [42]. The majority of these CubeSats are in a sunsynchronous orbit with an equator crossing time between 9:30 and 11:30 (local solar time) [42]. PlanetScope images have four spectral bands-blue (455-515 nm), green (500-590 nm), red (590-670 nm), and near-infrared (NIR) (780-860 nm). We used the Level-3B surface reflectance products that were atmospherically corrected by Planet Labs using the 6S radiative transfer model with ancillary data from Moderate Resolution Imaging Spectroradiometer (MODIS) [42,43]. Due to its high spatial resolution, the application of PlanetScope images includes change detection [44], crop monitoring [45], and vegetation detection [46]. Unlike the other data used in this study, PlanetScope data are not free, but through the Education and Research program of Planet Labs, special permission was given for free download. All the optical images used in this study were selected based on availability as there is dense cloud cover over the study locations during the rainy season. Table 2 describes the dates of acquisition of time series satellite data used in this study. As intimated earlier, only very few cloud-free optical images were found for the study area during the growing season, a situation that affected the complexity of classification in our study and the number of experiments to be conducted. Altogether, we acquired three Sentinel-2 images captured on 7th January and 23rd February 2020, two PlanetScope images from 23rd February (Thimalala) and 1st April 2020 (Edundu), and six Sentinel-1 radar time-series images between 18th January and 13th February 2020. Field data on crop types and land cover were collected, starting in late November 2019 up to the end of April 2020. A team of trained farmer researchers went to farmlands in the two village areas with GPS devices. On each field, the coordinates of the center and the boundaries of each farm were recorded in the GPS device and on a predesigned datasheet. The type of crop(s) and the cropping system (monocrop or mixed/intercrop) Remote Sens. 2021, 13, 700 6 of 21 of each field were also recorded. The coordinates of landcover data were also recorded in the predesigned datasheets. The location of ancillary data, including vegetation, anthills, and buildings within each farm was also recorded. Figure 2 shows photographs of crops in the study locations captured during the field data collection. In all, 1668 ground truth samples of different crop and land cover types were collected. These included cereals, legumes, tubers, and vegetables, as well as roads, settlements, forests, and bare lands. The coordinates of the samples were exported from the GPS devices and the .gpx file converted to feature datasets in a geodatabase in ArcGIS Pro 2.6.3. More land cover samples were digitized from Google Earth Pro and DigitalGlobe images to complement the field data. Seventy percent (70%) of ground truth data were used as training samples, and the remaining thirty (30%) were used for validation.
Field data on crop types and land cover were collected, starting in late Novembe 2019 up to the end of April 2020. A team of trained farmer researchers went to farmland in the two village areas with GPS devices. On each field, the coordinates of the center an the boundaries of each farm were recorded in the GPS device and on a predesigne datasheet. The type of crop(s) and the cropping system (monocrop or mixed/intercrop) o each field were also recorded. The coordinates of landcover data were also recorded i the predesigned datasheets. The location of ancillary data, including vegetation, anthill and buildings within each farm was also recorded. Figure 2 shows photographs of crop in the study locations captured during the field data collection. In all, 1668 ground trut samples of different crop and land cover types were collected. These included cereal legumes, tubers, and vegetables, as well as roads, settlements, forests, and bare lands. Th coordinates of the samples were exported from the GPS devices and the .gpx file con verted to feature datasets in a geodatabase in ArcGIS Pro 2.6.3. More land cover sample were digitized from Google Earth Pro and DigitalGlobe images to complement the fiel data. Seventy percent (70%) of ground truth data were used as training samples, and th remaining thirty (30%) were used for validation.

Radar Data Processing
The Sentinel-1 IW images were preprocessed in twelve steps ( Figure 3) using th SNAP toolbox. The data were first split into sub-swaths to focus on the two study location and reduce the size of the image to improve processing efficiency. Due to the reduce sensitivity to displacement gradient when the Sentinel-1 satellite is scanning, split-band width interferometry also allows for more multilooking than is possible with standar

Radar Data Processing
The Sentinel-1 IW images were preprocessed in twelve steps ( Figure 3) using the SNAP toolbox. The data were first split into sub-swaths to focus on the two study locations and reduce the size of the image to improve processing efficiency. Due to the reduced sensitivity to displacement gradient when the Sentinel-1 satellite is scanning, split-bandwidth interferometry also allows for more multilooking than is possible with standard interferometric phase in highly deformed areas such as the hilly landscape in northern Malawi. Splitting the images further improves accuracy in low coherence areas [36]. The orbit state vectors provided in the metadata of an SAR product are generally not accurate and can be refined with the precise orbit files. These files are available days-to-weeks after the generation of the product [47]. The orbit file operator in the SNAP toolbox automatically downloads the latest released orbit file and applies it to the image so that the image can be geocoded more precisely.
Malawi. Splitting the images further improves accuracy in low coherence areas [36]. orbit state vectors provided in the metadata of an SAR product are generally not accu and can be refined with the precise orbit files. These files are available days-to-weeks a the generation of the product [47]. The orbit file operator in the SNAP toolbox autom cally downloads the latest released orbit file and applies it to the image so that the im can be geocoded more precisely. The next preprocessing step was to calibrate the images. SAR calibration produ imagery with pixel values that correlate with the radar backscatter of the captured sc [48]. Applying the radiometric correction operator also produces images that can be c pared with other SAR images acquired with different sensors or acquired from the s sensor but at different times, in different orbits, or processed by different processors. ibrating the images used in this study is thus relevant, since a multitemporal analys to be performed [47]. Creating the multitemporal multiband image requires stacking eral Sentinel-1 data collected over the growing season and fusing them with optical d The calibration process was followed by debursting. Sentinel-1 IW SLC products con of one image per swath per polarization. The IW products that were used have th swaths. Each sub-swath image consists of a series of bursts. Each burst was processe a unique SLC image. To merge all these bursts into a single SLC image, the TOPSAR burst and Merge operator in the SNAP toolbox was applied. The polarimetric matrix g eration process was used to generate a covariance matrix from the Sentinel-1 C-band S data [49]. Dual polarimetric SAR sensors collect half of the total polarimetric informa involved in fully polarimetric imagery or quad polarization [50]. This implies that the olution cell at each time point is defined by a 2 × 2 covariance (C2) matrix that is obtai from C3 (representing the average polarimetric information extracted from a set of ne boring pixels). The resulting C2 matrix is represented by Equation (1)  The next preprocessing step was to calibrate the images. SAR calibration produces imagery with pixel values that correlate with the radar backscatter of the captured scene [48]. Applying the radiometric correction operator also produces images that can be compared with other SAR images acquired with different sensors or acquired from the same sensor but at different times, in different orbits, or processed by different processors. Calibrating the images used in this study is thus relevant, since a multitemporal analysis is to be performed [47]. Creating the multitemporal multiband image requires stacking several Sentinel-1 data collected over the growing season and fusing them with optical data. The calibration process was followed by debursting. Sentinel-1 IW SLC products consist of one image per swath per polarization. The IW products that were used have three swaths. Each sub-swath image consists of a series of bursts. Each burst was processed as a unique SLC image. To merge all these bursts into a single SLC image, the TOPSAR Deburst and Merge operator in the SNAP toolbox was applied. The polarimetric matrix generation process was used to generate a covariance matrix from the Sentinel-1 C-band SAR data [49]. Dual polarimetric SAR sensors collect half of the total polarimetric information involved in fully polarimetric imagery or quad polarization [50]. This implies that the resolution cell at each time point is defined by a 2 × 2 covariance (C 2 ) matrix that is obtained from C 3 (representing the average polarimetric information extracted from a set of neighboring pixels). The resulting C 2 matrix is represented by Equation (1).
Dual polarization imagery has only diagonal elements. As such, the matrix with off-diagonal components were set to zero and do not follow a complex Wishart distribution; but, the two diagonal blocks (1 by 1) follow a complex Wishart distribution [51,52]. Multilooking and polarimetric speckle filtering were applied to reduce speckle noise. The presence of speckle intensity fluctuations in SLC SAR imagery is the result of the nature of coherent image formation during SAR data collection [47]. Each radar resolution cell contains multiple scatters, each of which contributes to the overall signal returned from the resolution cell. The phase obtained from each scatter is effectively random because the radar wavelength is normally much smaller than the size of the resolution cell. The signals from each scatter may be summed according to the principle of superposition, resulting in constructive and destructive interference [53]. Low reflectivity occurs in cells where destructive interference dominates, while constructive interference dominates in cells where high reflectivity prevails. Consequently, the phenomenon of speckle occurs on the image. Multilooking, which is achieved by dividing the signal spectrum and then incoherently averaging the recovered sub-images of an SAR image, is widely used to reduce these speckles in conventional SAR signal processing [54]. The Boxcar polarimetric speckle filtering method with a 5 × 5 window size was used to further reduce speckle noise in the image while preserving the complex information of all the bands, enhancing interpretation, and improving their ability for quantitative analysis. The Boxcar filter is the simplest filter that locates similar pixels by a moving window with the predefined size. The filter reduces the speckle phenomenon while producing a blurring visual effect [55].
Polarimetric decomposition was used to separate different scattering contributions and provide information about the scattering process [56]. Polarimetric decompositions are techniques used to generate polarimetric discriminators that can be used for analysis, interpretation, and classification of SAR data [57]. The H/α dual-polarized decomposition of VV-VH dual-polarization with the window size of 5 × 5 was performed. H/α polarimetric decomposition is an entropy-based method based on the theory that the polarization scattering characteristics can be represented by the space of the entropy and the averaged scattering angle α employing the eigenvalue analysis of Hermitian matrices [52]. The H/α polarimetric decomposition method proposed by Cloude and Pottier [52] is widely used for land cover classification and object recognition. The method has good properties such as rotation invariance, irrelevance to specific probability density distributions, and covers the whole scattering mechanism space [58]. Only the 2 × 2 covariance matrix was derived, because the Sentinel-1 product is dual-polarized. The entropy and alpha images were derived from the H/α dual-polarized decomposition and used for the present analysis.
During range-Doppler terrain correction [59], a Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) was downloaded automatically from the SNAP toolbox to fulfill orthorectification [60]. Terrain-correcting the images is important because the discrepancies in topographical variations of a scene and the tilt of the satellite sensor distort distances in SAR images. Applying terrain correction compensates for these distortions and make features appear as close as possible to real-world features [59]. The final image was resampled into 3-m to conform with the optical PlanetScope data. A subset of all the images is clipped and a stack was then made of all the layers.

Optical Data Processing
For the Sentinel-2 and PlanetScope data, we applied atmospheric correction, topographic correction, resampling, band stacking, seamless mosaicking, and image sub-setting as preprocessing steps. The Sentinel-2 Level-2A bottom-of-atmosphere (BOA) products we acquired are already atmospherically corrected. Examples of preprocessed natural color and grayscale images of the radar and optical data of the study locations are shown in Figure 4.

Image Fusion
Images from different spectral bands have the same geometric information [61]. Based on this principle of satellite images, Gašparović et al. [62] developed an image fusion method known as the P + XS fusion, in which an image is perceived as a function whose sampling corresponds to the discrete matrix form of the image [63]. The P + XS method introduces the geometry information of a higher resolution image by aligning all edges of the higher resolution image with each lower resolution multispectral band. To obtain the spectral information for the fused image, the method assumes that images captured in different spectral bands share common geometric information and that the higher resolution image can be approximated as a linear combination of the high-resolution multispectral bands [64,65]. Guided by this principle, each Sentinel-2 band was fused with the corresponding high spatial resolution PlanetScope band that shares similar spectral characteristics. He et al. [65] observed that using the P + XS method can better preserve sharp discontinuities such as edges and object contours on an image. The objective of the fusion was to produce higher resolution Sentinel-2 multispectral fused images from the original low-resolution Sentinel-2 and high-resolution PlanetScope images.

Image Fusion
Images from different spectral bands have the same geometric information [61]. Based on this principle of satellite images, Gašparović et al. [62] developed an image fusion method known as the P + XS fusion, in which an image is perceived as a function whose sampling corresponds to the discrete matrix form of the image [63]. The P + XS method introduces the geometry information of a higher resolution image by aligning all The smoothing filter-based intensity modulation (SFIM) [66] was used to fuse the Sentinel-2 and PlanetScope images. The SFIM aims to produce fused images that have the highest spatial resolution in the same multispectral bands of the original low-resolution images. The lower spatial resolution Sentinel-2 bands were fused with the higher spatial resolution PlanetScope bands based on the principles of the P + XS method. The corresponding bands of the PlanetScope were alternatively used as the high spatial resolution band in place of panchromatic bands ( Table 3). The SFIM technique improves spatial details while keeping the spectral properties of the images. The digital number (DN) value of the fused band is defined as: where FUS i is the DN value of the fused band i, X i is the corresponding lower spatialresolution band I, resampled to the same high resolution as Y, Y represents the high spatial resolution panchromatic band, and Y represents the Y band being averaged by the low-pass filtering. For Thimalala, the Sentinel-2 images acquired on 07/01/2020 were fused with the PlanetScope images on 23/02/2020. For Edundu, the Sentinel-2 images acquired on 07/01/2020 and 23/02/2020 were fused with the PlanetScope images taken on 01/04/2020 separately. Figure 5 shows examples of subset image cropped from the originally preprocessed Sentinel-2 images and the fused (Sentinel-2 + PlanetScope) image. It can be observed from the images that the fused data are clearer than the original ( Figure 5).

Image Classification
Fourteen (14) crop types in total and eight (8) land cover classes identified in the two study locations were used as the schema of the classification. The crop type classes were identified through field data collection. The land cover types were identified based on field observations during ground truth data collection, as well as from very high resolution Google Earth Pro and DigitalGlobe images. The supervised random forest (RF) machine learning algorithm [67] in ENVI 5.3 was used for the classification. As a non-parametric method, RF benefits from including categorical and continuous datasets based on well-developed rules and does not require training data to come from a unimodal distribution [68]. An RF is generated through the creation of a series of decorrelated decision trees using bootstrapping, thus solving the problem of overfitting and producing accurate results [67]. Tuning parameters, such as the number of trees and the number of split candidate predictors, are generally chosen based on the out-of-bag (OOB) prediction error. RF uses the OOB samples for cross-validation, and once the OOB errors stabilize at a reasonably large number of trees, training can be concluded. In this study, the model we used in each experiment consisted of 500 trees and the number of features to split the nodes was set to half of the features within the corresponding input dataset. Accuracy of RF classification is often very high, even when based on many input features such as multitemporal images [69].

Image Classification
Fourteen (14) crop types in total and eight (8) land cover classes identified in the two study locations were used as the schema of the classification. The crop type classes were identified through field data collection. The land cover types were identified based on field observations during ground truth data collection, as well as from very high resolution Google Earth Pro and DigitalGlobe images. The supervised random forest (RF) machine learning algorithm [67] in ENVI 5.3 was used for the classification. As a non-parametric method, RF benefits from including categorical and continuous datasets based on well-developed rules and does not require training data to come from a unimodal distribution [68]. An RF is generated through the creation of a series of decorrelated decision trees using bootstrapping, thus solving the problem of overfitting and producing accurate results [67]. Tuning parameters, such as the number of trees and the number of split candidate predictors, are generally chosen based on the out-of-bag (OOB) prediction error. RF uses the OOB samples for cross-validation, and once the OOB errors stabilize at a reasonably large number of trees, training can be concluded. In this study, the model we used in each experiment consisted of 500 trees and the number of features to split the nodes was set to half of the features within the corresponding input dataset. Accuracy of RF classification is often very high, even when based on many input features such as multitemporal images [69].

Post-classification Processing
A 3 × 3 filter with the majority vote was applied to remove the salt-and-pepper effect in the classified images. Accuracy assessment was performed to evaluate the performance of each classification experiment using the ground truth data acquired from fieldwork, Google Earth Pro, and DigitalGlobe images based on Congalton [70]. The overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA), and the Kappa coefficients were computed for each of the five experiments to determine the bestperforming model for identifying crop types and land cover categories. The OA describes the proportion of pixels correctly classified, with 85% or more generally accepted as the threshold for a good classification. The PA is the map accuracy from the point of view of the map producer. It explains how often features on the ground are correctly shown on the classified map, or the probability that a certain land cover type on the ground is classified as such [71]. The UA, on the other hand, is the accuracy from the point of view of a map user. It describes how often the features identified on the map will be present on the ground [71]. The Kappa coefficient is the ratio of the agreement between the classifier output and reference data, and the probability that there is no chance agreement between the classified and the reference data [70]. Table 4 presents the results of the accuracy assessment for all five experiments. For the Sentinel-1-only experiment, the overall accuracies of both locations were lower than 50%. For the Sentinel-2-only experiment, the accuracies increased significantly to 72% and 74%, respectively, but were still lower than expected for an ideal classification. The water class was better extracted in the Sentinel-2-only experiment, compared with the Sentinel-1-only experiment. When using the fused image of Sentinel-2 and PlanetScope, the introduction of the relatively higher spatial resolution image of 3 m added more detail to the classes, though it failed to accurately identify groundnut crops. The overall accuracy increased to 76.03% for Thimalala and 84.12% for the Edundu location. Table 4. Comparison of the accuracy assessment results of the random forest classification.

No. Data Combination Thimalala Edundu
Overall Accuracy (%) Kappa Coefficient Overall Accuracy (%) Kappa Coefficient In the fourth experiment, we made a stack of the fused Sentinel-2 and PlanetScope images and the Sentinel-1 C 2 matrix (Sentinel-2 and PlanetScope fused image + Sentinel-1 C 2 matrix). The addition of the Sentinel-1 C 2 matrix generated minor improvements on the Sentinel-2 + PlanetScope experiment in experiment (iii), as observed in the improved accuracy statistics ( Table 4). The combination of optical and SAR data resulted in a slightly better increment (5%) in the overall accuracy for the Thimalala site, compared to the Edundu site (0.42%), but the accuracies were still lower than the accepted value for a good classification.
In the final experiment, a multiband image was created comprising the fused image from the fourth experiment and the H/α polarimetric decomposition (i.e., the fusion of Sentinel-2 + PlanetScope + Sentinel-1 C 2 matrix + H/α polarimetric decomposition). This final experiment produced results that outperformed all the other integration of images. The overall accuracies from this final experiment exceeded 85% in both locations, which is the threshold for a good classification. Overall, the accuracies improved by 13.53% and 11.7% for Thimalala and Edundu, respectively, compared with the Sentinel-2-only classification. Table 5 shows the PA and UA accuracies for the experiment (v). Crops such as bambara nuts (100%), groundnuts (90%), maize (100%), tobacco (87.5%), and tomato (100%) had high PAs with equally high UAs, meaning the method accurately identified the crops. The corresponding UAs also showed that the maps are useful for identifying these common crops on the landscape from the user's perspective. Crops including banana, onion, and sweet potato were very few in the landscape and only a few samples were obtained. As such, our model could not adequately assess these classes. This explains the nature of the PAs and UAs for these crop classes (Table 5). Land cover classes such as forest, shrubland, untarred road, and water also had high PAs, compared to settlement, bare rock, and tarred roads, which had low PAs. Corresponding UAs were within a similar range to the PAs. Table 5. Summary of classification accuracies of Thimalala using the fused image of Sentinel-2 and PlanetScope, Sentinel-1 C 2 matrix, and H/α polarimetric decomposition.

Land Cover and Crop Class
Producer's Accuracy (%) User's Accuracy (%)  Figures 6 and 7 show the maps generated for the various experiments for the two locations. Figures 6a and 7a show results for the Sentinel-1-only experiment for Thimalala and Edundu, respectively, while Figures 6b and 7b show the results for the Sentinel-2-only experiment for the two village areas. The results reflect the poor accuracies generated from the experiments (see Tables 4 and 5). Results presented in Figures 6c and 7c show that the various classes are more identifiable than in the foregoing experiments in both locations, while Figures 6d and 7d reflect the minor improvement of accuracies when moving from the experiment (iii) to (iv) (see Table 4). It can be seen in Figures 6e and 7e that the classification improved with the introduction of the H/α decomposition layer. In terms of the distribution and dominance of crop types, both Figures 6 and 7 show that maize is the dominant crop type in both study locations. Tobacco crops are also common in both areas and are observed to have, on average, larger sizes than most food croplands. Intercropped farmlands, which are mainly maize intercropped with other legumes and pumpkins, are also prominent on the landscape in both locations. Though the classification was able to identify intercropped fields, most likely maize with other crops that are equally tall, such as pigeon pea or cassava, it was, however, not able to identify which specific crop combinations (e.g., maize/bean and maize/soybean) constitute each intercropped farm. Croplands that were not classified into any of the individual categories were grouped into the mixed/intercropped class. As expected, croplands are mostly in low-lying valley areas. Shrublands and forests are the most dominant land cover types on hills and hillslopes in both locations. The forests are mostly surrounded by shrublands interspersed with some farmlands, especially in the Thimalala area. Overall

Discussion
The study reveals useful observations about crop type and land cover identification in heterogenous smallholder agricultural landscapes. The experiments provide evidence that to effectively identify crop types, cropping systems, and land cover using pixelbased classification, a combination of multitemporal satellite data from different sensors can be used. We show that the fusion of Sentinel-2 and PlanetScope optical images integrated with Sentinel-1 C 2 matrix and Sentinel-1 H/α polarimetric decomposition was the most effective in generating high-accuracy crop type and land cover maps ( Table 4). The degree of accuracy increased from the experiment (i) as the number of images integrated for the classification increased. Using only the Sentinel-1 images produced OAs lower than 50% for both locations, similar to observations made by Mercier et al. [72]. When the Sentinel-1 C 2 matrix was fused with Sentinel-2, the accuracies improved for both study locations. Fusing the two optical images further improved the classification results of the Sentinel-2-only by 3.95% and 10.04% for Thimalala and Edundu, respectively. Orynbaikyzy et al. [73] noted similarly that fusing Sentinel-1 and Sentinel-2 data improved the results of crop classification in Northern Germany, as did Mercier et al. [72] in Paragominas (Brazil). Gašparović et al. [62] observed a similar improvement when they integrated Sentinel-2 with PlanetScope for vegetation mapping and monitoring. There were only a few banana, onion, and sweet potato samples in the validation data because these crops are generally scarce in the agricultural landscape of both study locations. They are also often planted on relatively smaller plot sizes. A such, the PAs and UAs for these crop classes show that the assessment might have been less inaccurate ( Table 5). The foregoing explanation is one of the main limitations of classifying crop types in complex heterogeneous smallholder agriculture landscapes, as highlighted in other studies [8]. At the same time, the observation presents an opportunity for future studies to explore other approaches to capturing such rare crops using image classification.
The observation that fusing Sentinel-1 and the optical data improved crop type mapping is also consistent with findings of other studies that have combined radar and optical data for crop and land cover type mapping [22,65]. The finding provides further evidence of the contribution of dual-polarized Sentinel-1 data for accurate crop and land cover classification. We conclude that the 13.53% (Thimalala) and 11.7% (Edundu) improvements in overall accuracy on the Sentinel-2-only classification is because of the integration of the multitemporal datasets from all the various sensors with the Sentinel-1 H/α polarimetric decomposition. This finding suggests that adding more multitemporal images could further improve the OAs, PAs, and UAs, and using only a unitemporal image produces underwhelming results. Previous studies have also found that integrating the H/α polarimetric decomposition information with other data achieves better accuracy in complex agricultural landscapes than other classical methods [74][75][76]. Dual-polarized data such as Sentinel-1 are a valuable data resource for decomposing radar data to map crop types [73]. The poor accuracy from the Sentinel-1-only experiment (i) is, therefore, likely the result of the design of our experiment (i). The two bands of the Sentinel-1 IW mode compared to the several bands available from Sentinel-2 and PlanetScope and the use of the random forest algorithm may explain the poor accuracy from the Sentinel-1-only classification [77]. Quad-polarized data are known to perform better in identifying crop types when integrated with other data [78]. Future crop type mapping in the study locations should, therefore, integrate time series quad-polarized H/α polarimetric decomposition, if available, with available time series optical data to achieve better accuracies.
The combination of multitemporal, multispectral, and multisensor data to map crops in such heterogeneous landscapes suggests that in areas of high cloud cover where optical data collection is not feasible, combining data with different spectral characteristics, including radar and optical data, holds the potential for reliable mapping of crop types and land cover for building crop type inventories. The classification outcomes indicate that the multitemporal routine of fusing the quality high-resolution optical images (Plan-etScope) with radar data and the random forest classification approach outperformed the Sentinel-1-only data. Malawi, similar to many other locations in the tropics, has dense cloud cover during critical stages of crop growth but has limited coverage of high temporal resolution remote sensing satellites. This makes it difficult to acquire time series optical data for monitoring crop growth and mapping crop diversity. Clouds are often cited as the main advantage of using radar for mapping land cover. Being able to routinely map the landscapes in such resource-poor contexts using a combination of images from different sensors is critical to creating crop inventory maps by local officials to facilitate food security and environmental management decision-making [79]. The timing and quality of satellite observations play a crucial role in the accuracy of classification to identify crops [77]. Even though the seasonal satellite data we used contributed to the overall accuracy, only a few of the cloud-free optical images were obtained for the start of the growing season due to heavy cloud cover. As such, it was difficult for our optimal model to effectively separate the different crop combinations on intercropped farms, even though some intercropped fields were identified. The experiment (v) was able to identify the distribution of crops, with maize being the most common crop, a finding consistent with the maizification of Malawi narrative found in other studies, a phenomenon attributed to government agricultural policies [79]. Future studies can overcome the challenge of inaccurate identification of intercropped fields by using images acquired at the start of the growing season when the specific crop combinations on intercropped fields are visible to satellite sensors and can be captured and separated by image classification algorithms.
A combination of field ground truth data and samples from Google Earth Pro and DigitalGlobe images, as well as ancillary information, was used to train and validate the random forest algorithm to attain more than 85% accuracy in the experiment (v) ( Table 4). The accuracy attained suggests that using reference data from diverse sources can improve classification results, as observed in other studies [23,74]. This observation implies that the level of accuracy obtained in our experiments can be further improved with more training and validation samples. In future, crowdsourcing should be used to collect more field samples by training farmers to use smartphones and other spatial data collection applications and tools to geolocate and record farm and ancillary information for training and validating classifiers. Using crowdsourcing to collect crop type samples has been experimented with in other contexts and found to have contributed to getting large volumes of training data for mapping crop types more accurately [23]. Adding more training samples will also allow classification algorithms to identify the specific crops in intercropped fields.

Conclusions
The outcomes of the experiments conducted in this study highlight the importance of exploiting the capabilities of various satellite sensors to create high temporal resolution images for mapping crop types and land cover in smallholder agriculture areas where the landscapes tend to be more heterogeneous. Both individual Sentinel-1 and Sentinel-2 images failed to produce high-accuracy crop type and land cover maps. The overall accuracy and Kappa coefficients of the classification improved as the number of images increased and the spatial resolution improved. The fusion of Sentinel-1 C 2 matrix, Sentinel-2, and PlanetScope optical data with the Sentinel-1 H/α polarimetric decomposition outperformed all other combinations of images. Fusing the images created high temporal resolution data, with Sentinel-1 contributing the greatest number of images due to the ability of SAR to penetrate cloud cover. Though processing time may increase due to the high volume of data being integrated, the fact that acceptable accuracies were achieved is crucial. In the tropical areas of SSA where cloud cover is often dense during the growing season, this study demonstrates that computing the H/α polarimetric decomposition from cloud-penetrating radar data and fusing it with other high-and moderate-resolution optical data can be very cost-effective for developing large-scale crop inventory data. Several studies have shown how these current and emerging sensors can be harmonized to map crop types and land cover for use as decision support tools to facilitate food security decision-making. Our study makes a valuable contribution to the literature on image fusion for crop type mapping. Intercropping is predominant in most smallholder contexts in SSA, but our experiments could not adequately identify the individual crops on such intercropped fields due to data constraints. Since the H/α polarimetric decomposition image contributed to the improved accuracy of the experiments, future studies should apply H/α polarimetric decomposition from quad-polarized data, if available, with earlyseason optical data and more ground truth samples to better identify the mix of crops on intercropped farms. Given the current fast pace of sensor and image processing algorithm development, it is possible to apply the method explored in this study on an operational basis to develop crop inventory data at larger scales to guide decision-making for improved food security and environmental management in cloud-dense, resource-poor locations. The approach used in this study can be used by local agriculture decision-makers to map and monitor cropping patterns and land cover change dynamics over time to build land cover inventory datasets.