Cartes de l’étendue des terres cultivées en Afrique

Produits utilisés : crop_mask

Mots-clés : datasets; crop_mask, analysis; agriculture, cropland extent

Aperçu

Les organes directeurs africains se concentrent sur la nécessité de garantir les sources de nourriture nécessaires pour subvenir aux besoins de leurs populations. On estime que la production actuelle de cultures devra doubler d’ici 2050 pour répondre aux besoins futurs en matière de production alimentaire. Les produits de plus haut niveau basés sur les cultures qui peuvent aider à gérer l’insécurité alimentaire, tels que l’intensité d’arrosage des cultures, les types de cultures ou la productivité des cultures, nécessitent comme point de départ des cartes précises et exactes de l’étendue des terres cultivées indiquant où se trouvent les terres cultivées. Les cartes actuelles de l’étendue des terres cultivées sont soit inexactes, ont des résolutions spatiales grossières ou ne sont pas mises à jour régulièrement. Une carte précise, à haute résolution et régulièrement mise à jour des superficies des terres cultivées pour le continent africain est donc reconnue comme une lacune dans les services actuels de surveillance des cultures.

Les cartes de l’étendue des terres cultivées de Digital Earth Africa pour l’Afrique montrent l’emplacement estimé des terres cultivées pour la période de janvier à décembre 2019.

Pour une description complète des spécifications du produit, des résultats de validation et des méthodes utilisées pour développer les produits, consultez le document Cropland_extent_specifications.

Description

Ce bloc-notes vous montrera comment charger, tracer et effectuer une analyse simple à l’aide du produit d’étendue des terres cultivées. Les étapes sont les suivantes :

Listez les produits d’étendue des terres cultivées disponibles
Charger le produit « crop_mask »
Tracé des différentes mesures du crop-mask
Exemple d’analyse 1 : Identifier les tendances des cultures avec le NDVI
Exemple d’analyse 2 : Comparaison des superficies cultivées avec les ensembles de données sur la couverture terrestre mondiale
Inspectez différentes régions du masque de recadrage.

Pour un exemple plus détaillé de l’utilisation du produit d’étendue des terres cultivées, consultez les blocs-notes suivants :

Commencer

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Charger des paquets

Importez les packages Python utilisés pour l’analyse.

[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

Se connecter au datacube

Connectez-vous au datacube pour que nous puissions accéder aux données de DE Africa.

[2]:

dc = datacube.Datacube(app='cropland_extent')

Paramètres d’analyse

Cette section définit les paramètres d’analyse, notamment :

« lat, lon, buffer » : latitude/longitude centrales et taille de la fenêtre d’analyse pour la zone d’intérêt
time_period : période de temps à charger pour le masque de recadrage. Actuellement, seule une carte pour 2019 est disponible
resolution: la résolution en pixels à utiliser pour charger le crop_mask_<region>. La résolution native du produit est de 10 mètres, soit (-10,10)

L’emplacement par défaut se trouve dans une vallée largement cultivée au nord d’Addis-Abeba, en Éthiopie

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

Afficher l’emplacement sélectionné

La cellule suivante affichera la zone sélectionnée sur une carte interactive. N’hésitez pas à zoomer et dézoomer pour mieux comprendre la zone que vous allez analyser. En cliquant sur n’importe quel point de la carte, vous découvrirez les coordonnées de latitude et de longitude de ce point.

[4]:

display_map(lon_range, lat_range)

[4]:

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

Liste des produits d’étendue des terres cultivées disponibles dans Digital Earth Africa

Nous pouvons utiliser la fonctionnalité « list_measurements » de Datacube pour inspecter les produits d’étendue des terres cultivées disponibles dans le Datacube. Le tableau ci-dessous indique les noms de produits que vous pouvez utiliser pour charger des données et les mesures disponibles pour chaque produit. Les masques d’étendue des terres cultivées sont fournis avec trois mesures : « mask », « filtered » et « prob ».

Nous pouvons voir dans le tableau ci-dessous qu’il existe un seul produit « crop_mask » qui couvre l’ensemble du continent. Il est composé de nombreux produits « crop_mask_<region> » cousus ensemble. Nous pouvons utiliser les produits continentaux « crop_mask » pour la plupart des applications, bien que nous devions être conscients des limites des produits régionaux qui le composent. Nous explorerons cela plus en détail ci-dessous.

[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	mask	mask	uint8	1	255.0	[crop_mask, MASK]	NaN
	prob	prob	uint8	1	255.0	[crop_prob, PROB]	NaN
	filtered	filtered	uint8	1	255.0	[mode]	NaN
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

Chargement du produit d’étendue des terres cultivées

Dans cet exemple, nous allons charger le produit 'crop_mask'.

[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',
             **query).squeeze()
print(cm)

<xarray.Dataset> Size: 2MB
Dimensions:      (y: 1011, x: 773)
Coordinates:
    time         datetime64[ns] 8B 2019-07-02T11:59:59.999999
  * y            (y) float64 8kB 1.093e+06 1.093e+06 ... 1.083e+06 1.083e+06
  * x            (x) float64 6kB 3.922e+06 3.922e+06 ... 3.93e+06 3.93e+06
    spatial_ref  int32 4B 6933
Data variables:
    mask         (y, x) uint8 782kB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
    prob         (y, x) uint8 782kB 41 41 41 45 43 45 40 ... 17 18 24 22 23 22
    filtered     (y, x) uint8 782kB 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

Représentation graphique de l’étendue des terres cultivées

Ci-dessus, nous pouvons voir que les produits « crop_mask » contiennent trois mesures :

mask: Cette bande affiche les régions cultivées sous forme de carte binaire. Les valeurs de « 1 » indiquent la présence de cultures, tandis qu’une valeur de « 0 » indique l’absence de culture. Cette bande est une carte d’étendue des terres cultivées basée sur les pixels, ce qui signifie que la carte affiche la sortie brute de la classification Random Forest basée sur les pixels.
prob : cette bande affiche les probabilités de prédiction pour la classe « culture ». Comme ce service utilise un classificateur de forêt aléatoire, les probabilités de prédiction font référence au pourcentage d’arbres qui ont voté pour la classification de forêt aléatoire. Par exemple, si le modèle avait 200 arbres de décision dans la forêt aléatoire et que 150 de ces arbres ont voté « culture », la probabilité de prédiction est de 150 / 200 x 100 = 75 %. Le fait de fixer ce seuil à > 50 % produira une carte identique à « masque ».
filtered: Cette bande affiche les régions recadrées sous forme de carte binaire. Les valeurs de 1 indiquent la présence de cultures, tandis qu’une valeur de 0 indique l’absence de recadrage. Cette bande est une carte d’étendue des terres cultivées basée sur les objets où la bande mask a été filtrée à l’aide d’un algorithme de segmentation d’image (voir ce document pour plus de détails sur l’algorithme utilisé). Au cours de ce processus, les segments plus petits que 1 Ha (100 pixels de 10 m x 10 m) sont fusionnés avec les segments voisins, ce qui donne une carte où la plus petite région classée a une taille de 1 Ha. L’ensemble de données filtré est fourni en complément de la bande mask ; les petites erreurs de commission sont supprimées par le filtrage basé sur les objets, et l’effet “sel et poivre” typique de la classification des pixels est diminué.

Ci-dessous, nous allons tracer les trois mesures côte à côte :

[7]:

fig, axes = plt.subplots(1, 3, figsize=(24, 8))
cm.mask.where(cm.mask<255).plot(ax=axes[0], # we filter to <255 to omit missing data
                   cmap='Greens',
                   add_labels=False)

cm.filtered.where(cm.filtered<255).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();

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

Exemple d’application : identifier les tendances des cultures avec le NDVI

L’indice de végétation par différence normalisée (NDVI) peut être utilisé pour suivre la santé et la croissance des cultures à mesure qu’elles mûrissent. Ici, nous chargerons une série chronologique annuelle de données satellite Sentinel-2, calculerons le NDVI et masquerons la région avec le produit d’étendue des terres cultivées de DE Africa. Avec seulement les pixels de culture restants, nous pouvons résumer le cycle de croissance des cultures dans la région.

[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

/opt/venv/lib/python3.10/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dest = _reproject(

Filtering to 122 out of 158 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 122 time steps

Calculer le NDVI

[9]:

ds = calculate_indices(ds, 'NDVI', drop=True, satellite_mission='s2')

Dropping bands ['red', 'nir']

Rééchantillonner l’ensemble de données à des intervalles de temps mensuels

[10]:

ds = ds.NDVI.resample(time='MS').mean()

#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

Masque avec la carte de l’étendue des terres cultivées

En utilisant la mesure « filtrée », nous pouvons masquer l’ensemble de données avec le produit « crop_mask ».

[11]:

#Filter out the no-data pixels (255) and non-crop pixels (0) from the crop-mask 'filtered' band and
#mask using 'filtered' band.
ds = ds.where(cm.filtered == 1)

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

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

Résumer les tendances du NDVI pour les régions de culture avec un graphique linéaire

Ci-dessous, le graphique montre que le calendrier de culture commence en avril, le pic de la saison de croissance se produisant vers la fin de l’année en octobre/novembre.

[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

Exemple d’analyse 2 : comparaison avec les ensembles de données sur la couverture terrestre mondiale

Le cube de données de DE Africa contient deux ensembles de données de couverture terrestre mondiaux d’une résolution de 10 m, chacun contenant une classe de culture. Ces ensembles de données de couverture terrestre ont tous deux été produits à l’aide d’images Senitnel-2, les mêmes que le produit d’étendue des terres cultivées de DE Afric. Ci-dessous, nous allons charger les ensembles de données de couverture terrestre sur la même région que la carte de l’étendue des terres cultivées et comparer la zone classée comme culture. Vous pouvez utiliser la carte interactive tracée ci-dessous pour déterminer les zones où chaque ensemble de données réussit bien (ou mal) à classer les cultures.

Pour en savoir plus sur les ensembles de données de couverture terrestre, consultez le bloc-notes Landcover Classification.

Tout d’abord, nous devons ressaisir certains paramètres d’analyse. Vous remarquerez ci-dessous que nous avons défini la résolution à 60 m pour nous permettre de charger rapidement une grande région, bien que tous les ensembles de données soient stockés nativement à une résolution de 10 m.

La région par défaut est une zone délimitée au-dessus du lac Kyoga, en Ouganda. Cette région est largement recadrée.

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

Tracer une carte interactive

Si vous zoomez, vous pouvez examiner les classes de couverture terrestre dominantes dans la région (indice, il y a beaucoup d’agriculture !)

[14]:

display_map(lon_range, lat_range)

[14]:

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

Charger l’étendue des terres cultivées, le produit Landcover d’ESRI et le produit WorldCover de l’ESA

[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', measurements=['mask'], **query, time='2019').squeeze()

#load the two landcover datasets for the year 2020.
lulc_esa = dc.load(product='esa_worldcover_2020', time="2020", measurements='classification', like=cm.geobox).squeeze()
lulc_esri = dc.load(product='io_lulc_v2', time="2020-07", measurements='classification', like=cm.geobox).squeeze()

Reclasser les ensembles de données de couverture terrestre en une image binaire de cultures/non-cultures pour permettre une comparaison simple avec la carte de l’étendue des terres cultivées de l’Afrique du Nord

[16]:

esri_crops = xr.where(lulc_esri['classification']==5, 1, 0)
esa_crops = xr.where(lulc_esa['classification']==40, 1, 0)

Filtrez les pixels sans données (255) de la bande « masque » du masque de recadrage.

[17]:

cm_mask = cm.mask.where(cm.mask != 255)

Tracer l’étendue du recadrage des trois ensembles de données

[18]:

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

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_41_0.png

Calculer la surface recadrée dans chaque produit

Dans cet exemple, vous verrez que le produit de terres cultivées de DE Africa cartographie beaucoup plus de cultures que les deux produits mondiaux de couverture terrestre. Dans cet exemple, le produit DE Africa est plus précis que les ensembles de données de couverture terrestre. Cependant, sur les zones humides du sud-ouest du lac Kyoga, le produit DE Afria cartographie à tort certaines zones humides comme des cultures, alors que les ensembles de données de couverture terrestre ne le font pas. Bien qu’il existe des régions où les cartes d’étendue des terres cultivées de DE AFrica comportent des erreurs, elles sont en général beaucoup plus précises que n’importe lequel des ensembles de données de couverture terrestre actuellement existants sur l’Afrique.

[19]:

pixel_length = query["resolution"][1]
m_per_km = 1000  # conversion from metres to kilometres
area_per_pixel = pixel_length**2 / m_per_km**2

[20]:

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

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_44_0.png

Limites des masques de culture régionaux

Comme nous l’avons vu dans le tableau des ensembles de données, le « crop_mask » est composé de nombreux produits régionaux assemblés. Nous pouvons charger et utiliser le produit continental « crop_mask » pour la plupart des applications. Cependant, nous pouvons observer des artefacts inhabituels aux limites des masques de culture régionaux qui méritent d’être pris en compte. Nous examinerons cela ci-dessous.

Tout d’abord, nous allons définir et inspecter une zone aux frontières de l’Ouganda, du Rwanda et de la RD Congo. Cela forme également la frontière entre « crop_mask_central » et « crop_mask_eastern ». Nous pouvons inspecter l’étendue de chaque produit régional sur la plate-forme de cartes de l’Afrique Digital Earth <https://maps.digitalearth.africa/>.

[21]:

lat, lon = -0.83, 29.58

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)

[22]:

display_map(lon_range, lat_range)

[22]:

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

[23]:

query = {
    'x': lon_range,
    'y': lat_range,
    'resolution':resolution
}

Charger les masques de culture régionaux

Ci-dessous, nous chargeons le produit continental en plus des masques de culture centraux et orientaux. Nous le faisons en utilisant l’argument « region » dans « dc.load ».

[24]:

cm = dc.load(product='crop_mask',
             measurements=['mask'],
             time='2019',
             **query).squeeze()

cm_east = dc.load(
             product='crop_mask',
             region='eastern',
             measurements=['mask'],
             time='2019',
             **query).squeeze()

cm_central = dc.load(
             product='crop_mask',
             region='central',
             measurements=['mask'],
             time='2019',
             **query).squeeze()

Produits régionaux de la parcelle

[25]:

fig, axes = plt.subplots(1, 3, figsize=(24, 8))
cm.mask.where(cm.mask<255).plot(ax=axes[0], # we filter to <255 to omit missing data
                   cmap='Greens',
                   add_labels=False)

cm_east.mask.where(cm_east.mask<255).plot(ax=axes[1],
                   cmap='Greens',
                   add_labels=False)

cm_central.mask.where(cm_central.mask<255).plot(ax=axes[2],
                   cmap='Greens',
                   add_labels=False)

axes[0].set_title('Crop mask- continental')
axes[1].set_title('Crop mask- eastern')
axes[2].set_title('Crop mask- central');

plt.tight_layout();

../../../_images/sandbox_notebooks_Datasets_Cropland_extent_52_0.png

Informations Complémentaires

Licence : Le code de ce carnet est sous licence Apache, version 2.0 <https://www.apache.org/licenses/LICENSE-2.0>. Les données de Digital Earth Africa sont sous licence Creative Commons par attribution 4.0 <https://creativecommons.org/licenses/by/4.0/>.

Contact : Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack Open Data Cube <http://slack.opendatacube.org/>`__ ou sur le GIS Stack Exchange en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment ici). Si vous souhaitez signaler un problème avec ce bloc-notes, vous pouvez en déposer un sur Github.

Version de Datacube compatible :

[26]:

print(datacube.__version__)

1.8.19

Dernier test :

[27]:

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

[27]:

'2024-11-22'