Cartographie de l’étendue des arbres urbains
Produits utilisés : gm_s2_annual
Aperçu
Les espaces verts urbains, et en particulier les arbres, jouent un rôle essentiel dans l’amélioration de la qualité de vie en ville en fournissant des services écosystémiques tels que la purification de l’air, la régulation de la température et des espaces de loisirs pour les habitants. Une cartographie et un suivi précis de l’étendue des arbres urbains sont essentiels pour les urbanistes et les gestionnaires de l’environnement afin de garantir un développement urbain durable. La cartographie du couvert arboré urbain permet d’analyser leur potentiel en tant que solutions naturelles pour atténuer les chaleurs extrêmes, comme illustré ici <https://www.mapping-africa-transformations.org/green-spaces-nbs-heat-waves/>`__.
Ce travail s’inscrit dans le cadre d’une initiative collaborative entre Digital Earth Africa et le Club du Sahel et de l’Afrique de l’Ouest (CSAO) de l’OCDE visant à développer des outils et des méthodologies pour la durabilité urbaine en Afrique. En exploitant l’imagerie satellite de Sentinel-2, ce partenariat vise à fournir des informations fiables et fondées sur des données pour la cartographie et le suivi de la végétation urbaine dans les villes africaines. L’imagerie satellite, notamment celle issue de capteurs haute résolution comme Sentinel-2, constitue un outil précieux pour la cartographie et le suivi de la végétation en zone urbaine. Sentinel-2 fournit des données multispectrales d’une résolution adaptée à la détection et à l’analyse des arbres urbains à l’échelle de la ville. Ces données, combinées à des techniques de classification avancées, peuvent contribuer à la délimitation précise des espaces verts urbains.
Ce notebook complète les travaux de l’OCDE/CSAO, qui fournit des données sur les agglomérations urbaines (pour 2015 et 2020) ainsi que des indicateurs d’espaces verts, basés sur les jeux de données WorldCover de l’ESA pour 2020 et 2021, via le jeu de données Africapolis <https://docs.digitalearthafrica.org/en/latest/data_specs/Africapolis_urban_specs.html>. Cette analyse s’appuie sur ces données en proposant une méthode plus détaillée de cartographie et de suivi de l’étendue des arbres urbains à l’aide des données Sentinel-2, offrant ainsi un aperçu plus approfondi de la végétation urbaine au fil du temps. En intégrant l’analyse des séries chronologiques, nous évaluons l’évolution du couvert arboré dans les agglomérations urbaines, contribuant ainsi à une compréhension globale de la durabilité urbaine.
Dans ce notebook, nous mettons en œuvre un algorithme pour cartographier l’étendue des arbres urbains à l’aide du « Digital Earth Africa Sentinel-2 Annual GeoMAD <https://docs.digitalearthafrica.org/en/latest/data_specs/GeoMAD_specs.html> » , intégrant une analyse des séries chronologiques pour évaluer les changements dans la couverture arborée sur plusieurs années.
Description
Ce notebook montre comment charger et visualiser l’étendue des arbres urbains. Il couvre :
Définir la zone d’intérêt (AOI) : La méthode par défaut calcule l’étendue de l’arbre urbain pour la zone d’intérêt (AOI) que vous définissez. Cette AOI permet de sélectionner les agglomérations urbaines dans le jeu de données Africapolis <https://docs.digitalearthafrica.org/en/latest/data_specs/Africapolis_urban_specs.html>`__. Cependant, si vous ne souhaitez pas utiliser d’agglomération et préférez utiliser uniquement la AOI que vous avez définie, vous pouvez commenter les cellules 4 à 7. Cela vous permettra de travailler exclusivement avec la zone spécifiée sans sélectionner d’agglomération.
Charger les données Sentinel-2 GeoMAD : récupérez et chargez les données annuelles Sentinel-2 GeoMAD couvrant la zone d’intérêt.
Calculer le NDVI : Calculez l’indice de végétation par différence normalisée (NDVI) pour mettre en évidence les zones avec une végétation saine susceptibles de contenir des arbres.
Appliquer le seuillage Otsu pour la détection des arbres : utilisez le seuillage Otsu à plusieurs niveaux sur les valeurs NDVI pour segmenter les arbres urbains des autres types de couverture végétale.
Calculer la superficie des arbres : Calculez la superficie totale couverte d’arbres en kilomètres carrés en comptant le nombre de pixels classés comme arbres au fil du temps.
Graphique de série chronologique : Tracez une série chronologique montrant la superficie totale des arbres couverte au fil du temps pour analyser les tendances et les changements dans la couverture arborée au cours de la période d’étude.
Charger des paquets
Importez les packages Python utilisés pour l’analyse.
[1]:
import datacube
import rioxarray
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.patches as mpatches
from odc.geo.geom import Geometry
from skimage.filters import threshold_multiotsu
from deafrica_tools.plotting import rgb, display_map
from deafrica_tools.areaofinterest import define_area
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.spatial import xr_rasterize, xr_vectorize
from deafrica_tools.load_africapolis import get_africapolis
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.
[2]:
dc = datacube.Datacube(app="OECD_SWAC_Urban_tree_extent")
Télécharger ou sélectionner une ou plusieurs agglomérations urbaines
La cellule suivante vous permet de définir la zone d’intérêt pour l’analyse de l’étendue des arbres urbains. #### Sélectionner l’emplacement Pour définir la zone d’intérêt, deux méthodes sont disponibles :
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.
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 Kinshasa, capitale de la République démocratique du Congo.
Cette zone d’intérêt (AOI) peut être utilisée directement pour l’analyse de l’étendue des arbres urbains. Cependant, dans les cellules suivantes, vous verrez comment cette même AOI peut servir de point de départ pour sélectionner des agglomérations urbaines dans le jeu de données Africapolis si vous préférez travailler avec une agglomération spécifique. Notez que les lignes de code pour la sélection depuis Africapolis (lignes 4 à 7) sont incluses par défaut. Si vous ne souhaitez pas utiliser d’agglomération et préférez travailler uniquement avec l’AOI définie, commentez simplement les lignes 4 à 7 en ajoutant le symbole « # » au début de chacune d’elles.
[3]:
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-4.3615, lon=15.3088, buffer=0.05)
# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')
#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)
# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])
display_map(lon_range, lat_range)
[3]:
Alternative : Sélectionner les polygones d’agglomération urbaine AFRICAPOLIS dans la zone d’intérêt
Si vous préférez travailler avec une agglomération urbaine du jeu de données AFRICAPOLIS plutôt qu’avec l’AOI précédemment définie, le code ci-dessous propose cette alternative. Il chargera les agglomérations de l’AOI définie précédemment.
Si vous préférez travailler avec une agglomération urbaine de l’ensemble de données AFRICAPOLIS au lieu d’utiliser la zone d’intérêt (AOI) par défaut, vous pouvez utiliser le code ci-dessous pour charger les agglomérations dans l’AOI.
Par défaut, le code sélectionne les agglomérations du jeu de données Africapolis au sein de l’AOI. Cependant, si vous ne souhaitez pas utiliser l’agglomération, vous pouvez choisir d’utiliser l’AOI précédemment définie. Pour utiliser cette option, commentez simplement les lignes 4 à 7, qui sont responsables du chargement des données Africapolis et de la sélection de l’agglomération, en ajoutant un symbole « # » au début de ces lignes.
Cela vous permettra soit de procéder à une analyse basée sur l’agglomération, soit de vous en tenir à l’AOI d’origine que vous avez définie.
Dans les étapes suivantes, nous sélectionnerons l’agglomération de Kinshasa à partir de l’ensemble de données Africapolis <https://docs.digitalearthafrica.org/en/latest/data_specs/Africapolis_urban_specs.html>`__ pour notre analyse.
[4]:
# Create a bounding box from study area coordinates
bbox = (lon_range[0], lat_range[0], lon_range[1], lat_range[1])
# Load Africapolis data.
agglomerations = get_africapolis(bbox=bbox, layer='Africapolis_2020')
# Print list of urban agglomerations and their green urban green space indicators
agglomerations.head()
[4]:
| id | agglosID | agglosName | ISO3 | Pop2020 | p_Urban_green_space_WC2020 | p_Tree_cover_WC2020 | p_Urban_green_space_WC2021 | p_Tree_cover_WC2021 | Longitude | Latitude | color_HEX | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Africapolis_2020.2158 | 2158 | Brazzaville | COG | 2031742 | 43.5 | 13.5 | 45.8 | 19.5 | 15.245336 | -4.220031 | #f1b78c | MULTIPOLYGON (((15.2308 -4.3103, 15.231 -4.310... |
| 1 | Africapolis_2020.4858 | 4858 | Kinshasa | COD | 9850000 | 41.0 | 17.6 | 40.6 | 20.2 | 15.296038 | -4.410428 | #f8dca0 | MULTIPOLYGON (((15.2286 -4.3333, 15.2286 -4.33... |
[5]:
# Explore the urban agglomeration polygons located within the bounding box in an interactive map
agglomerations.explore(tiles="https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
color=agglomerations['color_HEX'],
attr='http://mt0.google.com/vt/lyrs=y&hl=en&x={x}&y={y}&z={z}')
[5]:
[6]:
# Edit the name below to select the area of interest
region_name = "Kinshasa"
geopolygon_gdf = agglomerations[agglomerations["agglosName"] == region_name]
# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])
# Explore the urban agglomeration polygons located within the bounding box
geopolygon_gdf.explore(tiles="https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
color=agglomerations['color_HEX'],
attr='http://mt0.google.com/vt/lyrs=y&hl=en&x={x}&y={y}&z={z}')
[6]:
Calculer l’étendue des arbres urbains dans l’agglomération urbaine
Les paramètres de sélection des données annuelles GeoMAD de Sentinel-2 sont définis pour extraire des informations spécifiques à l’agglomération urbaine sélectionnée. Seules les quatre bandes (rouge, verte, bleue et proche infrarouge (nir)) sont sélectionnées, car elles sont essentielles au traçage RVB et au calcul du NDVI. Les bandes RVB visuelles (rouge, verte et bleue) servent à la représentation visuelle de la zone urbaine, tandis que la bande « nir » est essentielle au calcul du NDVI, un indice de végétation clé utilisé pour évaluer la présence et la densité de la végétation, comme les arbres urbains.
[7]:
#Create a query object
query = {
'x': lon_range,
'y': lat_range,
'resolution': (-10, 10),
'output_crs':'EPSG:6933',
'time': ('2017', '2022'),
}
#load Sentinel 2 data
ds = dc.load(product=['gm_s2_annual'],
measurements=['red', 'green', 'blue', 'nir'],
**query)
print(ds)
<xarray.Dataset> Size: 546MB
Dimensions: (time: 6, y: 2899, x: 3920)
Coordinates:
* time (time) datetime64[ns] 48B 2017-07-02T11:59:59.999999 ... 202...
* y (y) float64 23kB -5.476e+05 -5.476e+05 ... -5.766e+05
* x (x) float64 31kB 1.463e+06 1.463e+06 ... 1.502e+06 1.502e+06
spatial_ref int32 4B 6933
Data variables:
red (time, y, x) uint16 136MB 1200 1228 1157 1089 ... 347 336 347
green (time, y, x) uint16 136MB 1207 1237 1189 1157 ... 642 588 566
blue (time, y, x) uint16 136MB 1024 1062 1029 992 ... 361 351 356
nir (time, y, x) uint16 136MB 2546 2587 2626 ... 3262 3323 3498
Attributes:
crs: epsg:6933
grid_mapping: spatial_ref
Un géopolygone représente les limites et non la forme réelle, car il est conçu pour représenter l’étendue de l’entité géographique cartographiée, plutôt que sa forme exacte. En d’autres termes, le géopolygone est utilisé pour définir la limite extérieure de la zone d’intérêt, plutôt que les entités et caractéristiques internes.
Il est important de découper les données selon la forme exacte de la zone d’intérêt, car cela permet de garantir que les données utilisées sont pertinentes pour la zone d’étude spécifique. Bien qu’un géopolygone fournisse des informations sur la limite de l’entité géographique représentée, il ne reflète pas nécessairement la forme ou l’étendue exacte de la zone d’intérêt.
[8]:
#Rasterise the area of interest polygon
aoi_raster = xr_rasterize(gdf=geopolygon_gdf, da=ds, crs=ds.crs)
#Mask the dataset to the rasterised area of interest
ds = ds.where(aoi_raster == 1)
# Plot the RGB for a sinle image (first image)
rgb(ds.isel(time=0))
L’indice de végétation par différence normalisée (NDVI) est calculé pour évaluer la santé et la densité de la végétation. Largement utilisé, il quantifie la différence entre la lumière proche infrarouge (NIR) et la lumière rouge réfléchie par la végétation. Cette différence est normalisée afin de réduire les variations atmosphériques et liées aux capteurs. La formule du NDVI est la suivante :
Les valeurs NDVI varient de -1 à 1, où les valeurs les plus élevées (proches de 1) indiquent une végétation saine et dense, et les valeurs inférieures (proches de -1) indiquent des surfaces non végétalisées.
[9]:
# Calculate the chosen vegetation proxy index and add it to the loaded data set
ds = calculate_indices(ds=ds, index='NDVI', satellite_mission='s2', drop=False)
Le seuil NDVI inférieur pour la végétation a été déterminé en analysant la distribution des valeurs NDVI dans l’image cartographiée. L’hypothèse était que les pixels dont le NDVI était inférieur à la médiane moins la distance entre la médiane et le 95e percentile de la distribution étaient des pixels non végétaux (par exemple, des bâtiments, des routes ou des lacs situés en forêt) (Ottosen et al. 2020). Ce seuil NDVI initial permet de distinguer les zones à forte végétation, notamment les prairies, les arbustes et les arbres. Afin de différencier spécifiquement les arbres urbains des autres types de végétation, un seuillage Otsu est appliqué ultérieurement. Cette méthode affine la classification pour distinguer les zones arborées des autres types de végétation.
[10]:
# Loop through each time step
for time_step in ds.time:
# Select data for the current time step
ds_time = ds.sel(time=time_step)
# Compute NDVI statistics for the current time step
ndvi_median = ds_time.NDVI.median(dim=['x', 'y']).values
ndvi_95th = np.nanpercentile(ds_time.NDVI.values, 95)
# Determine the NDVI threshold for the current time step
ndvi_threshold = (ndvi_median - (ndvi_95th - ndvi_median)).item()
# Apply the NDVI threshold and mask out non-tree areas
tree_extent = ds_time.NDVI.where(ds_time.NDVI > ndvi_threshold)
# Plot NDVI for the current time step with thresholded tree extent
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
ds_time.NDVI.plot(ax=axes[0], cmap='RdYlGn', vmin=-1, vmax=1)
tree_extent.plot(ax=axes[1], cmap='YlGn', vmin=0, vmax=1)
# Add titles to each subplot
axes[0].set_title(f'NDVI for {pd.to_datetime(time_step.values).year}')
axes[1].set_title(f'Vegetated area for {pd.to_datetime(time_step.values).year}')
plt.tight_layout()
plt.show()
Afin de différencier spécifiquement les arbres urbains des autres types de végétation, la méthode multi-niveaux d’Otsu est appliquée aux valeurs NDVI. Cette méthode statistique détermine automatiquement les valeurs seuils optimales pour classer les intensités des pixels en plusieurs catégories distinctes. Cette approche est idéale pour segmenter la végétation en différentes classes de densité selon le NDVI, qui varie de -1 à 1, les valeurs les plus élevées indiquant une végétation plus dense.
Un histogramme des valeurs NDVI est généré, et la méthode détermine deux seuils, chacun représenté par des lignes verticales pointillées rouges. Ces seuils divisent les données NDVI en trois catégories :
Végétation basse : Zones avec une végétation minimale, comme les infrastructures urbaines ou les terres arides.
Végétation modérée : Zones à densité de végétation intermédiaire, y compris les arbustes et les prairies.
Végétation haute : Zones à végétation dense, ciblant spécifiquement la couverture arborée, qui est l’objectif principal de la cartographie des arbres urbains.
En séparant les valeurs NDVI dans ces catégories, cette méthode améliore la précision de l’identification des zones couvertes d’arbres, en les distinguant des prairies ou des arbustes, affinant ainsi l’analyse globale des espaces verts urbains.
L’étendue des arbres urbains est visualisée en appliquant le seuil d’Otsu aux données NDVI. Plus précisément, les valeurs NDVI dépassant le deuxième seuil d’Otsu (otsu_thresholds[1]) sont sélectionnées pour représenter les zones à forte densité de végétation, susceptibles d’être des arbres urbains.
[11]:
# Prepare a list to store tree area data
tree_area_data = []
# Loop through each time step
for time_step in ds.time:
# Select data for the current time step
ds_time = ds.sel(time=time_step)
# Compute NDVI statistics for the current time step
ndvi_median = ds_time.NDVI.median(dim=['x', 'y']).values
ndvi_95th = np.nanpercentile(ds_time.NDVI.values, 95)
# Determine the NDVI threshold for the current time step
ndvi_threshold = (ndvi_median - (ndvi_95th - ndvi_median)).item()
# Filter NDVI values above the threshold and flatten them
ndvi_values = ds_time.NDVI.where(ds_time.NDVI > ndvi_threshold).values.flatten()
remaining_ndvi = ndvi_values[~np.isnan(ndvi_values)] # Exclude NaN values
# Compute multi-level Otsu thresholds
otsu_thresholds = threshold_multiotsu(remaining_ndvi, classes=3)
# Create a binary mask for tree extent using the second Otsu threshold
tree_extent = ds_time.NDVI > otsu_thresholds[1]
# Calculate the total area covered by trees for the current time step
# Get the spatial resolution (pixel size) in meters
resolution_x, resolution_y = ds_time.rio.resolution() # May return negative values, take the absolute value
pixel_area_m2 = abs(resolution_x * resolution_y) # Calculate the pixel area in square meters
pixel_area_km2 = pixel_area_m2 / 1_000_000 # Convert to square kilometers
# Count the number of pixels classified as trees
tree_pixels = np.sum(tree_extent)
# Calculate the total area covered by trees in km²
tree_area_km2 = tree_pixels * pixel_area_km2
# Append the result to the list with the corresponding time step
tree_area_data.append({'date': pd.to_datetime(time_step.values).year, 'tree_area_km2': tree_area_km2})
# Plot the histogram and tree extent as before
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# Plot histogram of NDVI values
axes[0].hist(remaining_ndvi, bins=50, edgecolor='k', alpha=0.7)
for threshold in otsu_thresholds:
axes[0].axvline(threshold, color='r', linestyle='--', label=f'Threshold {threshold:.2f}')
axes[0].set_xlabel('NDVI Value')
axes[0].set_ylabel('Frequency')
axes[0].set_title(f'Distribution of NDVI Values with Multi-level Otsu Thresholds ({pd.to_datetime(time_step.values).year})')
axes[0].legend()
# Plot the tree extent
tree_extent.plot(ax=axes[1], cmap='Greens', add_colorbar=False)
axes[1].set_title(f'Tree Extent for {pd.to_datetime(time_step.values).year}')
axes[1].set_xlabel('Longitude')
axes[1].set_ylabel('Latitude')
plt.tight_layout()
plt.show()
Calculer et analyser la superficie du couvert forestier
Détermine la quantité de terre couverte par les arbres dans la zone d’intérêt, fournissant des informations sur l’étendue de la couverture arborée et sa proportion par rapport à la zone d’intérêt globale.
[12]:
# Convert the list of results into a DataFrame for easier analysis
tree_area_df = pd.DataFrame(tree_area_data)
# Now you can visualize the time series of tree area
plt.figure(figsize=(10, 5))
plt.plot(tree_area_df['date'], tree_area_df['tree_area_km2'], marker='o')
plt.title('Time series of tree area covered (km²)')
plt.xlabel('Date')
plt.ylabel('Area covered by trees (km²)')
plt.xticks(rotation=45)
plt.tight_layout()
plt.grid()
plt.show()
Informations Complémentaires
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
[13]:
print(datacube.__version__)
1.8.20
Dernier test :
[14]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[14]:
'2025-01-31'