Cropland extent maps for Africa¶
Products used: crop_mask_eastern, crop_mask_western, crop_mask_northern, crop_mask_sahel, crop_mask_southern, crop_mask_southeast, crop_mask_central, crop_mask_indian_ocean
Keywords: datasets; crop_mask_eastern, datasets; crop_mask_western, datasets; crop_mask_northern, datasets; crop_mask_sahel, datasets; crop_mask_southern, datasets; crop_mask_southeast, datasets; crop_mask_central, datasets; crop_mask_indian_ocean analysis; agriculture, cropland extent
Background¶
A central focus for governing bodies in Africa is the need to secure the necessary food sources to support their populations. It has been estimated that the current production of crops will need to double by 2050 to meet future needs for food production. Higher level crop-based products that can assist with managing food insecurity, such as cropping watering intensities, crop types, or crop productivity, require as a starting point precise and accurate cropland extent maps indicating where cropland occurs. Current cropland extent maps are either inaccurate, have coarse spatial resolutions, or are not updated regularly. An accurate, high-resolution, and regularly updated cropland area map for the African continent is therefore recognised as a gap in the current crop monitoring services.
Digital Earth Africa’s cropland extent maps for Africa show the estimated location of croplands in the following countries for the period of January to Decemeber 2019:
crop_mask_eastern
: Tanzania, Kenya, Uganda, Ethiopia, Rwanda, and Burundicrop_mask_western
: Nigeria, Benin, Togo, Ghana, Cote d’Ivoire, Liberia, Sierra Leone, Guinea, and Guinea-Bissaucrop_mask_northern
: Morocco, Algeria, Tunisia, Libya, and Egyptcrop_mask_sahel
: Mauritania, Senegal, Gambia, Mali, Burkina Faso, Niger, Chad, Sudan, South Sudan, Somalia, Djibouti, and Eritreacrop_mask_southern
: South Africa, Namibia, Botswana, Lesotho, and Eswanticrop_mask_southeast
: Zimbabwe, Zambia, Mozambique, and Malawicrop_mask_central
: Angola, Democratic Republic of the Congo, Congo, Gabon, Cameroon, Equatorial Guinea, and Central African Republiccrop_mask_indian_ocean
: Madagascar, Mauritius, Reunion, and Comoros
For a full description of the product specifications, validation results, and methods used to develop the products, see the Cropland_extent_specifications document.
Description¶
This notebook will show you how to load, plot, and conduct a simple analysis using the cropland extent product. The steps are as follows:
List the available cropland extent products
Load the
crop_mask_eastern
productPlotting the different measurements of the crop-mask
Example analysis 1: Identifying crop trends with NDVI
Example analysis 2: Comparison of cropped area with global land cover datasets
For a more detailed example of using the cropland extent product, see the following notebooks:
Getting started¶
To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.
Load packages¶
Import Python packages that are used for the analysis.
[1]:
%matplotlib inline
import datacube
import xarray as xr
import matplotlib.pyplot as plt
from deafrica_tools.plotting import display_map
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.datahandling import load_ard
Connect to the datacube¶
Connect to the datacube so we can access DE Africa data.
[2]:
dc = datacube.Datacube(app='cropland_extent')
Analysis parameters¶
This section defines the analysis parameters, including:
lat, lon, buffer
: center lat/lon and analysis window size for the area of interesttime_period
: time period to load for the crop mask. Currently, only a map for 2019 is availableresolution
: the pixel resolution to use for loading thecrop_mask_<region>
. The native resolution of the product is 10 metres i.e.(-10,10)
The default location is in a extensivley cultivated valley north of Addis Ababa, Ethiopia
[3]:
lat, lon = 8.5615, 40.691
buffer = 0.04
time_period = ('2019')
resolution=(-10, 10)
#join lat, lon, buffer to get bounding box
lon_range = (lon - buffer, lon + buffer)
lat_range = (lat + buffer, lat - buffer)
View the selected location¶
The next cell will display the selected area on an interactive map. Feel free to zoom in and out to get a better understanding of the area you’ll be analysing. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point.
[4]:
display_map(lon_range, lat_range)
[4]:
List cropland extent products available in Digital Earth Africa¶
We can use datacube’s list_measurements
functionality to inspect the cropland extent products available in the datacube. The table below shows the product names you can use to load data, and the measurements available for each product. The cropland extent masks come with three measurements: mask
, filtered
, and prob
.
[5]:
dc_measurements = dc.list_measurements()
dc_measurements.filter(like='crop_mask', axis=0)
[5]:
name | dtype | units | nodata | aliases | flags_definition | ||
---|---|---|---|---|---|---|---|
product | measurement | ||||||
crop_mask_central | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_eastern | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_indian_ocean | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_northern | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_sahel | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_southeast | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_southern | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN | |
crop_mask_western | mask | mask | uint8 | 1 | 0.0 | [crop_mask, MASK] | NaN |
prob | prob | uint8 | 1 | 0.0 | [crop_prob, PROB] | NaN | |
filtered | filtered | uint8 | 1 | 0.0 | [mode] | NaN |
Loading the cropland extent product¶
In this example, we will load the 'crop_mask_eastern'
product.
[6]:
# generate a query object from the analysis parameters
query = {
'time': time_period,
'x': lon_range,
'y': lat_range,
'resolution':resolution
}
# now load the crop-mask using the query
cm = dc.load(product='crop_mask_eastern',
**query)
print(cm)
<xarray.Dataset>
Dimensions: (time: 1, y: 1011, x: 773)
Coordinates:
* time (time) datetime64[ns] 2019-07-02T11:59:59.999999
* y (y) float64 1.093e+06 1.093e+06 ... 1.083e+06 1.083e+06
* x (x) float64 3.922e+06 3.922e+06 3.922e+06 ... 3.93e+06 3.93e+06
spatial_ref int32 6933
Data variables:
mask (time, y, x) uint8 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
prob (time, y, x) uint8 41 42 41 45 43 45 40 ... 17 18 24 21 23 22
filtered (time, y, x) uint8 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
Attributes:
crs: EPSG:6933
grid_mapping: spatial_ref
Plotting the cropland extent¶
Above, we can see the crop_mask_<region>
products contains three measurements:
mask
: This band displays cropped regions as a binary map. Values of1
indicate the presence of crops, while a value of0
indicates the absence of cropping. This band is a pixel-based cropland extent map, meaning the map displays the raw output of the pixel-based Random Forest classification.prob
: This band displays the prediction probabilities for the ‘crop’ class. As this service uses a random forest classifier, the prediction probabilities refer to the percentage of trees that voted for the random forest classification. For example, if the model had 200 decision trees in the random forest, and 150 of the trees voted ‘crop’, the prediction probability is 150 / 200 x 100 = 75 %. Thresholding this band at > 50 % will produce a map identical tomask
.filtered
: This band displays cropped regions as a binary map. Values of1
indicate the presence of crops, while a value of0
indicates the absence of cropping. This band is an object-based cropland extent map where themask
band has been filtered using an image segmentation algorithm (see this paper for details on the algorithm used). During this process, segments smaller than 1 Ha (100 10m x 10m pixels) are merged with neighbouring segments, resulting in a map where the smallest classified region is 1 Ha in size. Thefiltered
dataset is provided as a complement to themask
band; small commission errors are removed by object-based filtering, and the ‘salt and pepper’ effect typical of classifying pixels is diminished.
Below, we will plot the three measurements side-by-side:
[7]:
fig, axes = plt.subplots(1, 3, figsize=(24, 8))
cm.mask.plot(ax=axes[0],
cmap='Greens',
add_labels=False)
cm.filtered.plot(ax=axes[1],
cmap='Blues',
add_labels=False)
cm.prob.where(cm.prob<255).plot(ax=axes[2],
cmap='magma',
add_labels=False)
axes[0].set_title('"Mask": pixel-based cropland extent')
axes[1].set_title('"Filtered": object-based cropland extent')
axes[2].set_title('"Prob": Probabilities of cropland');
plt.tight_layout();

Example application: Identifying crop trends with NDVI¶
The Normalised Difference Vegetation Index (NDVI) can be used to track the health and growth of crops as they mature. Here will load an annual time series of Sentinel-2 satellite data, calculate the NDVI, and mask the region with DE Africa’s cropland extent product. With only the cropping pixels remaining, we can summarise the growth cycle of the crops within the region.
[8]:
ds = load_ard(dc=dc,
products=['s2_l2a'],
time=('2019-01', '2019-12'),
measurements=['red', 'nir'],
mask_filters=(['opening',5],['closing', 5]), #improve cloud mask
min_gooddata=0.35,
like=cm.geobox
)
Using pixel quality parameters for Sentinel 2
Finding datasets
s2_l2a
Counting good quality pixels for each time step
Filtering to 114 out of 145 time steps with at least 35.0% good quality pixels
Applying morphological filters to pq mask (['opening', 5], ['closing', 5])
Applying pixel quality/cloud mask
Loading 114 time steps
Calculate NDVI¶
[9]:
ds = calculate_indices(ds, 'NDVI', drop=True, collection='s2')
Dropping bands ['red', 'nir']
Resample the datset to monthly time-steps¶
[10]:
ds = ds.NDVI.resample(time='MS').mean()
#plot the result
ds.plot.imshow(col='time', col_wrap=6, vmin=0, vmax=0.9);

Mask with the cropland extent map¶
Using the measurement filtered
, we can mask the dataset with the crop_mask_eastern
product.
[11]:
#Mask using 'filtered' band
ds = ds.where(cm.squeeze().filtered)
ds.plot.imshow(col='time', col_wrap=6, vmin=0, vmax=0.9);

Summarise the trends in NDVI for the croppping regions with a line plot¶
Below, the plot shows the cropping schedule starts in April, with the peak of the growing season occuring towards the end of the year in October/November.
[12]:
ds.mean(['x', 'y']).plot.line(marker='o', figsize=(10,5))
plt.title('Trend in NDVI for cropping region');

Example Analysis 2: Comparison with global landcover datasets¶
Indexed within the DE Africa’s datacube are two global 10m resolution landcover datasets that each contain a cropping class. These landcover datasets were both produced using Senitnel-2 images, the same as DE Afric’s cropland extent product. Below, we will load the landcover datasets over the same region as the cropland extent map and compare the area classified as croppping. You can use the interactive map plotted below to determine the areas where each datasets does well (or poorly) at classifying crops.
To learn more about the landcover datasets see the Landcover Classification notebook.
First, we need to re-enter some analysis parameters. You’ll notice below we’ve set the resolution to 60m to allow us to load a large region quickly, though all datasets are natively stored at 10m resolution.
The default region is a bounding box over Lake Kyoga, Uganda. This region is extensively cropped.
[13]:
lat, lon = 1.4315, 33.1207
buffer = 1.0
resolution=(-60, 60) #resample to load larger area
#join lat, lon, buffer to get bounding box
lon_range = (lon - buffer, lon + buffer)
lat_range = (lat + buffer, lat - buffer)
Plot an interactive map¶
If you zoom in you can examine the dominant land cover classes in the region (hint, there’s a lot of agriculture!)
[14]:
display_map(lon_range, lat_range)
[14]:
Load the cropland extent, ESRI’s Landover product, and ESA’s WorldCover product¶
[15]:
# generate a query object from the analysis parameters
query = {
'x': lon_range,
'y': lat_range,
'resolution':resolution
}
# now load the crop-mask using the query
cm = dc.load(product='crop_mask_eastern',measurements=['mask'], **query, time='2019').squeeze()
#load the two landcover datasets
lulc_esa = dc.load(product='esa_worldcover', measurements='classification', like=cm.geobox).squeeze()
lulc_esri = dc.load(product='io_lulc', measurements='classification', like=cm.geobox).squeeze()
Reclassify the landcover datasets to a binary crop/non-crop image to allow a straightforward comparison with DE Africa’s cropland extent map
[16]:
esri_crops = xr.where(lulc_esri['classification']==5, 1, 0)
esa_crops = xr.where(lulc_esa['classification']==40, 1, 0)
Plot the cropping extent of the three datasets¶
[17]:
fig,ax=plt.subplots(1,3, sharey=True, figsize=(21,7))
cm['mask'].plot.imshow(ax=ax[0], add_colorbar=False)
esri_crops.plot.imshow(ax=ax[1], add_colorbar=False)
esa_crops.plot.imshow(ax=ax[2], add_colorbar=False)
ax[0].set_title('DE Africa Cropland Extent')
ax[1].set_title('ESRI Land Cover, cropland class')
ax[2].set_title('ESA Worldcover, cropland class')
plt.tight_layout();

Calculate the cropped area in each product¶
In this example, you’ll see that DE Africa’s cropland product maps a lot more cropping than the two global landcover products. In this example, the DE Africa product is more accurate than the landcover datasets. However, over wetlands in the south-west of Lake Kyoga DE Afria’s product incorrectly maps some of the wetlands as cropping, while the landcover datasets do not. While there are regions where DE AFrica’s cropland extent maps have errors, in general they are much more accurate than any of the currently existing landcover datasets over Africa.
[18]:
pixel_length = query["resolution"][1]
m_per_km = 1000 # conversion from metres to kilometres
area_per_pixel = pixel_length**2 / m_per_km**2
[19]:
cm_area = cm.mask.sum().data * area_per_pixel
esri_area = esri_crops.sum().data * area_per_pixel
esa_area = esa_crops.sum().data * area_per_pixel
label = ['DEAfr crop mask', 'ESRI crops', 'ESA crops']
plt.bar(label, [cm_area,esri_area,esa_area])
plt.ylabel("Area (Sq. km)")
plt.title('Cropped Area Comparison');

Additional information¶
License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Africa data is licensed under the Creative Commons by Attribution 4.0 license.
Contact: If you need assistance, please post a question on the Open Data Cube Slack channel or on the GIS Stack Exchange using the open-data-cube
tag (you can view previously asked questions here). If you would like to repoart an issue with this notebook, you can file one on
Github.
Compatible datacube version:
[20]:
print(datacube.__version__)
1.8.6
Last Tested:
[21]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[21]:
'2022-03-09'