Climate variability—combined with a lack of resources, social and political instability, pest outbreaks, and other contributing factors—have led to food-insecurity events throughout sub-Saharan Africa, compromising the lives and livelihoods of the most vulnerable populations (Devereux, 2009; Samasse et al., 2018; Funk, 2021). Homegrown food production is a function of crop area and crop yield, but these components are difficult to assess because agricultural statistics in sub-Saharan Africa are known to be inaccurate due to poor organization and data analysis (Devereux, 2009; Carletto et al., 2015). They are also quite coarse, being reported at the administrative unit level for ground-based statistics, and generally, 300–1,000 m pixels for satellite-based maps (Carletto et al., 2015; Samasse et al., 2018). The recent availability of 10 m Sentinel-2 data in Google Earth Engine (GEE) allows for efficient processing of high spatial resolution data, making high spatial resolution crop area maps over large areas feasible (Chivasa et al., 2017; Samasse et al., 2018; Jin et al., 2019; Amani et al., 2020; Karlson et al., 2020; Kerner et al., 2020; Masiza et al., 2020; Tseng et al., 2020; Verde et al., 2020).

The high temporal (5-days repeat) resolution of Sentinel-2 data allows for the improved observation—and differentiation—of crop and natural vegetation phenology, as well as a higher likelihood of minimizing cloud impacts on the time series (Misra et al., 2020). The high spatial resolution allows for fewer mixed pixels in these landscapes characterized by smallholder farms, which result in mosaics of fields, forests, and pastures, often heterogeneously mixed at even the 30 m Landsat scale, and certainly at the 500 m MODIS scale (Ozdogan and Woodcock, 2006; Samasse et al., 2018; Jin et al., 2019; Misra et al., 2020). The major drawback of small pixel sizes is huge data set size, complicating both data storage/transfer and computational requirements. The use of GEE reduces these requirements because: (1) the data sets are already loaded into GEE; they do not have to be ordered, downloaded, and stored locally, and (2) processing can be done on Google's server cloud, effectively bringing supercomputing to the average user, for free (Gorelick et al., 2017).

Remotely sensed data have been used in agricultural applications since satellites were launched in the early 1970s (Hammond, 1975; MacDonald and Hall, 1980). Opening the Landsat archive to free access in 2008, quasi-daily MODIS data, and high temporal resolution Sentinel-2 data have allowed for techniques that leverage phenology to map crops and produce yield estimates (Lobell, 2013; Wang et al., 2020). The phenology of crops may differ from that of natural areas in many ways, due to growth form, irrigation, harvest, field management, and other factors. As the season progresses, crop areas will show bare soil due to plowing or clearing to prepare the area, followed by a steady increase in green vegetation as crops grow, then some vegetation die-back as crops mature, particularly in seasonal rainfed agriculture areas. Finally, there is a rapid decrease in vegetation amount due to harvest, leading to bare fields with a mix of soil and crop residue. In contrast, natural areas may show different trends. For instance, forested areas may stay green year-round; shrublands may green up earlier than crops or grasslands and stay green longer due to established root systems; grasslands may show similar timing in green up to crops in rainfed agricultural areas, but dry down would be more drawn out due to a lack of harvest, and bare soil likely would not be exposed.

Most of the studies using crop phenology to map crop area use either raw bands, vegetation indices (VIs) such as the Normalized Difference Vegetation Index (NDVI, Rouse et al., 1973), or a combination of the two (Samasse et al., 2018; Jin et al., 2019; Amani et al., 2020; Karlson et al., 2020; Kerner et al., 2020; Masiza et al., 2020; Tseng et al., 2020; Verde et al., 2020). In this study we took a different approach. Spectral mixture analysis (SMA, Roberts et al., 2002) decomposes the signal of a pixel into percentages/fractions of the scene components making up the pixel. These scene components are termed endmembers (EMs), and generally consist of the spectrally distinct constituents: green vegetation (GV), non-photosynthetic vegetation (NPV), soil, and shade. SMA works well on data sets that include broad spectral coverage, such as Landsat or MODIS. The addition of bands on the red edge for Sentinel-2 in comparison to earlier sensors further increases the confidence that SMA produces physically meaningful EM fractions. GV is generally highly correlated with NDVI and other VIs focusing on vegetation greenness [e.g., enhanced vegetation index (EVI, Huete et al., 1997)]. The other three EM fractions provide unique information. This is important, for instance, because a drop in NDVI in a pixel can result from either browning of vegetation (leading to a mixture of GV and NPV) or a partial crop harvest (leading to a mixture of GV and newly exposed soil), and while NDVI cannot distinguish between the two events, SMA can.

Malawi normally receives enough precipitation to produce most of the maize and other crops required to feed its people; there has not been a famine since 2001–2002 (Devereux 2009). However, what makes Malawi an excellent case study are two existing national reference data sets of crop area: (1) district-wide statistics from the Malawi Department of Agriculture, and (2) from a map based on manual interpretation of landcover on a 1 km grid by Gray Tappan of the U.S. Geological Survey (USGS). We also compare our results to two global data sets, from IIASA-IFPRI (Fritz et al., 2015) and a protype map from the European Space Agency (ESA, ESA-CCI, 2021).

In this study, we used 2018 Sentinel-2 data to map crop area in central and southern Malawi (6752347 ha) at a 10 m pixel resolution. The basic approach was to manually identify reference data points using very high spatial resolution imagery (<1 m), randomly divide that data set into training and testing, extract phenologically-based predictor variables from Sentinel-2 for the training data set, input those predictor variables into random forests (RF) to map crop area, and evaluate classification accuracy, both with the test data set, and in relation to crop area from the two Malawi reference data sets and the two global data sets.

The EM fraction trajectories reflect changes due to plant phenology. To highlight the separability of the different cover types, Figure 4 shows sample EM trajectories for a crop field and a natural shrubland area, located within 100 m of each other. The GV EM trajectory for the natural area follows precipitation well; green up begins in October, peak greenness occurs between January and March except for 2017, when the rainfall peak was a month late (Figure 1) and the GV peak was in April, and then there is a steady decline as the landscape dries down during the 6 months long dry season (Figure 4). Crop greenup is delayed a few months compared to shrub greenup, the timing of peak greenness is later in the season, and the peak is more pronounced than that of the natural area. Finally, there is an extended period of near zero GV values between harvest and planting the following growing season. Negative EM values, while not physically meaningful, occur because SMA is a matrix algebra transform with the constraint that fractions sum to 1.0, and can be interpreted as the absence of that EM for that time step. For this particular example the natural area demonstrated a higher GV fraction than the crop fields throughout 2018, but that is not always the case. NPV for the natural area has an inverse trend to GV, the new green leaves are dominant at the beginning of the season so NPV is low, then the leaves senesce, and woody shrub material and senesced ground cover become more exposed to the satellite, leading to an increase in NPV. A similar pattern is present during the growing season for the crop field, although it appears that harvest residue was left on the field in 2017 and 2018 as NPV remains high during the initial increase in GV for those years. The soil fractions stayed relatively consistent for the two areas in Figure 4, showing no evidence of plowing in this particular area, which is consistent with the observation, based on NPV, that harvest residue was left on the soil. Soil fraction was higher for crop areas, likely due to the presence of some exposed soil between crop rows. NPV and soil show more noise than GV, likely due to residual cloud contamination. Shade shows a similar pattern to GV as the landscape absorbs more light (high shade) when it is highly vegetated, and reflects more (low shade) when the ground is bare. Some locations experienced persistent cloud cover during parts of the rainy season, resulting in missing data in the satellite analysis and corresponding gaps in the plots shown in Figure 4.

We only present output for the monthly predictor variables (the second method of data distillation), as accuracies were slightly better, and the predictor variable data set was half the size of the bimonthly data set, thus it was more convenient. The overall accuracies when compared with independent test data for the six tiles with the most training data (Table 1) were on the order of 85% (Table 3). We also show data for classes re-coded to crop/no crop (Table 4), which is the more common way that crop area classifications are presented. Accuracies increased 5% for five of the six tiles, with three tiles exhibiting accuracy over 90%. Further analysis of the results in Tables 3, 4 show that the user's accuracy for crop tends to be a bit higher than the producer's accuracy, indicating a slight under-identification of agriculture in the model results, and an expected underestimate of the overall cropped area by a small amount. The biggest source of confusion was between crops and shrubs, though more crops were classified as shrub than vice versa (Table 3). It may be that the diversity in crop types and or farming intensity is more variable than shrub types, so the crop class is more heterogeneous and thus members of the class are more likely to appear to behave like shrubs.

Table 5 presents the variable importance values for each of the predictor variables for each of the six tiles in Table 1, with the top ten variables for each tile in bold. There are a number of interesting features. GV variables were generally most important. The median GV at the beginning of the growing season ranked high, while the median GV for the peak of the growing season (February) showed much lower importance. Natural areas green up earlier than crop areas, with separability being reduced at peak greenness (Figure 4). This also points to the advantage of monthly vs. bi-monthly variables as January and February values would have been grouped together, reducing the signal. Natural areas also remain green longer due to enhanced water availability due to established root systems, and this is reflected in high importance scores for GV variables as the season progresses. NPV showed lower importance than GV, though January was in the top ten for three tiles. January soil was also highly important. Also, soil in May, June, and July, while only in top ten importance for the NWW tile, tended to have relatively high importance. Hence, soil was important early in the growing season and during the harvest season when soil may be exposed in crop areas while natural vegetation is largely comprised of GV and NPV (Figure 4). Shade variables were second to GV in importance, with July and August being the most important months for four of the tiles. At this time of the year natural vegetation would show canopy shading due to uneven vegetation heights both within individual plants and between neighboring plants while harvested crop area would be relatively flat with higher albedo. For the remaining two tiles January was important, for similar reasons—young crops would be too small to cast much in the way of shadows compared to larger natural vegetation. The annual summary variables were generally not of high importance.

We mosaicked the maps from the individual tiles to generate a crop area map for all of southern and central Malawi. Seam lines were not noticeable, likely because the classification models were accurate, so there were not differences across tile boundaries. Table 6 presents crop area estimates for the different data sources. Agreement in percent crop area at the district level (comparing sets of 22 values) is high for RF, Malawi government, and Tappan USGS data sets. R² for fractional crop area at the district level between RF and Malawi government is 0.74, between RF and USGS it is 0.79, and between Malawi government and USGS it is 0.53 (Table 6). It is interesting that the relationships between RF and Malawi government and RF and USGS are similar, but the relationship between Malawi government and USGS is lower. The primary differences between the RF map and Malawi government are possible double-counted fields and human error in gathering data, both data sources cover the entire region. Both RF and USGS rely on manual image interpretation to label landcover, but RF uses all of the imagery whereas USGS examined imagery on a 1 km grid. The Malawi government data and USGS map have less in common. The ESA-CCI map shows good agreement with the RF Sentinel-2 map (R² of 0.63), however Figure 5 suggests that the ESA map systematically overpredicts crop area in southern Malawi.

Our accuracies compare favorably with the three published studies examining crop area in sub-Saharan African countries at 10 m resolution that we are aware of Kerner et al. (2020), Tseng et al. (2020), and Jin et al. (2019). Kerner et al. (2020) obtained 83% accuracy when classifying crop/no crop in the country of Togo. Jin et al. (2019) obtained 85% accuracy for crop/no crop in Kenya and Tanzania. Tseng et al. obtained 86% accuracy in Kenya. However, our accuracy values are based on many more samples, making them more robust.

Aside from the increase in accuracy, the main benefit of utilizing SMA for crop area mapping is that you use the entire electromagnetic spectrum as measured by Sentinel-2, reprojected into a more intuitive/interpretable GV, NPV, soil, shade data space. Jin et al. (2019) also used a suite of greenness measures for prediction, as well as raw bands, but used additional information on vegetation structure from radar data from Sentinel-1. NPV and particularly shade can give information on vegetation structure (Roberts et al., 2002). We explored the use of Sentinel-1 data but ultimately decided against it due to some irregularities in image registration between data from the two sensors, and because of the presence of shadows/data gaps in the radar data behind taller features, due to radar being side-looking.

We chose to use TOA data rather than surface reflectance because it allows us to do temporal analyses of crop area over 5 years (2016–2020) rather than three (2018–2020). For instance, 2016 and 2017 were low rainfall years in Kenya and we could study the effect of the drought on crop area. TOA data contain atmospheric scattering that primarily affects the visible bands. For instance, for surface reflectance data, the blue band (490 nm) reflectance is usually similar to red reflectance (665 nm) for GV, whereas for TOA data there is additional blue light (Table 2). Spatio-temporal variability in atmospheric scattering should incur noise when using TOA vs. surface reflectance data, however as atmospheric scattering primarily only affects three of the ten bands it may not have had a strong effect on EM fractions, and our classification accuracies remained high. Other sources of noise include residual cloud contamination and occasional geolocation errors (particularly considering the heterogeneity of small-scale farming areas). For instance, Tremas et al. (2015) found geolocation errors for Sentinel-2 data of 12.5 m, or over 1 pixel, which could cause an agricultural area pixel to have the EM values of a neighboring forest pixel at certain timesteps, affecting monthly median and annual summary values.

Fritz et al. (2015) state that the global overall accuracy of the IIASA-IFPRI map is 82.4%, but it appears to be less accurate in Malawi, as the range in values is less than the other maps, and agreement at the district scale was low. Our finding that the IIASA-IFPRI map differed from the three other reference data sets and our classification map is in line with the findings of Samasse et al. (2018). They compared eight landcover maps in West African countries, and found that the coarser resolution maps, including IIASA-IFPRI, performed much worse than 30 m maps using Landsat data. Interestingly, the ESA-CCI 20 m map, using Sentinel-2 data, also performed worse than the Landsat-based maps with accuracies 10–40% points lower. Samasse et al. (2018) suggested the techniques and West Africa training data used in the ESA-CCI product needed to be re-examined. Alkhalil et al. (2020) also evaluated the ESA-CCI map in West Africa. They found extremely low producer's accuracies (over-prediction of crops) for three polygons in the Sahel (0.07, 0.34, 0.03) and less, but still some over-prediction for three polygons closer to the coast (producer's accuracies of 0.61, 0.72, 0.56), leading them to declare the ESA-CCI product was not an acceptable crop mask. Our research also showed overprediction of crop area by ESA-CCI when compared with the RF Sentinel-2 map, particularly in southern Malawi (Table 6, Figure 5).

RF are known as a greedy classifier, hence, a large number of manually interpreted points were identified—a total of 20,848. For comparison to other studies, central and southern Malawi comprise 6,752,347 ha. Kerner et al. (2020) mapped percent crop area in Togo (5,678,500 ha) using a different algorithm and hence, a fraction of the reference data-−1,319 crop or no crop points from within the country, and a global data set of 35,866 crop/no crop points. Tseng et al. (2020) used a similar algorithm and the same global 35,866 points, as well as 14,080 crop/no crop points for the much larger country of Kenya (58,036,700 ha). Jin et al. (2019) used RF to map crop/no crop in both Kenya and Tanzania (94,730,300 ha), 4,509 points were identified in Kenya, and 4,140 in Tanzania. The country areas are much larger than that of Malawi; however, the agricultural portion of the countries is much smaller than these numbers, making the comparison less straightforward, because large parts of each country could be excluded as potentially arable. In fact, latitude and longitude were included as predictor variables in Jin et al. (2019). Furthermore, it is not known exactly how many fewer points would have been needed if our analysis had been crop/no crop instead of crop, sparse, shrub, tree, and water. Another factor is we used roughly 50% of the data for training, 50% for testing while Jin et al. (2019) and Kerner et al. (2020) trained on 80% of the data and Tseng et al. (2020) trained on 90% of the data. Hence, they leveraged their points differently, and had fewer points for robust accuracy assessment. While our method may be less efficient and/or an excessive amount of reference points may have been identified, our accuracies are based on far more data points, which increases confidence.

This study generated a crop area map for southern and central Malawi with very high crop/no crop classification accuracies between 87 and 93%. The approach combined some of the oldest (air photo interpretation) and newest (cloud computing) remote sensing techniques. Overall, we find the results presented here to be a promising result for the potential to use GEE and RF algorithms to produce high-quality cropped area estimates for smallholder farms in rainfed agriculture regimes. EM phenology predictor variables and variable importance measures combine to produce interpretable, non-“black box” results. Accuracies for the two most populous classes, crop and shrub, are based on 50% of the reference data points so are extremely robust. A more accurate assessment of the cropped area in a region can help better understand the dynamics impacting food production and food security in vulnerable regions of sub-Saharan Africa. The method can be repeated using the same training data to perform the classification over multiple years which would be useful for the identification of changing landscape dynamics and for seasons when weather, agricultural inputs (seeds or fertilizer), or labor impacted the amount of cropped area for a particular season, over a region, relatively simply. The results presented here can also be utilized to develop an agricultural mask that can be used to better focus agroclimatic analysis over wide areas. All these different benefits could have a positive impact on the ability to anticipate, assess, and mitigate the impacts of cropped area on food security over sub-Saharan Africa.

Publicly available datasets were analyzed in this study. This data can be found here: Publicly available Sentinel-2 data were analyzed in this study, this data can be found in Google Earth Engine.

SP and GH conceptualized the research, revised, and improved the manuscript. SP defined the methods, analyzed the findings, composed the figures, and drafted the manuscript. Both authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.