PlanetScope imagery provided by NICFI
The PlanetScope constellation consists of multiple launches of groups of individual CubeSats (more than 130), and collectively delivers images at a spatial resolution of 3 m for the entire land/water surface of the Earth on a daily revisit time58. Through NICFI, we accessed the freely available high-resolution (4.77-m) satellite imagery mosaic, which was assembled from the original images acquired by the PlanetScope CubeSats. The preprocessed surface reflectance mosaics have been optimized with minimized effects of atmosphere and sensor characteristics by conducting calibration using Landsat data59,60. The mosaics cover the geographic area between 30 degrees latitude South and 30 degrees latitude North and were provided at two different cadences: a bi-annual product from December 2015 to August 2020 and a monthly product from September 2020 onwards. We visually inspected the cloud presence in the mosaics for the Sahel region. The bi-annual product is cloud-free with limited thin cirrus presence for either the periods of June–November or May–December. Image resolution blurring effects can occur inconsistently across the study area.
The mapping of Sahelian cropped and fallowed fields requires a cloud-free image acquired from a period of the second half of the growing season, approximately from mid-September26, when cropped fields are harvested, and separable from fallow40. Therefore, the PlanetScope bi-annual mosaic of June–November is ideal to use for this purpose. We downloaded the surface reflectance mosaics to cover the 200 mm to 600 mm isotype of the Sahel region (13N–17N, 25W–50E) for the year 2019. In addition to the four multispectral bands (blue, green, red, and near-infrared), the NDVI was generated as it is indispensable for separating cropped fields and fallow40. The resulting five-band image was passed as input data for the mapping of cropped and fallowed fields.
Mapping cropped and fallowed fields with NICFI and deep learning
Only recently, the mapping of Sahelian cropped and fallowed fields has been studied using the time series of Sentinel-240. Unlike traditional machine learning models employed by Tong et al., 202040, deep convolutional neural networks, which use convolution layers in addition to spectral information, enable the exploitation of contextual features in the image and thus outperform the traditional algorithms when the target has distinct shapes61–63. Cropped fields commonly have geometric straight-line boundaries created deliberately by human. These straight boundaries are the key textural feature to identify cropland. In the Sahel, in addition to the sectioned fields divided by lines into smaller segments, there are also large continuous fields without clear line segments. Our manually delineated training sets capture both types of fields. In addition, manured fields, excluded in previous studies because of their spectral inseparability from fallow, were included in the training data of this study. These fields have higher productivity and are located in close proximity of villages. Unarable lands, such as rocky terrain, hard pan outcrops, and rangelands, are incorporated into the training dataset as non-cropped and non-fallowed land, consequently, they are not mapped.
The manually mapped cropped fields were then used to train a supervised deep learning model, which was subsequently used to predict cropped fields for the entire study area. We also manually delineated fallowed fields and trained another deep learning model. Fallowed fields, characterized by herbaceous vegetation during fallow years and periodic weeding during cropped years, necessitated the use of multi-year images to distinguish them from savannas, which maintain herbaceous vegetation on consecutive years. Criteria for the visual identification of fallow were: (1) history of being actively cropped and (2) location adjacent to cropped fields in the same cultivation year, which separates fallow from the natural rangelands, which would by indistinguishable from fallow based on spectral properties alone. The mapped cropped fields are included as an additional input band to the fallow model, to enhance the capacity of the model to learn the adjacent features of cropped field and fallow, thus reducing the probability of unintendedly including savannah areas in the mapping. The delineation of cropped and fallowed fields excludes woody hedges occurring around fields in some of the Sahel regions, and the mapping was conducted within sandy soil (Extended Data Fig. 9) using the soil texture class39.
We employed the published U-Net architecture-based framework of satellite image segmentation49,64, after modifying it to accept multiple channels (five channels of R, G, B, NIR and NDVI for cropped field mapping, and six channels of R, G, B, NIR, NDVI and cropped field for fallow mapping). Tversky loss function was used to deal with class imbalance in segmentation. Alpha and beta were set to 0.4 and 0.6 and the patch size to 512 x 512 pixels to enable the model to understand a larger extent of imagery, which allows a better segmentation of continuous cropped field.
In total, we manually created 2107 (525 km2) cropped polygons and 895 (113 km2) fallow polygons within 389 training frames. We conducted an independent test dataset by randomly extracting 65 grid cells of 1x1 km and then randomly generating 5 points within each grid. For each point, we first assessed whether the point is within cropped fields in 2019 by visualizing the Surface Reflectance mosaics of 2019. We did not rely on an existing cropland mask for the current study to define/constrain cropland extent. Rather, the deep learning model also learns the difference between cropland and other land covers characterizing savannah ecosystems, such as bare land, grassland, shrub land etc. In the absence of an authoritative cropland mask, whether a point was within fallow, or part of the other, natural land cover types cannot be determined by a single-year image but requires its history over multiple years. To determine a fallow record, we visualized each of the annual Surface Reflectance mosaics over 4 years (2016–2019) —beyond which an uncultivated point would be defined as savannah. If a sample point was not cropped in 2019 but had been cropped in the previous years, we would attribute fallow to the sample point. The open-source Planet QGIS Plugin was used to facilitate multiple-year validation.
On the independent test data, we achieved a per-pixel user’s accuracy of 0.97 for actively cropped fields (Extended Data Table. 1), ensuring a highly reliable analysis when using it as a mask to extract information on soil quality and tree cover of cropped fields. As for the fallow areas, the user’s and producer’s accuracy are 0.71 and 0.75 respectively, largely due to instances of confusion with savannah in the agricultural frontiers, where small fields are newly cultivated and surrounded by large proportion of uncleared savannahs over the four consecutive years. The accuracy of mapped cropland (sum of actively cropped and fallowed fields) is 0.85, which is 10% higher than that of currently published maps40.
Spatial cultivated fraction
We define cultivated fraction as the ratio between cropped fields and the total area of cropped and fallowed fields for 1x1 km grid cells by conducting the following steps.
We first derived the percentage cover of cropped and fallowed field from 0 to 100% at a spatial resolution of 100 m using the fine resolution maps of cropped and fallowed field via upsampling. The cultivated fraction (CF) is defined as the ratio between cropped field cover to the total field cover of crops and fallows, multiplied by 100 (Fig. 1b). The Ruthenberg’s R-value65, which is a measure of land use intensity, is commonly applied to shifting cultivation systems by calculating the ratio between the length of active cultivation to the total length of active cultivation and fallow. However, our crop and fallow maps present a snapshot of 2019 and lack the temporal dimension necessary to allow an assessment of the changes in the fallow duration. We therefore adopted the ‘space-for-time’ approach by spatially aggregating the 100 m cultivated fraction map into 1 km of gridded information presenting the regional-level cultivated fraction to substitute temporal steps by cultivated fraction classes in space. When doing aggregation on the 100 m cultivated fraction map, we ran a sliding window of 1 km around each 100 m pixel and calculated the mean of cultivated proportion within the window to include the influence of neighboring situation to account for regional cultivated fraction (Fig. 1b). We defined five spatial cultivated fraction classes accordingly by equally dividing the range of cultivated fractions (Extended Data Table. 2). To maintain the focus on the agricultural environment, we excluded the grid cells with less than half the area covered by cropped and fallowed fields together, which are often located on the agricultural frontier and are plagued by a lower mapping accuracy due to a lower model transferability when applied to a very different landscape. The remaining cells are concentrated in spatially contiguous areas comprising historical agrarian cultures as shown in Fig. 3: the Serer’s in western Senegal, the Mandinka’s in western Mali, the Dogon’s in the Seno plain in Mali, the Zarma’s in western Niger, and the Hausa’s in central Niger.
To assess the applicability of the 'space-for-time' approach, we conducted visual comparisons spanning 33 years (1987–2019) in both an intensely cropped agroforestry parkland region and a marginal land area (Extended Data Fig. 10–11). We employed median Landsat imagery covering June-November with a 30-meter spatial resolution to map cropped and fallowed fields in Senegal for the year 1987. We adapted our 2019 training datasets to match the available Landsat imagery, being the only historical satellite-based record for the region. We developed two deep learning models—one for mapping cropped fields and another for identifying fallowed fields. It is important to note that the spatial and temporal resolutions of Landsat imagery may not yield pixel-level precision, but can effectively illustrate dominant land use patterns at a regional scale. In the intensively cropped parkland of the Serer region in western Senegal, fallowed fields diminished from 1987 to 2019, in particular the southeast, while the intensively cropped parkland in the west maintained its fallow extent (Extended Data Fig. 10). The marginal agricultural land adjacent to the savannas had longer fallow periods, long enough to allow fields to turn to savannas. This region is characterized by a low rate of cultivated land (Extended Data Fig. 11).
It's important to note that the 'space-for-time' approach does not imply consistent agricultural land use dynamics over the years, especially at the local level, where farming practices can vary between villages, farmers and over time. Instead, it highlights the enduring presence of signature intensified agricultural landscapes such as the dense parkland in western Senegal and the expansion of cropland in marginal areas, which have persisted for decades. Our methodology relies on these enduring features of cropping systems in the archetypical regions, and the cultivated fraction, crop growth rate and tree cover observed in recent years reflects the cumulative outcome of long-term farming practices.
Mapping crop growth rate using Sentinel-2
Ideally, we would utilize the daily PlanetScope time series for the rate of growth for cropped fields. However, daily PlanetScope data are not provided through the NICFI program and would have been very costly for our regional analysis. Instead, a Sentinel-2 5-day time series was used to extract the growth rate. We accessed the Sentinel-2 Level-2 Surface Reflectance dataset for the full calendar year 2019 via Google Earth Engine, which hosts the global Sentinel-2 L2 data downloaded from scihub and processed using sen2cor. The S2 cloud probability dataset that was created with the sentinel2-cloud-detector library was used as an initial step to remove clouds and cirrus. NDVI was further calculated for each image. In addition, the NDVI was used to remove any residual unmasked atmospheric artifacts where a decrease of NDVI greater than 0.02 between consecutive images indicated the presence of atmospheric disturbance. This was done iteratively for two times to couple with a series of two consecutive drops after one good observation. Then we reduced the raw NDVI time series into a 15-day composite time series by averaging to mitigate the spectral noise caused by scene-specific illumination conditions and residual cloud contaminations, and also to form a time series of regular time step. Lastly, gap filling was achieved by replacing missing values by the mean of the previous and next observation. The slopes of 2019 were calculated between NDVI values at adjacent time-steps with a unit of NDVI change per day. Consequently, one slope value was derived for every 15-day time-step. The maximum slope value of the 15-day slopes was used as the indicator of annual agricultural productivity for the cropped field.
d0 is the initial date, and t equals 15 in this study.
We investigated the correlation between millet yields measured in the field in 201954 and the corresponding crop growth rate (R-value = 0.75). It's important to note that the millet was harvested along transects designed in 1994, which may not perfectly align with cropped fields from different years due to rotational practices. Utilizing the cropped field of the analysis year to mask the crop growth rate proves beneficial for aligning the harvest data with the crop growth rate (Extended Data Fig. 12 b–d). In addition to crop growth rate, we calculated the annual maximum NDVI within the extent of cropped fields as an alternative proxy of annual agricultural productivity. The patterns of cultivated fraction across the five classes (very low to very high) remain consistent along the rainfall gradient when either crop growth rate or maximum NDVI is utilized (Extended Data Fig. 13).
Individual tree mapping using PlanetScope imagery
Traditional tree cover products do not detect trees on fields, which requires the mapping of trees as objects from very high-resolution images. Recently, a previous study49 has mapped all individual trees for the West African Sahel using sub-meter resolution images and deep learning. Here, we used a modified version using the same method but applied on PlanetScope images from 2019 and covering the entire African continent49,64. We used the raw NICFI PlanetScope level 2 images (rather than the NICFI basemaps) at a spatial resolution of 3 m. The difference in spatial resolution as compared to sub-meter resolution images limits the detection of very small trees, however, their contribution to the overall tree cover is marginal. The product maps scattered trees as individuals, allowing for detailed statistics on the cover of farmland trees. For the comparison with lower-resolution metrics such as cultivated fraction, crop growth rate and soil quality, we then aggregated the very high-resolution tree cover map to calculate accurate tree cover of cropped fields for the grid cells of 1x1 km.
Cropping system archetypes
We characterized cropping system archetypes using three metrics: cultivated fraction, crop growth rate, and tree cover of cropped fields, representing distinct clusters based on these metrics. To determine the optimal number of archetypes/clusters, we computed the silhouette score, a measure of how similar an object is to its own cluster compared to other clusters, for varying numbers, ranging from 2 to 11. The silhouette score, indicating the separation quality of archetypes/clusters, peaked at seven (Extended Data Fig. 14). Our analysis included six, seven and eight clusters, and we calculated the standard deviation of coverage percentages to ensure reliability.
Agro-ecological control metrics
We obtained the layers of soil texture class, soil organic carbon (SOC), soil pH (pH) and soil extractable Phosphorus (P), Potassium (K) and soil total organic Nitrogen (N) for 0–20 cm depth from the newly released 30 m map of African soil properties and nutrients39. We used the sandy soil class from the soil texture layer to restrict our study area within arable sandy soil extent, because the easy-to-cultivate sandy soil is preferable for agricultural use. The log-transformed soil parameters were converted back to their original unit using the provided formulas. The soil quality index (SQI) was used in our analysis to combine the soil chemical properties and nutrients, following Mukherjee and Lal (2014)55. The available soil layers made it possible to construct the SQI that focuses on nutrient storage capacity. The other two components of soil quality related to root development and water storage were not considered in our analysis. We constructed a principal component analysis based SQI in order to avoid a subjective definition of the weights. Principle components with eigenvalues greater than 0.5 were retained. Then the soil parameters with 10% highest non-highly correlated eigen vectors were chosen (Extended Data Fig. 15), which were SOC of the first component and pH of the second component. The SOC was min-max normalized to 0–1 and the pH was transformed to 0–1, where higher values represented closer to a pH of 7.2. The SQI was constructed by weighting the transformed SOC and pH by the variance of the corresponding component (Extended Data Table. 2).
Mean annual rainfall was derived using CHIRPS daily rainfall data from 1982–2020. The study area is defined as within rainfall isohyets 200–600 mm in the Sahel. Terrain elevation was from the Shuttle Radar Topography Mission (SRTM version 3) digital elevation data at a resolution of 1 arc-second (approximately 30m). Population density was from the Worldpop dataset for the year 2019. Larger cities and settlements were excluded by applying the total cropland area as a mask, such that it represents human population density in rural areas. We further aggregated the ~ 100 x 100m grid cells into 0.01° grid cells by summing the population per grid cell. Grid cells with > 50 persons km− 2 are defined as densely populated cell and the others as sparsely populated cells. We also used the arbitrary threshold of 25–75 persons km− 2 to test the sensitivity of the analysis on the associations between cultivated fraction, population size and crop growth rate (Extended Data Fig. 7–8).
A gradient boosting regression model was built to understand the associations between cultivated fraction and the selected agro-ecological control metrics, including rainfall, terrain elevation, soil quality, tree cover, population density, and crop growth rate. We used Grid Search to determine an optimal set of hyper-parameters for the gradient boosting model. Estimates of feature importance from the model were reported as relative importance of the control metrics to cultivated fraction.