Cropland extent maps for Africa

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:

  1. List the available cropland extent products

  2. Load the crop_mask_eastern product

  3. Plotting the different measurements of the crop-mask

  4. 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

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