Approche OCDE/CSAO d’estimation des espaces verts urbains et du couvert arboré

• Produits utilisés : esa_worldcover_2020, esa_worldcover_2021

**Keywords**: :index:`datasets; esa_worldcover_2020`, :index:`datasets; esa_worldcover_2021`

Aperçu

Ce notebook présente l’approche du Club du Sahel et de l’Afrique de l’Ouest (CSAO, OCDE) pour estimer les indicateurs d’espaces verts urbains et de couvert arboré dans les agglomérations urbaines africaines. Les espaces verts urbains renforcent la résilience et la durabilité d’une ville au changement climatique en atténuant les phénomènes météorologiques extrêmes et les menaces graduelles comme la sécheresse et l’érosion. Ils réduisent la pollution atmosphérique, séquestrent le carbone, améliorent la qualité de l’eau et soutiennent la biodiversité et la santé physique et mentale des populations. Une analyse du CSAO met en évidence ces avantages pour les villes africaines et révèle des tendances continentales intéressantes.

L’OCDE/CSAO a créé le jeu de données Africapolis <https://docs.digitalearthafrica.org/en/latest/data_specs/Africapolis_urban_specs.html>__ afin de fournir une base de données standardisée et géospatiale indispensable sur la dynamique de l’urbanisation en Afrique, dans le but de rendre les données urbaines africaines comparables entre les pays et dans le temps. Les données d’Africapolis s’appuient sur un vaste inventaire de recensements de l’habitat et de la population, de registres électoraux et d’autres sources officielles de données démographiques, remontant parfois au début du XXème siècle.

Ces ensembles de données (africapolis_2015 et africapolis_2020) sont partagés avec des indicateurs d’espaces verts urbains, c’est-à-dire le pourcentage de la zone d’agglomération urbaine qui a une couverture verte (arbres, arbustes et herbe), ainsi que le pourcentage de couverture arborée. Ces indicateurs ont été estimés à l’aide des informations sur la couverture terrestre des bases de données ESA World Cover v100 (2020) et ESA World Cover v200 (2021). L’ensemble de données ESA World Cover fournit des informations sur les types de couverture terrestre, permettant le calcul des fractions de couverture terrestre telles que les espaces verts urbains, les zones bâties, les terres ouvertes et les eaux libres. Les statistiques zonales permettent d’établir des histogrammes de fréquence des classes d’occupation du sol au sein de chaque agglomération, puis de les traiter pour calculer ces indicateurs. L’espace vert urbain est calculé comme la somme de la superficie des classes couvert arboré, arbustaie et prairie divisée par la superficie totale de l’agglomération. Le pourcentage de couvert arboré est calculé en divisant la superficie de la classe d’occupation arborée par la superficie totale de l’agglomération. Pour plus d’informations et d’analyses issues de cet indicateur, consultez la plateforme « Mapping Territorial Transformations in Africa » <https://www.mapping-africa-transformations.org/climate/>.

Ce notebook présente la méthodologie utilisée par le SWAC pour calculer les indicateurs d’espaces verts urbains et de couvert arboré pour l’ensemble de données Africapolis. Il permet également aux utilisateurs de reproduire cette analyse pour leurs propres agglomérations urbaines ne faisant pas partie de l’ensemble de données Africapolis. Pour vérifier si une agglomération spécifique est déjà incluse dans l’ensemble de données Africapolis, consultez le carnet « Africapolis_urban_agglomerations <../Datasets/Africapolis_urban_agglomerations.ipynb> ».

Référence

Anderson, B., J.E. Patiño Quinchia, R. Prieto-Curiel (2022), « Renforcer la résilience des villes africaines au changement climatique : le rôle des espaces verts », Notes ouest-africaines, n° 37, Éditions OCDE, Paris. https://doi.org/1787/3303cfb3-en

Description

Ce cahier est conçu pour calculer les indicateurs relatifs aux espaces verts urbains en Afrique. Il suit les étapes suivantes.

1. Sélectionner/charger un emplacement pour l’analyse.

2. Calculer la fréquence d’occurrence de la couverture terrestre dans chaque agglomération urbaine d’intérêt

3. Calculez les fractions de couverture terrestre et les regroupezs en catégories (c.-à-d. espaces verts, couverture arborée).r

Charger des paquets

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

import datacube
import rasterio
import numpy as np
import pandas as pd
import geopandas as gpd

from tqdm.auto import tqdm
import matplotlib.pyplot as plt
from datacube.utils import geometry
from odc.dscache.tools import tiling
from shapely.geometry import Polygon, shape
from odc.geo.geom import Geometry

from deafrica_tools.spatial import xr_rasterize
from deafrica_tools.areaofinterest import define_area
from deafrica_tools.plotting import display_map, plot_lulc

Se connecter au datacube

Connectez-vous au datacube pour que nous puissions accéder aux données de DE Africa. Le paramètre « app » est un nom unique pour l’analyse qui est basé sur le nom du fichier du notebook.

dc = datacube.Datacube(app="OECD_SWAC_Urban_green_space")

Paramètres d’analyse

Sélectionnez l’emplacement

La cellule suivante définit les paramètres qui définissent la zone d’intérêt sur laquelle effectuer l’analyse. #### Sélectionner l’emplacement Pour définir la zone d’intérêt, deux méthodes sont disponibles :

1. En spécifiant la latitude, la longitude et la zone tampon, ou en séparant les zones tampons de latitude et de longitude, cette méthode vous permet de définir une zone d’intérêt autour d’un point central. Vous pouvez saisir la latitude centrale, la longitude centrale et une valeur de zone tampon en degrés pour créer une zone carrée autour du point central. Par exemple, « lat = 10.338 », « lon = -1.055 » et « buffer = 0.1 » sélectionneront une zone avec un rayon de 0,1 degré carré autour du point avec les coordonnées « (10.338, -1.055) ».

   Vous pouvez également fournir des valeurs de tampon distinctes pour la latitude et la longitude d’une zone rectangulaire. Par exemple, « lat = 10.338 », « lon = -1.055 » et « lat_buffer = 0.1 » et « lon_buffer = 0.08 » sélectionneront une zone rectangulaire s’étendant sur 0,1 degré au nord et au sud et 0,08 degré à l’est et à l’ouest à partir du point « (10.338, -1.055) ».

   Pour des temps de chargement raisonnables, définissez la mémoire tampon sur « 0,1 » ou moins.

2. By uploading a polygon as a GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.

Pour utiliser l’une de ces méthodes, vous pouvez décommenter la ligne de code concernée et commenter l’autre. Pour commenter une ligne, ajoutez le symbole "#" avant le code que vous souhaitez commenter. Par défaut, la première option qui définit l’emplacement à l’aide de la latitude, de la longitude et du tampon est utilisée.

Si vous exécutez le notebook pour la première fois, conservez les paramètres par défaut ci-dessous. Cela illustrera le fonctionnement de l’analyse et fournira des résultats significatifs. L’exemple porte sur une zone autour de Kumasi, dans le sud du Ghana, le même emplacement que celui sélectionné dans le notebook du jeu de données « Africapolis_urban_agglomerations <../Datasets/Africapolis_urban_agglomerations.ipynb> », car ce bloc-notes illustre l’approche utilisée pour calculer les indicateurs d’espaces verts et de couvert arboré dans le jeu de données Africapolis.

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=6.6975, lon=-1.6286, buffer=0.2)

# Method 2: Use a polygon as a GeoJSON or Esri Shapefile. This can include multiple polygons.
# aoi = define_area(vector_path='aoi.shp')

#Create a geopolygon and geodataframe of the area of interest
geopolygon_list = [Geometry(feature['geometry'], crs="EPSG:4326") for feature in aoi['features']]
agglomeration = gpd.GeoDataFrame(geometry=geopolygon_list, crs="EPSG:4326")
agglomeration['agglosID'] = range(1, len(agglomeration) + 1)

# Get the latitude and longitude range of the geopolygon
lat_range = (agglomeration.total_bounds[1], agglomeration.total_bounds[3])
lon_range = (agglomeration.total_bounds[0], agglomeration.total_bounds[2])
# display_map(x=lon_range, y=lat_range)

agglomeration.explore(tiles="https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
                attr='http://mt0.google.com/vt/lyrs=y&hl=en&x={x}&y={y}&z={z}')

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

Effectuer des statistiques zonales à l’aide d’ESA Worldcover (2020 et 2021)

Cette section du notebook effectue des statistiques zonales sur les polygones d’agglomération urbaine à l’aide des données ESA WorldCover de 2020 et 2021. Chaque polygone est traité individuellement à l’aide d’une boucle : les données pertinentes pour chaque année sont extraites, masquées à la limite du polygone, puis les histogrammes de fréquence des classes d’occupation du sol au sein de cette zone sont calculés. Ces histogrammes fournissent des informations sur la répartition des différents types d’occupation du sol au sein de l’agglomération. Les résultats sont stockés pour analyse et comparaison ultérieures.

Définir les paramètres

Définissez les paramètres clés nécessaires au traitement des données ESA WorldCover. La variable « classes » représente les différents types d’occupation du sol selon les classifications ESA WorldCover.

# Define resolution, measurements, and classes
resolution = (-10, 10)
measurements = ['classification']
classes = [10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]

# Define the years of interest
years = ['2020', '2021']

# Zonal stats results container
zs_results = []

Définir la fonction de traitement des données de l’ESA

Cette fonction extrait les données d’occupation du sol ESA WorldCover dans les polygones d’agglomération spécifiés. Elle effectue une opération statistique zonale qui calcule la distribution des classes d’occupation du sol au sein de chaque polygone.

def process_esa_data(query, row, idx, year):
    try:
        # Load the ESA data for the current polygon and year
        ds_esa = dc.load(product=f'esa_worldcover_{query["time"]}', **query).squeeze()

        # Rasterize the area of interest polygon
        aoi_raster = xr_rasterize(
            gdf=gpd.GeoDataFrame(geometry=[row['geometry']], crs="EPSG:4326"),
            da=ds_esa,
            crs=ds_esa.crs
        )

        # Mask the dataset to the rasterized area of interest
        ds_esa = ds_esa.where(aoi_raster == 1)

        # Calculate frequency histogram of land cover classes
        lc_array = ds_esa.classification.values.flatten()
        lc_classes, lc_counts = np.unique(lc_array, return_counts=True)
        lc_frequency = dict(zip(lc_classes, lc_counts))

        # Remove NaN values if present
        lc_frequency.pop(np.nan, None)

        return lc_frequency

    except Exception as e:
        print(f"Error processing polygon {idx} for year {year}: {e}")
        return {}

Itérer sur les polygones d’agglomérations

Pour chaque polygone d’agglomération, extrayez les données ESA WorldCover pour 2020 et 2021. Les résultats sont stockés dans une liste contenant des dictionnaires, où chaque entrée correspond à un polygone d’agglomération.

# Iterate over each row in the GeoDataFrame with a progress bar
for idx, row in tqdm(agglomeration.iterrows(), total=agglomeration.shape[0]):
    # Define query parameters based on the polygon's geometry
    geom = geometry.Geometry(geom=row['geometry'], crs="EPSG:4326")

    # Initialize a dictionary to hold histograms for each year
    yearly_histograms = {}

    # Process each year
    for year in years:
        query = {
            'time': year,  # Select the year of interest
            'geopolygon': geom,
            'resolution': resolution,
            'output_crs': 'EPSG:6933',
            'measurements': measurements
        }

        # Load and process ESA WorldCover data
        yearly_histograms[year] = process_esa_data(query, row, idx, year)

    # Append results to zs_results
    zs_results.append({'agglosID': row.agglosID, 'geometry': row['geometry'], 'histograms': yearly_histograms})

Convertir les résultats en DataFrame

Les résultats collectés à l’étape précédente sont convertis en un DataFrame Pandas. Cela facilite le traitement et l’analyse des statistiques d’occupation du sol par agglomération urbaine.

# Convert results to DataFrame
zs_results_df = pd.DataFrame(zs_results)

# Display first few rows
print(zs_results_df.head())

   agglosID                                           geometry  \
0         1  POLYGON ((-1.4286 6.4975, -1.4286 6.8975, -1.8...

                                          histograms
0  {'2020': {10.0: 9281028, 20.0: 2103586, 30.0: ...

Traiter les données de l’histogramme dans des colonnes séparées et calculer les fractions

Les résultats des statistiques zonales de WorldCover de l’ESA permettent de calculer les fractions des classes d’occupation du sol au sein de chaque polygone. Ces fractions sont calculées en divisant chaque classe par le nombre total, puis les résultats sont triés et sélectionnés pour chaque polygone.

Traiter les histogrammes dans des colonnes distinctes pour chaque année

Chaque ligne de « zs_results_df » contient un dictionnaire d’histogrammes de classe de couverture terrestre pour 2020 et 2021. Cette étape extrait les valeurs d’histogramme pour chaque année séparément et les convertit en colonnes DataFrame.

def process_histogram_per_year(row, year):
    """Extract histogram values for a given year and return as a Series."""
    histogram = row['histograms'].get(year, {})
    return pd.Series(histogram)

# Process histograms for both years (2020 and 2021)
histograms_2020_df = zs_results_df.apply(lambda row: process_histogram_per_year(row, '2020'), axis=1).fillna(0).astype(int)
histograms_2021_df = zs_results_df.apply(lambda row: process_histogram_per_year(row, '2021'), axis=1).fillna(0).astype(int)

Réindexer les DataFrames pour inclure toutes les classes de couverture terrestre

Les histogrammes peuvent ne pas inclure toutes les classes d’occupation du sol possibles, notamment si certaines classes sont absentes dans certaines zones. Nous nous assurons ici que toutes les classes sont présentes en réindexant les DataFrames et en remplissant les valeurs manquantes avec des zéros. Les colonnes sont renommées pour plus de clarté.

# Reindex to include all classes (fill missing with 0)
histograms_2020_df = histograms_2020_df.reindex(columns=classes, fill_value=0)
histograms_2021_df = histograms_2021_df.reindex(columns=classes, fill_value=0)

# Rename columns for clarity
histograms_2020_df.rename(columns={k: f'lc_{k}_count_2020' for k in classes}, inplace=True)
histograms_2021_df.rename(columns={k: f'lc_{k}_count_2021' for k in classes}, inplace=True)

Combiner les données de l’histogramme avec le DataFrame d’origine

Maintenant que les histogrammes sont correctement formatés, nous les fusionnons avec le DataFrame d’origine. Seuls « agglosID » et « geometry » du DataFrame d’origine sont conservés avec les données d’histogramme traitées.

zs_ESA_LC = pd.concat([zs_results_df[['agglosID', 'geometry']], histograms_2020_df, histograms_2021_df], axis=1)

Calculer les nombres totaux et les fractions

Pour chaque année, le nombre total de pixels classés est calculé. La fraction de chaque classe d’occupation du sol est ensuite calculée en pourcentage du nombre total de pixels. Une multiplication par 100 convertit ces fractions en pourcentages.

# Calculate total counts and fractions for both years
for year in ['2020', '2021']:
    zs_ESA_LC[f'total_count_{year}'] = zs_ESA_LC.filter(like=f'_count_{year}').sum(axis=1)
    for k in classes:
        zs_ESA_LC[f'lc_{k}_fraction_{year}'] = zs_ESA_LC[f'lc_{k}_count_{year}'] / zs_ESA_LC[f'total_count_{year}']

# Multiply fractions by 100 to convert them into percentages
for year in ['2020', '2021']:
    for k in classes:
        zs_ESA_LC[f'lc_{k}_fraction_{year}'] *= 100

Organiser et formater le DataFrame final

Pour plus de clarté, nous avons explicitement ordonné les colonnes de fractions pour les deux années. Le DataFrame final inclut « agglosID », « geometry » et toutes les valeurs de fractions, arrondies à deux décimales.

# Explicitly order the fraction columns for each year
fraction_columns_ordered_2020 = [f"lc_{k}_fraction_2020" for k in classes]
fraction_columns_ordered_2021 = [f"lc_{k}_fraction_2021" for k in classes]

# Combine ordered columns into the final DataFrame
zs_ESA_LC_selected = zs_ESA_LC[['agglosID', 'geometry'] + fraction_columns_ordered_2020 + fraction_columns_ordered_2021].fillna(0).round(2)

# Display the final DataFrame
zs_ESA_LC_selected.head()

	agglosID	geometry	lc_10_fraction_2020	lc_20_fraction_2020	lc_30_fraction_2020	lc_40_fraction_2020	lc_50_fraction_2020	lc_60_fraction_2020	lc_70_fraction_2020	lc_80_fraction_2020	...	lc_20_fraction_2021	lc_30_fraction_2021	lc_40_fraction_2021	lc_50_fraction_2021	lc_60_fraction_2021	lc_70_fraction_2021	lc_80_fraction_2021	lc_90_fraction_2021	lc_95_fraction_2021	lc_100_fraction_2021
0	1	POLYGON ((-1.4286 6.4975, -1.4286 6.8975, -1.8...	47.45	10.75	12.6	0.28	19.72	8.86	0.0	0.32	...	6.39	15.34	0.01	26.09	0.41	0.0	0.35	0.02	0.0	0.0

1 rows × 24 columns

Comparaison de la fraction de couverture terrestre pour 2020 et 2021

Nous comparons les fractions de couverture terrestre pour différentes classes entre 2020 et 2021. Chaque agglomération est analysée et les fractions de ces types de couverture terrestre sont calculées pour les deux années.

Le graphique à barres ci-dessous compare ces fractions pour chaque agglomération, les barres de 2020 étant bleu ciel et celles de 2021, saumon. Ce graphique illustre visuellement l’évolution de la couverture terrestre entre les deux années, mettant en évidence les changements notables dans l’utilisation des terres ou les tendances environnementales.

# List of land cover class names for the x-axis
land_cover_labels = [
    "Trees",     # lc_10_fraction
    "Shrubland", # lc_20_fraction
    "Grassland", # lc_30_fraction
    "Cropland",  # lc_40_fraction
    "Built-up",  # lc_50_fraction
    "Barren / Sparse Vegetation", # lc_60_fraction
    "Snow and Ice", # lc_70_fraction
    "Open Water", # lc_80_fraction
    "Herbaceous Wetland", # lc_90_fraction
    "Mangroves",  # lc_95_fraction
    "Moss and Lichen", # lc_100_fraction
]

# Bar width and x-axis positions
bar_width = 0.4  # width of each bar
index = np.arange(len(land_cover_labels))  # positions for the bars on the x-axis

# Iterate through the rows in the final DataFrame and plot histograms
for idx, row in zs_ESA_LC_selected.iterrows():
    agglosID = row['agglosID']

    # Extract the counts for 2020 and 2021
    lc_counts_2020 = [row[f'lc_{k}_fraction_2020'] for k in classes]
    lc_counts_2021 = [row[f'lc_{k}_fraction_2021'] for k in classes]

    # Create a figure for the combined plot
    plt.figure(figsize=(12, 6))

    # Plot for 2020 (shifted to the left)
    plt.bar(index - bar_width / 2, lc_counts_2020, bar_width, alpha=0.7, color='lightgreen', label='2020')

    # Plot for 2021 (shifted to the right)
    plt.bar(index + bar_width / 2, lc_counts_2021, bar_width, alpha=0.7, color='lightcoral', label='2021')

    # Add titles and labels
    plt.title(f'Land cover fractions comparison for 2020 & 2021 - Agglomeration {agglosID}')
    plt.xlabel('Land cover fraction')
    plt.ylabel('Percentage (%)')
    plt.xticks(index, land_cover_labels, rotation=45, ha='right')
    plt.legend()

    # Display the plot
    plt.tight_layout()
    plt.show()

Calculer les indicateurs d’espaces verts et de couvert arboré

Calculer des indicateurs indirects pour les espaces verts et le couvert arboré en fonction des fractions de classes spécifiques de couverture terrestre. L’indicateur « Espaces verts urbains » est obtenu en additionnant les fractions des types de couverture terrestre associés à la végétation (arbres, arbustes et prairies) pour 2020 et 2021. L’indicateur « Couvert arboré » se concentre spécifiquement sur la fraction de terre couverte par les arbres. Ces indicateurs sont calculés pour chaque agglomération et exprimés en pourcentages afin de faciliter la comparaison entre les différentes zones urbaines.

# Calculate specific fractions for 2020 and convert to percentages
zs_ESA_LC['p_Tree_cover_WC2020'] = (zs_ESA_LC['lc_10_fraction_2020']).round(1)
zs_ESA_LC['p_Urban_green_space_WC2020'] = (zs_ESA_LC[['lc_10_fraction_2020', 'lc_20_fraction_2020', 'lc_30_fraction_2020']].sum(axis=1)).round(1)

# Calculate specific fractions for 2021 and convert to percentages
zs_ESA_LC['p_Tree_cover_WC2021'] = (zs_ESA_LC['lc_10_fraction_2021']).round(1)
zs_ESA_LC['p_Urban_green_space_WC2021'] = (zs_ESA_LC[['lc_10_fraction_2021', 'lc_20_fraction_2021', 'lc_30_fraction_2021']].sum(axis=1)).round(1)

# Select and rename columns
green_space_indicators = zs_ESA_LC[['agglosID',
                             'p_Urban_green_space_WC2020', 'p_Tree_cover_WC2020',
                             'p_Urban_green_space_WC2021', 'p_Tree_cover_WC2021', 'geometry']].copy()

# Convert the resulting DataFrame into a GeoDataFrame to handle spatial data, setting the geometry column appropriately and defining the CRS as "EPSG:4326"
green_space_indicators = gpd.GeoDataFrame(green_space_indicators, geometry='geometry', crs="EPSG:4326")
# Display final DataFrame
green_space_indicators.head()

	agglosID	p_Urban_green_space_WC2020	p_Tree_cover_WC2020	p_Urban_green_space_WC2021	p_Tree_cover_WC2021	geometry
0	1	70.8	47.4	73.1	51.4	POLYGON ((-1.4286 6.4975, -1.4286 6.8975, -1.8...

Visualiser les indicateurs d’espaces verts urbains et de couverture arborée

Visualisez la comparaison des espaces verts urbains et du couvert arboré entre 2020 et 2021 pour chaque agglomération sélectionnée. Les visualisations comprennent une carte affichant les limites des agglomérations, des cartes superposées des classifications ESA WorldCover pour les deux années afin de montrer l’évolution de la couverture terrestre, et un graphique à barres groupées comparant le pourcentage d’espaces verts urbains et de couvert arboré au cours des deux années.

# Create a figure with dynamic subplots (3 columns: map, stacked WC, and bar)
num_agglomerations = len(green_space_indicators)
fig = plt.figure(figsize=(20, 10 * num_agglomerations))  # Set figure size

# Iterate over selected agglomerations and their corresponding subplot axes
for i, (index, agglomeration) in enumerate(green_space_indicators.iterrows()):
    # Subplot layout (3 columns: map, stacked WC, and bar)
    # Map subplot (spans 2 rows)
    ax_map = plt.subplot2grid((num_agglomerations * 2, 3), (i * 2, 0), rowspan=2)
    # ESA WorldCover 2020 (top half of the second column)
    ax_wc_2020 = plt.subplot2grid((num_agglomerations * 2, 3), (i * 2, 1))
    # ESA WorldCover 2021 (bottom half of the second column)
    ax_wc_2021 = plt.subplot2grid((num_agglomerations * 2, 3), (i * 2 + 1, 1))
    # Bar chart subplot (spans 2 rows)
    ax_bar = plt.subplot2grid((num_agglomerations * 2, 3), (i * 2, 2), rowspan=2)
    # Extract the agglomeration as a GeoDataFrame
    agglomeration_gdf = green_space_indicators.iloc[i:i+1]

    # Plot the agglomeration polygon on the map axis
    ax_map.set_title(f"Agglomeration {agglomeration['agglosID']}")
    agglomeration_gdf.plot(ax=ax_map, edgecolor='black')
    ax_map.set_xlabel("Longitude")
    ax_map.set_ylabel("Latitude")

    # Get the top-left corner of the bounding box
    bbox = agglomeration_gdf.total_bounds
    top_left_x, top_left_y = bbox[0], bbox[3]

    # Prepare data for the bar chart
    labels = ['Urban green space', 'Tree cover']
    wc2020 = [
        agglomeration['p_Urban_green_space_WC2020'],
        agglomeration['p_Tree_cover_WC2020']
    ]
    wc2021 = [
        agglomeration['p_Urban_green_space_WC2021'],
        agglomeration['p_Tree_cover_WC2021']
    ]

    # Plot grouped bar chart
    bar_width = 0.35
    index = np.arange(len(labels))  # Position of the bars on the x-axis
    ax_bar.bar(index - bar_width/2, wc2020, bar_width, label='2020', color='lightgreen')
    ax_bar.bar(index + bar_width/2, wc2021, bar_width, label='2021', color='lightcoral')

    # Add labels and title to the bar chart
    ax_bar.set_xlabel('Green space types')
    ax_bar.set_ylabel('Percentage (%)')
    ax_bar.set_title(f"Green space comparison")
    ax_bar.set_xticks(index)
    ax_bar.set_xticklabels(labels)
    ax_bar.legend()

    # Load ESA WorldCover Data for each agglomeration
    query = {
        'y': (agglomeration_gdf.total_bounds[1], agglomeration_gdf.total_bounds[3]),
        'x': (agglomeration_gdf.total_bounds[0], agglomeration_gdf.total_bounds[2]),
        'resolution': (10, -10),
        'output_crs': 'epsg:6933',
        'measurements': 'classification'
    }

    # Load ESA WorldCover data for both 2020 and 2021
    esa_wc = dc.load(product=['esa_worldcover_2020', 'esa_worldcover_2021'], **query).squeeze()

    # Rasterize the area of interest polygon for masking ESA WorldCover data
    aoi_raster = xr_rasterize(gdf=agglomeration_gdf, da=esa_wc, crs=esa_wc.crs)

    # Mask ESA WorldCover data based on the agglomeration boundary
    esa_wc_2020_masked = esa_wc.sel(time="2020").squeeze()['classification'].where(aoi_raster == 1)
    esa_wc_2021_masked = esa_wc.sel(time="2021").squeeze()['classification'].where(aoi_raster == 1)

    # Plot ESA WorldCover 2020 (top part of the stack)
    plot_lulc(esa_wc_2020_masked, product="ESA", legend=True, ax=ax_wc_2020)
    ax_wc_2020.set_title("ESA WorldCover 2020")

    # Plot ESA WorldCover 2021 (bottom part of the stack)
    plot_lulc(esa_wc_2021_masked, product="ESA", legend=True, ax=ax_wc_2021)
    ax_wc_2021.set_title("ESA WorldCover 2021")

# Add a title for the entire figure
fig.suptitle(f"Urban green space and tree cover comparison", fontsize=18, fontweight='bold')

# Adjust layout and show the plot
plt.tight_layout()
plt.subplots_adjust(top=0.90)
plt.show()

Interprétation du graphique final

La figure générée fournit une analyse comparative des espaces verts urbains et du couvert arboré au sein des agglomérations sélectionnées. Voici comment interpréter chaque composante :

1. Vue cartographique (première colonne)

   • Affiche la limite de l’agglomération pour fournir un contexte spatial.

   • Aide à identifier l’étendue des zones urbaines analysées pour les espaces verts et la couverture arborée.

2. Couverture terrestre ESA WorldCover (deuxième colonne)

   • Deux graphiques empilés montrent la classification de la couverture terrestre pour 2020 et 2021.

   • Tout changement notable dans l’espace vert ou la couverture arborée au fil du temps peut être inspecté visuellement.

3. Comparaison du graphique à barres (troisième colonne)

   • Affiche le pourcentage d’espaces verts urbains et de couverture arborée en 2020 et 2021.

   • Observations clés :

     • Ces augmentations indiquent des améliorations de la végétation et des efforts potentiels de reboisement.

     • Les diminutions peuvent suggérer une expansion urbaine, une déforestation ou des changements d’utilisation des terres.

     • Des valeurs cohérentes peuvent indiquer une couverture végétale stable.

Conclusions et recommandations

Ce carnet présente la méthodologie utilisée par le Club du Sahel et de l’Afrique de l’Ouest (CSAO, OCDE) pour calculer les indicateurs d’espaces verts urbains et de couvert arboré des agglomérations africaines. Pour les agglomérations non incluses dans la base de données Africapolis, ce carnet offre un outil précieux pour reproduire l’analyse et obtenir des données comparables. L’examen de ces indicateurs permet de comparer l’état de la verdure urbaine entre différentes régions et périodes, contribuant ainsi à éclairer les politiques, l’urbanisme et les stratégies environnementales, conformément aux objectifs de durabilité et de résilience climatique.

Pour les agglomérations non incluses dans l’ensemble de données Africapolis, consultez le notebook« Africapolis_urban_agglomerations <../Datasets/Africapolis_urban_agglomerations.ipynb> » pour vérifier si votre agglomération est disponible. Pour des données plus récentes et plus précises sur la couverture arborée, pensez à utiliser le notebook « Urban_tree_extent <Urban_tree_extent.ipynb> », qui utilise la géomatique Sentinel-2 de DE Africa pour analyser les variations annuelles de la couverture arborée, fournissant ainsi des informations plus récentes que les ensembles de données WorldCover de l’ESA de 2020 et 2021.

---

Licence Le code de ce notebook 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 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 DE Africa <https://digitalearthafrica.slack.com/> ou sur le GIS Stack Exchange <https://gis.stackexchange.com/questions/ask?tags=open-data-cube> en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment `ici <https://gis.stackexchange.com/questions/tagged/open-data-

Si vous souhaitez signaler un problème avec ce notebook, vous pouvez en déposer un sur Github.

Version de Datacube compatible

print(datacube.__version__)

1.8.20

Dernier test :

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

'2025-01-31'