Cropland extent maps for Africa

Products used: crop_mask_eastern, crop_mask_western, crop_mask_northern

Keywords: datasets; crop_mask_eastern, datasets; crop_mask_western, datasets; crop_mask_northern, 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 Eastern, Western, and Northern Africa show the estimated location of croplands in the countries for the period of January to Decemeber 2019:

crop_mask_eastern: Tanzania, Kenya, Uganda, Ethiopia, Rwanda, and Burundi
crop_mask_western: Nigeria, Benin, Togo, Ghana, Cote d’Ivoire, Liberia, Sierra Leone, Guinea, and Guinea-Bissau
crop_mask_northern: Morocco, Algeria, Tunisia, Libya, and Egypt

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 aconduct a simple analysis using the cropland extent product. The steps are as follows:

List the available cropland extent products
Load the crop_mask_eastern product
Plotting the different measurements of the crop-mask
Example application: Identifying crop trends with NDVI

For a more detailed example of using the cropland extent product, see the Cropland seasonal vegetation anomalies or Phenology_optical or Phenology_radar 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

/env/lib/python3.8/site-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
  warnings.warn(
/env/lib/python3.8/site-packages/datacube/storage/masking.py:7: DeprecationWarning: datacube.storage.masking has moved to datacube.utils.masking
  warnings.warn("datacube.storage.masking has moved to datacube.utils.masking",

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 interest
time_period: time period to load for the crop mask. Currently, only a map for 2019 is available
resolution: the pixel resolution to use for loading the crop_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]:

Make this Notebook Trusted to load map: File -> Trust Notebook

List cropland extent products available in Digital Earth Africa

We can use datacube’s list_measurements functionality to inspect the cropland extent products that are available. The table below shows the product name that we will 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_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_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, x: 773, y: 1011)
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 of 1 indicate the presence of crops, while a value of 0 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 to mask.
filtered: This band displays cropped regions as a binary map. Values of 1 indicate the presence of crops, while a value of 0 indicates the absence of cropping. This band is an object-based cropland extent map where the mask 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. The filtered dataset is provided as a complement to the mask 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>0).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();

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_18_0.png

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 some satellite data, calculate NDVI, and then 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'],
              like=cm.geobox
             )

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Applying pixel quality/cloud mask
Loading 145 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').median()

#plot the result
ds.plot.imshow(col='time', col_wrap=6, vmin=0, vmax=0.9);

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_24_0.png

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);

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_26_0.png

Summarise the trends in NDVI for the croppping regions with a line plot

Below, the plot shows that in this region, the cropping schedule starts in April, with the peak of the growing season occuring towards the end of the year around October/November.

[12]:

ds.mean(['x', 'y']).plot.line(marker='o', figsize=(10,5))
plt.title('Trend in NDVI for cropping region');

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_28_0.png

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:

[13]:

print(datacube.__version__)

1.8.5

Last Tested:

[14]:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

[14]:

'2021-09-29'