Sophisticated cropping patterns mapping in Kenya using vegetation indices and phenology from Sentinel-2 and Sentinel-1 radar backscatter data

The proportion of area under various crops at a given point in time, known as a cropping pattern, plays an essential role in determining the level of agricultural production. In this study, cropping patterns of three sub-counties in Murang’a County, a typical African smallholder farming area in Kenya, were mapped. Specifically, we compared the performance of eight classification scenarios for mapping cropping patterns; namely using (i) only Sentinel-2 reflectance bands (S2), (ii) S2 and S2 derived vegetation indices (VIs); (iii) S2 and S2 vegetation phenology (VP); (iv) S2 and Sentinel-1 radar backscatter data (S1); (v) S2, VIs, and S1; (vi) S2, VP, and S1; (vii) S2, VIs and VP, and (viii) S2, VIs, VP and S1. Reference data of the dominant cropping patterns and non-croplands were collected. The guided regularized random forest (GRRF) algorithm was used to select the optimum variables and to perform the respective classification for each scenario. The most accurate result of the overall accuracy of 93.16% was attained from the scenario (viii) S2, VIs, VP, and S1. The McNemar’s test of significance did not show significant differences (p≤0.05) among the tested scenarios. Our study demonstrated the strength of GRRF and the synergetic advantage of S2 and S1 derivatives to map cropping patterns in a heterogeneous landscape where high resolution imagery are inaccessible. Our cropping pattern mapping approach can be used in other sites of relatively similar agro-ecological conditions.


Introduction
Cropping patterns and food security in Africa are highly affected by a myriad of factors including climate change and variability, and high pests and disease infestation levels, thus causing low production rates [1]. The proportion of area under various crops at a given point in time, known as a cropping pattern [2], plays an essential role in determining the level of agricultural production, particularly for small-scale farmers in Africa [3]. Despite the importance of coherent cropland information, studies have indicated that the small-scale structures (< 0.1 ha) and fragmentation of these cropping patterns in Africa, triggered by their high intra-and inter-seasonal dynamics, prohibit their accurate detection and characterization [4]. This complexity is further aggravated by the individual farmers' crop selection for planting in a particular season. The choice of crop for planting 3. Methodology Figure 2 illustrates the methodological approach used to map the cropping patterns in the study area using different remotely sensed variables extracted from S2 and S1 imageries.

Remotely sensed data
The remotely sensed data utilized in this study included multi-date reflectance bands, vegetation indices, phenological metrics from S2, and backscatter from S1. These freely available datasets of S2 and S1 were selected to emphasize the various parts of the electromagnetic spectrum for accurate cropping pattern mapping. 3

.1.1. Sentinel-2
Freely available multi-date S2A level 1C imagery (n=126 scenes) captured during four seasons from 10 December 2017 to 15 December 2018 (Table 1) with cloud cover less than 20% were acquired from the European space agency (ESA) Copernicus open access hub [41]. S2 has a swath width of 290 km 2 providing 13 spectral bands that range in pixel size from 10, 20, and 60 m across the visible, near-infrared, and the shortwave infrared spectrum [41], as shown in Table 2. The S2 images were atmospherically corrected using the Sen2cor module in the sentinel application platform (SNAP) toolbox [41]. In SNAP, we also performed cloud masking, resampling of the S2 20 m spectral bands to 10 m spatial resolution, layer stacking, and computation of the median pixel image for each season. We then used the Murang'a sub-county boundaries to subset the respective median images to our study area. Coastal and aerosol (band 1), water vapor (band 9) and cirrus (band 10) bands ( Table 2) were excluded from this analysis as they mostly contribute to atmospheric and geophysical parameters [41], which were not the focus of our current study.

. Vegetation indices
VIs are a combination of spectral characteristics of two or more wavelength bands that indicate the relative abundance of vegetation components such as chlorophyll and water contents [42]. The VIs described in Table 3 were selected and derived for this study from the median images of our multi-date S2 imagery. Composite images of each of the VIs (n=8) were then created. These VIs were selected for cropping patterns mapping as they are the most robust indices that mimic vegetation seasonality and are known to reduce residual contamination due to atmospheric noise and soil background [26]. The multi-date and multi-season aspects in the VIs can also cater for NDVI anisotropic effects [43].  [52]. The VP variables as described by Araya [53] were derived in the present study and are summarized in Table 4. The VP variables were simulated from the multi-date NDVI curve of the S2 images of each of the four seasons (Table 1). A composite image for the VP variables (n=15) was created and used in the cropping pattern classification experiment.

No.
Phenological variable Definition of the normalized difference vegetation index (NDVI) curve and physiological description 1 Onset_value The NDVI value at the start of the growth (seedling growth stage) 2 Onset_time The time when the growth onset is achieved 3 Max_value The maximum NDVI value in the season 4 Max_time The time when the max_value is attained (anthesis growth stage) 5 Offset_value The NDVI value at the end of the season 6 Offset_time The time when growth offset is attained (senescence growth stage) 7 LengthGS The length of the growing season 8 BeforeMaxT The length of time between onset and max_value 9 AfterMaxT The length of time between max_value and offset 10 GreenUpSlope The rate of increase in NDVI value between onset and offset 11 BrownDownSlope The rate of decrease in NDVI value between max_value and offset 12 TINDVI The area under the NDVI curve between onset and offset 13 TINDVIBeforeMax The area under the NDVI curve between onset and max_value 14 TINDVIAfterMax The area under the NDVI curve between max_value and offset 15 TINDVIAsymmetry The difference between BeforeMaxTINDVI and AfterMaxTINDVI  [41] and used in this analysis ( Table 5). The acquisition dates of S1 imagery were slightly different from those of S2 due to image availability in the ESA archive. S1 sensor provides C-band synthetic aperture radar (SAR) images in both singular and dual-polarization with a repeat cycle of 12 days [41]. Acquisition of these images are in four modes i.e., stripmap (SM), interferometric wide swath (IW), extra-wide swath (EW), and wave (WV) with different processing levels i.e. Level-0, Level-1 (Single Look Complex-(SLC), ground range detected-(GRD)) and Level-2 [41]. In this study, Level-1 S1 products of GRD and IW were used. The S1 images were dual-polarized in vertical transmit and vertical receive (VV) and vertical transmit and horizontal receive (VH) mode. The pre-processing procedures of S1 images were performed in the SNAP toolbox [41]. Based on Filipponi [54], these processes included applying the precise orbit file (provides accurate satellite position and velocity information), radiometric calibration using calibration beta, and speckle filtering using the lee sigma filter. Similarly, we computed the terrain correction using 20 m shuttle radar topography mission (SRTM) data and resampled the data to 10 m pixel size. These processed images were then stacked to produce per season median pixel value of the VV and VH image bands that were then subset to the extent of the study area.

Field data collection
A stratified random sampling approach was used for the field data collection, using three vegetation intensity classes (i.e. low, medium, and high) as strata at a landscape scale using the K-means clustering method [40]. The K-means clustering method groups similar data points together and uncovers underlying patterns [55]. The stratified random sampling was done to provide a comprehensive representation of the cropping patterns within different vegetation intensity in the study area at a landscape scale. Field reference data for the cropping patterns and non-croplands were collected during the period ranging from 13 December 2018 to 19 December 2018, corresponding to the short wet season. The cropping patterns classes included monocrop avocado, monocrop coffee, monocrop maize, monocrop tea, monocrop pineapple, mixed crop avocado, and mixed crop maize. The non-cropland classes included water bodies, built-up area, grassland, shrubland, and natural vegetation. The non-cropland classes were included due to the heterogeneity of the landscape in the study area. Reference points were field sampled randomly within each of the three vegetation intensity strata, ensuring that all classes had more than thirty data entries well spread across the strata. The onfield data of each class were collected as points (i.e. pixels), which were then converted to homogenous units (i.e. polygons) using on-screen digitization on high-resolution Google Earth imagery [56]. The pixel values within the homogenous polygons were then extracted to be utilized for the cropping pattern predictions as summarized in Table 6.

Predictor variables selection and classification
The GRRF algorithm was used to select the most important variables in eight individual classification scenarios (Table 7). GRRF algorithm works with using a concept similar to that of the random forest (RF) method. However, RF uses an adaptable number of predictor variables and does not adhere to the assumption of a normal distribution of the predictor variables [57], whereas GRRF additionally uses the importance scores from a preliminary RF to guide the feature selection process of regularized random forest (RRF). In RRF, a feature is selected by building only one ensemble where the feature importance is evaluated on a part of the training data at each node [33]. However, the limitation of RRF feature selection is that many features can share the same information gain at a node with a small number of instances and a large number of features, causing a likelihood of RRF to select a feature that is not strongly relevant. To counter for the limitations of RRF, GRRF assigns a penalty coefficient to each feature by changing the importance coefficient of gamma value (γ ∈ (0, 1)), which controls the weight of normalized importance. This ensures that the most relevant variables are retained while still ensuring accurate classification [58]. Thus, in this study γ=0.7 was used in the "CoefReg" function in the "caret" package in R software to select the most important variables [59]. The mean decrease in accuracy score value depicts the level of importance of the variable. We set the number of decision trees grown (ntree) to 500 default value and the number of variables used to split the trees (mtry) to the square root of the number of variables, which is also a default setting [57]. Furthermore, we used the flagReg=1, which initiates the penalization of the variables at each node and to determine the most important variables with a tune length of 3. The tune length parameter controls the algorithm to attempt different values for the main parameter in the prediction. The eight individual classification scenarios (C1-C8) combined different multi-date S2 and S1 variables ( Table 7) and their counterparts that combined only the selected S2 and S1 variables i.e. C1selected-C8selected were then carried out using the GRRF algorithm.

Classification accuracy assessment
We used a 2-fold cross-validation method to train and validate the performance of the classification scenarios [60]. In this approach, each fold had a chance to be a training and testing set. The performance of each classification scenario was then evaluated using the area under class method [61]. The area under class method establishes the unbiased total area of each class since it includes the area of the map omission error of each class, leaving out the commission error. The producer's accuracy (PA), user's accuracy (UA), and overall accuracy (OA) were then computed from the classification error matrix of estimated area proportions, taking into account the proportion of each class in the study area. Moreover, a class-wise accuracy assessment was performed for each class using the F1-score criterion [62]. F1-score is a measurement that balances the difference between PA and UA for each class i through the formulation of the harmonic mean of PA and UA as shown in Equation 1.
The unbiased estimated area of each class and a confidence interval of 95% [61] were also calculated from the most accurate scenario. Additionally, a McNemar's chi-square test [63] was carried out to test for any statistically significant differences (p ≤ 0.05) among the cropping patterns mapping results from the classification scenarios.

Variable selection using the guided regularized random forest (GRRF) algorithm
S2 reflectance bands across the four seasons in all the GRRF scenarios (Table 8) dominated in the most relevant bands selected. Specifically, the S2 bands of NNIR, SWIR3, RE1, Red, Green, and Blue were selected. Additionally, NDWI and GNDVI were the most relevant VIs selected from the eight indices investigated. The most relevant VP variables selected were the Offset_value, TINDVIBeforeMax, and TINDVIAfterMax, whereas the S1 bands of VV and VH were selected. The specific variables selected in each of the eight scenarios are shown in Figure 3.

Mapping accuracy assessment
The OA of all the classification scenarios was above 90.00% as shown in Figure 5. Among the cropping patterns, the combined scenarios C1-C8 (Table 7) showed the highest PA (100.00%) and F1score (99.97%) for monocrop pineapple while the highest UA (100%) was recorded for monocrop tea. On the other hand, the lowest PA (86.60%) was shown for monocrop coffee from scenario S2 and VP (C3), the lowest UA (69.63%) for mixed crop avocado from scenario S2, VP and V1s (C7), and the lowest F1-score (80.86%) for mixed crop avocado from scenario S2, VP, and VIs (C7).
Regarding the GRRF selected scenarios (C1selected-C8selected), the highest PA (100.00%) was shown for monocrop pineapple from the GRRF selected variables of S2 (C1selected), S2 and VIs (C2selected), S2 and VP (C3selected), and S2 and S1 (C4selected). Monocrop tea had the highest UA (100.00%) from the GRRF variables of S2 and VIs (C2selected), and S2 and VP (C3selected), while the highest F1-score (99.97%) Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 26 October 2020 doi:10.20944/preprints202010.0517.v1 was shown for monocrop pineapple from the GRRF selected variables of S2 and VIs (C2selected). In addition, the lowest PA (82.29%) was predicted for mixed crop maize in the GRRF selected variables of S2 (C1selected). The lowest UA (67.02%) and F1-score (78.30%) were recorded for mixed crop avocado from the GRRF selected variables of scenario S2 and V1s (C2selected). In as much as the non-croplands were not a focus of this present study, water was mapped with class-wise accuracy of 99.00% across all scenarios. A summary of the class-wise accuracies of all scenarios is shown in Figures 6 and 7. However, the McNemar's test for significance showed no statistical significance (p ≤ 0.05) difference between the combined scenarios and GRRF selected scenarios (Table 8).

Area estimation
The unbiased area estimate calculated from the best GRRF scenario i.e. C8selected among the cropping patterns revealed that mixed crop avocado had the greatest acreage (13,940.18 ha) while the lowest acreage in the cropping patterns was the monocrop avocado (1,929.92 ha) class. In noncroplands, shrubland had the greatest acreage (21,594.77 ha) while water had the lowest acreage (1,167.00 ha). The acreages of other cropping patterns and non-croplands are summarized in Figure  8.

Discussion
This study illustrated the utility of multi-date and multi-sensor variables for classifying cropping patterns for a heterogeneous agro-natural landscape in Kenya. We tested the performance of eight classification scenarios that resulted from combining different vegetation indices and vegetation phenological variables from S2 spectral and S1 SAR backscatter imagery. Seven cropping patterns were classified. The specific combination of the selected variables in each scenario using the GRRF algorithm was also analyzed. We utilized satellite-based data acquired over four seasons (i.e. short and long dry as well as short and long wet season) that were assumed to have captured all the changes in inter-and intra-vegetation dynamics. The use of optical satellite time-series data increased the chances of acquiring cloud-free imagery [64,65]. Persistent cloud coverage can considerably affect the quality of optical S2 imagery [66]. S2 and S1 have varied revisit periods at 5 and 12 days, respectively, and these differences were accounted using seasonal median composites.
In general, the overall accuracies achieved were above 90% for all scenarios (combined and selected variables using GRRF). However, the higher accuracies of all combined scenarios (C1-C8) have been treated with caution because of the effects of multicollinearity that can lead to overfitting [58]. We thus placed more emphasis on the performance of the GRRF selected scenarios (C1selected-C8selected) which have the capabilities of reducing multidimensionality and the expected multicollinearity by maintaining the most relevant variables for the analysis. Therefore, the most important variables selected across the eight combined scenarios were S2 bands of NNIR, SWIR3, RE1, Red, Green and Blue; VI bands of NDWI and GNDVI; VP bands of Offset_value, TINDVIBeforeMax and TINDVIAfterMax; and S1 bands of VV and VH. The selection of the S2 bands could be explained by the narrowness of the NIR waveband region at 865 nm. The narrow NIR spectral region is known to be less contaminated by water vapor and represents the NIR plateau for vegetation while also being sensitive to some soil chemical properties [41]. Moreover, the reflectance in the SWIR3 band is sensitive to the leaf internal structure [68] whilst the red edge band (RE1), in addition to visible bands of red, green, and blue, are also sensitive to vegetation chlorophyll and other biochemical contents [69][70][71].
Furthermore, the VIs were also useful in specific crop type identification [72] by measuring the photosynthetic size of specific plant canopies that could improve the individual cropping patterns Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 26 October 2020 doi:10.20944/preprints202010.0517.v1 classification [73]. The selection of the VIs i.e., NDWI and GNDVI can be attributed to the sensitivity of the NIR band in the NDWI to leaf internal structure and leaf dry matter content while the SWIR3 band in the index is sensitive to the vegetation water content and the spongy mesophyll structure in vegetation canopies reflecting biochemical metrics of vegetation [47]. On the other hand, the GNDVI is more sensitive to the chlorophyll content of the plant since it is composed of the green channel instead of the red band [50]. In addition, time-series phenology variables are found to be the best in differentiating temporal and spectral variability of crop growth [74]. The selection of the VP variables of Offset_value, TINDVIBeforeMax, and TINDVIAfterMax could be explained by the sensitivity of Offset_value to land use and land cover differences, while TINDVIBeforeMax and TINDVIAfterMax describe the pre-and post-anthesis stages, respectively which can differ among different plants [53]. Regarding S1 radar backscatter data, the inter-channel phase information enabled by dualpolarization, i.e. VV and VH, allowed the enhanced analysis of the backscattering properties of the target areas of our study hence improved the classification [75]. Nonetheless, S1 data may be sensitive to soil moisture in the background of sparsely vegetated areas [76]. The synergetic importance of combining the S2 bands, S2 derived VIs, VP and S1 radar backscatter data was also demonstrated in this study in that the GRRF selected bands of scenario S2, VIs, VP, and S1 had the highest overall accuracy (OA=93.16%) compared to the rest of the GRRF selected band scenarios. The GRRF selected bands in S2, VIs, VP and S1 scenario were TINDVIAfterMax, NDWI, SWIR3, NNIR, RE1, Red, and Green. S2 bands' dominance confirmed the strength and importance of the raw spectral bands of S2 in discriminating vegetation [77].
We mapped cropping patterns that included monocrop avocado, mixed crop avocado, monocrop maize, mixed crop maize, monocrop pineapple, and monocrop coffee, grown in a very complex, heterogeneous, and dynamic agro-natural setup in Murang'a County, Kenya. In the field, we observed that most of the mixed fields were small-scale avocado and small-scale maize fields that are mixed with other crops like the common bean, banana, and macadamia. We speculate that the farmers intend to maximize the profit of their land by mixing the crops in the same piece of the land, particularly with the uncertainty in rainfall trend, and invasion of insect pests like the fall armyworm that could damage their maize crop [78]. High class-wise UA, PA, and F1-scores were observed in the mapping of monocropping patterns of avocado, coffee, tea, and pineapple. This could be due to the high spectral uniformity of the monocropping patterns, which also resulted in less intra class variability. However, mixed crop avocado and mixed crop maize classes had F1-scores below 90.00% compared to other cropping pattern classes in our study. This can be explained by high inter-class variability associated with the different vegetation composition including non-croplands such as natural vegetation, grasslands, and shrublands, hence a high rate of misclassification within the farms [23]. Nonetheless, the mixed crop avocado class had the greatest acreage (13,940.18 ha) compared to other crops, which could be attributed to the popularity of avocado farming by small-scale farmers in Murang'a County due to its increasing export value in Kenya. Furthermore, we assessed the performance of our classification experiments using an error matrix that takes into account the estimated area proportion of each class [61]. This is a robust and reliable classification accuracy assessment method compared to the traditional classification confusion matrix that considers the number of mapped instances of each class. Understanding the variability of cropping patterns and non-croplands can also contribute to the sustainable management of crop insect pests and pollinators that can be influenced by landscape structure [40].
Although non-croplands such as grassland, shrubland, natural vegetation, and water were not the focus of this present study, it is worth noting that relatively high UA, PA, and F1-scores were reported for mapping these classes across all (combined and GRRF) classification scenarios. The insignificant statistical differences from the McNemar's test between the combined scenarios and the GRRF selected variables could be explained by the use of related samples (i.e. the same training data used in the classification of all the scenarios) [63]. The use of the same training data in all the classification scenarios was necessary in this study to enable the unbiased comparison of the performance of the classification scenarios. Also, the insignificant differences among the performance of our classification scenarios indicate the novelty and robustness of our variable selection Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 26 October 2020 doi:10.20944/preprints202010.0517.v1 experiment using GRRF, in such a way that a few selected numbers of variables performance were not statistically different to all combined variables. The mapping of cropping patterns is of profound importance in understanding the sustainability of food systems and how they are affected by climate and also for modeling the abundance and spread of crop insect pests and disease and pollinators [19,79]. Furthermore, cropping patterns are valuable parameters in many crop modeling frameworks that estimate crop production as a function of many biotic and abiotic factors [80].

Conclusions
This study investigated the contribution of time-series S2 reflectance bands, its derivatives of vegetation indices and phenology metrics, and S1 radar backscatter data in the classification of cropping patterns in three sub-counties of Murang'a County, Kenya using eight scenarios and GRRF machine learning algorithm. The GRRF selected scenarios that consisted of S2 bands, their derivatives of vegetation indices and phenology metrics, and S1 radar backscatter data accurately (OA > 93.15%) mapped cropping pattern in a complex dynamic heterogeneous landscape. Our study also provided information on cropping pattern acreage that is necessary for crop management, and different other crop produce harvesting and marketing plans. Consequently, our approach can be used to map cropping patterns in other sites of similar agroecological conditions when high resolution imagery are not readily available.