Déterminer l’étendue saisonnière des plans d’eau avec Sentinel-2

  • Produits utilisés : s2_l2a,

Mots clés : données utilisées ; sentinelle-2, eau ; étendue, analyse ; séries chronologiques, indice de bande ; MNDWI, visualisation ; animation

Aperçu

Les Nations Unies ont prescrit 17 « Objectifs de développement durable » (ODD). Ce carnet tente de suivre l’indicateur ODD 6.6.1 - changement dans l’étendue des écosystèmes liés à l’eau. L’indicateur 6.6.1 comporte 4 sous-indicateurs : > i. L’étendue spatiale des écosystèmes liés à l’eau > ii. La quantité d’eau contenue dans ces écosystèmes > iii. La qualité de l’eau dans ces écosystèmes > iv. La santé ou l’état de ces écosystèmes

Ce carnet se concentre principalement sur le premier sous-indicateur : l’étendue spatiale.

Description

Le cahier montre comment :

  1. Charger des données satellite sur le plan d’eau d’intérêt

  2. Calculer l’indice d’eau MNDWI

  3. Rééchantillonner la série chronologique du MNDWI en fonction des médianes saisonnières

  4. Générer une animation de la série chronologique de l’étendue de l’eau

  5. Calculer et tracer une série chronologique de l’étendue saisonnière des eaux (en kilomètres carrés)

  6. Trouvez les étendues d’eau minimales et maximales dans la série chronologique et tracez-les.

  7. Comparez deux périodes de temps désignées et indiquez où l’étendue de la masse d’eau a changé.

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 matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import geopandas as gpd
from IPython.display import Image
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

from odc.geo.geom import Geometry
from deafrica_tools.datahandling import load_ard
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.plotting import display_map, xr_animation
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.spatial import xr_rasterize
from deafrica_tools.areaofinterest import define_area

Se connecter au datacube

Activez la base de données Datacube, qui fournit des fonctionnalités pour le chargement et l’affichage des données d’observation de la Terre stockées.

[2]:

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

Configurer un cluster Dask

Dask peut être utilisé pour mieux gérer l’utilisation de la mémoire et effectuer l’analyse en parallèle. Pour une introduction à l’utilisation de Dask avec Digital Earth Africa, consultez le Dask notebook.

Remarque : nous vous recommandons d’ouvrir la fenêtre de traitement Dask pour afficher les différents calculs en cours d’exécution ; pour ce faire, consultez la section Tableau de bord Dask en Afrique de l’Ouest du Dask notebook.

Pour activer Dask, configurez le cluster de calcul local à l’aide de la cellule ci-dessous.

[3]:

create_local_dask_cluster()

Client

Client-880a55b6-476f-11f0-8241-a22f78693c99

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /user/nanaboamah89@gmail.com/proxy/8787/status

Cluster Info

Paramètres d’analyse

La cellule suivante définit les paramètres qui définissent la zone d’intérêt et la durée de réalisation de l’analyse.

Les paramètres sont :

  • « lat » : la latitude centrale à analyser (par exemple 10,338).

  • « lon » : la longitude centrale à analyser (par exemple -1,055).

  • lat_buffer : Le nombre de degrés à charger autour de la latitude centrale.

  • lon_buffer : Le nombre de degrés à charger autour de la longitude centrale.

  • start_year et end_year : la plage de dates à analyser (par exemple ('2017', '2020').

Sélectionnez 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. Cette méthode nécessite que vous saisissiez la latitude centrale, la longitude centrale et la valeur de la zone tampon en degrés carrés autour du point central que vous souhaitez analyser. 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 7a62516fca0a46f9b342ab576098dba5 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 bloc-notes pour la première fois, conservez les paramètres par défaut ci-dessous. Cela montrera comment fonctionne l’analyse et fournira des résultats significatifs. L’exemple couvre une partie du lac Sulunga. Tanzanie.

[4]:

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-5.9460, lon=35.5188 , buffer=0.03)

# 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])


# Define the start year and end year
start_year = '2017'
end_year = '2021-05'

Visualisez la zone d’intérêt sur une carte interactive

La cellule suivante affichera la zone sélectionnée sur une carte interactive. La bordure rouge représente la zone d’intérêt de l’étude. Effectuez un zoom avant et arrière pour mieux comprendre la zone d’intérêt. En cliquant n’importe où sur la carte, les coordonnées de latitude et de longitude du point cliqué s’affichent.

[5]:

display_map(lon_range, lat_range)

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

Charger des données satellite masquées par des nuages

Le code ci-dessous va créer un dictionnaire de requêtes pour notre région d’intérêt, puis charger les données du satellite Sentinel-2. Pour plus d’informations sur le chargement des données, consultez le bloc-notes Loading data notebook.

[6]:

#Create a query object
query = {
    'x': lon_range,
    'y': lat_range,
    'resolution': (-20, 20),
    'output_crs':'EPSG:6933',
    'time': (start_year, end_year),
    'dask_chunks':{'time':1,'x':500,'y':500}
}

#load Sentinel 2 data
ds = load_ard(dc=dc,
             products=['s2_l2a'],
             measurements=['green','swir_1'],
             mask_filters=[("opening", 3),("dilation", 2)], #improve cloud mask
             group_by='solar_day',
             **query)

print(ds)

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Applying morphological filters to pq mask [('opening', 3), ('dilation', 2)]
Applying pixel quality/cloud mask
Returning 174 time steps as a dask array
<xarray.Dataset> Size: 154MB
Dimensions:      (time: 174, y: 382, x: 290)
Coordinates:
  * time         (time) datetime64[ns] 1kB 2019-01-04T08:01:41 ... 2021-05-28...
  * y            (y) float64 3kB -7.534e+05 -7.534e+05 ... -7.61e+05 -7.61e+05
  * x            (x) float64 2kB 3.424e+06 3.424e+06 ... 3.43e+06 3.43e+06
    spatial_ref  int32 4B 6933
Data variables:
    green        (time, y, x) float32 77MB dask.array<chunksize=(1, 382, 290), meta=np.ndarray>
    swir_1       (time, y, x) float32 77MB dask.array<chunksize=(1, 382, 290), meta=np.ndarray>
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Découpez les ensembles de données selon la forme de la zone d’intérêt

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.

[7]:

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

Calculer l’indice d’eau MNDWI

[8]:

# Calculate the chosen vegetation proxy index and add it to the loaded data set
ds = calculate_indices(ds=ds, index='MNDWI', satellite_mission='s2', drop=True)

Dropping bands ['green', 'swir_1']

Rééchantillonner les séries chronologiques

En raison de nombreux facteurs (par exemple, des nuages obscurcissant la région, une couverture nuageuse manquante dans la couche fmask), les données seront incomplètes et bruyantes. Ici, nous allons rééchantillonner les données pour nous assurer que nous travaillons avec une série chronologique cohérente.

Pour ce faire, nous rééchantillonnons les données en fonction des pas de temps saisonniers à l’aide de médianes.

Ces calculs prendront plusieurs minutes car nous exécuterons .compute(), déclenchant toutes les tâches que nous avons planifiées ci-dessus et mettant les tableaux en mémoire.

[9]:

%%time
sample_frequency="QS-DEC"  # quarterly starting in DEC, i.e. seasonal

#resample using medians
print('calculating MNDWI seasonal medians...')
mndwi = ds['MNDWI'].resample(time=sample_frequency).median().compute()

calculating MNDWI seasonal medians...

/opt/venv/lib/python3.12/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dest = _reproject(
/opt/venv/lib/python3.12/site-packages/dask/utils.py:78: RuntimeWarning: All-NaN slice encountered
  return func(*args, **kwargs)

CPU times: user 4.49 s, sys: 363 ms, total: 4.85 s
Wall time: 1min 18s

Graphique à facettes des pas de temps MNDWI

[10]:

mndwi.plot(col='time', col_wrap=4, cmap='RdBu', vmax=1, vmin=-1);
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_sentinel_2_24_0.png

Animation de séries chronologiques

Dans la cellule suivante, nous traçons l’ensemble de données que nous avons chargé ci-dessus sous forme de GIF d’animation, en utilisant la fonction `xr_animation <../Frequently_used_code/Animated_timeseries.ipynb>`__. Le output_path sera enregistré dans le répertoire où se trouve le script et vous pouvez modifier les noms pour éviter l’écrasement des fichiers.

[11]:

out_path = 'water_extent.gif'

xr_animation(ds=mndwi.to_dataset(name='MNDWI'),
             output_path=out_path,
             bands = ['MNDWI'],
             show_text = 'Seasonal MNDWI',
             interval=500,
             width_pixels=300,
             show_colorbar=True,
             imshow_kwargs={'cmap':'RdBu','vmin': -0.5, 'vmax': 0.5},
             colorbar_kwargs={'colors': 'black'}
            )

# Plot animated gif
plt.close()
Image(filename=out_path)

Exporting animation to water_extent.gif

[11]:

<IPython.core.display.Image object>

Calculer la surface par pixel

Le nombre de pixels peut être utilisé pour la surface du plan d’eau si la surface en pixels est connue. Exécutez la cellule suivante pour générer les constantes nécessaires à la réalisation de cette conversion.

[12]:

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

Calculer l’étendue de l’eau

Calcule la surface des pixels classés comme eau (si MNDWI est > 0, alors eau)

[13]:

water = mndwi.where(mndwi > 0, np.nan)
area_ds = water.where(np.isnan(water),1)
ds_valid_water_area = area_ds.sum(dim=['x', 'y']) * area_per_pixel

Tracer des séries chronologiques saisonnières de « l’année de début » à « l’année de fin »

[14]:

plt.figure(figsize=(18, 4))
ds_valid_water_area.plot(marker='o', color='#9467bd')
plt.title(f'Observed Seasonal Area of Water from {start_year} to {end_year}')
plt.xlabel('Dates')
plt.ylabel('Waterbody area (km$^2$)')
plt.tight_layout()
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_sentinel_2_32_0.png

Déterminer l’étendue d’eau minimale et maximale

La cellule suivante extrait l’étendue minimale et maximale de l’eau de l’ensemble de données à l’aide des fonctions « min » et « max », nous ajoutons ensuite les dates à un « xarray.DataArray ».

[15]:

min_water_area_date, max_water_area_date =  min(ds_valid_water_area), max(ds_valid_water_area)
time_xr = xr.DataArray([min_water_area_date.time.values, max_water_area_date.time.values], dims=["time"])

print(time_xr)

<xarray.DataArray (time: 2)> Size: 16B
array(['2019-09-01T00:00:00.000000000', '2021-03-01T00:00:00.000000000'],
      dtype='datetime64[ns]')
Dimensions without coordinates: time

Tracez les dates auxquelles se produisent les étendues d’eau minimales et maximales

Tracez le pixel classé de l’eau pour les deux dates où nous avons l’étendue minimale et maximale de l’eau de surface.

[16]:

area_ds.sel(time=time_xr).plot.imshow(col="time", col_wrap=2, figsize=(14, 6));
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_sentinel_2_36_0.png

Comparer deux périodes

Les cellules suivantes déterminent l’étendue maximale de l’eau pour deux années différentes.

  • baseline_year : l’année de référence pour l’analyse

  • analysis_year : l’année à comparer à l’année de référence

[17]:

baseline_time = '2019-03-01'
analysis_time = '2020-03-01'

baseline_ds, analysis_ds = ds_valid_water_area.sel(time=baseline_time, method ='nearest'),ds_valid_water_area.sel(time=analysis_time, method ='nearest')

Un nouveau dataArray est créé pour stocker la nouvelle date de « l’étendue d’eau maximale » pour les deux années

[18]:

time_xr = xr.DataArray([baseline_ds.time.values, analysis_ds.time.values], dims=["time"])

Traçage

Tracez l’étendue de l’eau du produit MNDWI pour les deux périodes choisies.

[19]:

area_ds.sel(time=time_xr).plot(col="time", col_wrap=2, robust=True, figsize=(10, 5), cmap='viridis', add_colorbar=False);
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_sentinel_2_42_0.png

Calcul de la variation pour les deux périodes désignées

Les cellules ci-dessous calculent la quantité d’eau gagnée, perdue et stable pour les deux périodes

[20]:

# The two period Extract the two periods(Baseline and analysis) dataset from
ds_selected = area_ds.where(area_ds == 1, 0).sel(time=time_xr)

analyse_total_value = ds_selected[1]
change = analyse_total_value - ds_selected[0]

water_appeared = change.where(change == 1)
permanent_water = change.where((change == 0) & (analyse_total_value == 1))
permanent_land = change.where((change == 0) & (analyse_total_value == 0))
water_disappeared = change.where(change == -1)

La cellule ci-dessous calcule la superficie de l’étendue de l’eau pour la perte d’eau, le gain d’eau, l’eau permanente et le sol

[21]:

total_area = analyse_total_value.count().values * area_per_pixel
water_apperaed_area = water_appeared.count().values * area_per_pixel
permanent_water_area = permanent_water.count().values * area_per_pixel
water_disappeared_area = water_disappeared.count().values * area_per_pixel

Traçage

Les variables de l’eau sont tracées pour visualiser le résultat

[22]:

water_appeared_color = "Green"
water_disappeared_color = "Yellow"
stable_color = "Blue"
land_color = "Brown"

fig, ax = plt.subplots(1, 1, figsize=(10, 10))

ds_selected[1].plot.imshow(cmap="Pastel1",
                                       add_colorbar=False,
                                       add_labels=False,
                                       ax=ax)
water_appeared.plot.imshow(
    cmap=ListedColormap([water_appeared_color]),
    add_colorbar=False,
    add_labels=False,
    ax=ax,
)
water_disappeared.plot.imshow(
    cmap=ListedColormap([water_disappeared_color]),
    add_colorbar=False,
    add_labels=False,
    ax=ax,
)
permanent_water.plot.imshow(cmap=ListedColormap([stable_color]),
                            add_colorbar=False,
                            add_labels=False,
                            ax=ax)

plt.legend(
    [
        Patch(facecolor=stable_color),
        Patch(facecolor=water_disappeared_color),
        Patch(facecolor=water_appeared_color),
        Patch(facecolor=land_color),
    ],
    [
        f"Water to Water {round(permanent_water_area, 2)} km2",
        f"Water to No Water {round(water_disappeared_area, 2)} km2",
        f"No Water to Water: {round(water_apperaed_area, 2)} km2",
    ],
    loc="lower left",
)

plt.title("Change in water extent: " + baseline_time + " to " + analysis_time);
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_sentinel_2_48_0.png

Prochaines étapes

Revenez à la section « Paramètres d’analyse », modifiez certaines valeurs (par exemple « latitude », « longitude », « start_year », « end_year ») et relancez l’analyse. Vous pouvez utiliser la carte interactive dans la section « Afficher l’emplacement sélectionné » pour trouver de nouvelles valeurs de latitude et de longitude centrales en effectuant un panoramique et un zoom, puis en cliquant sur la zone pour laquelle vous souhaitez extraire les valeurs d’emplacement. Vous pouvez également utiliser Google Maps pour rechercher un emplacement que vous connaissez, puis renvoyer les valeurs de latitude et de longitude en cliquant sur la carte.

Modifiez également l’année dans la section « Comparer deux périodes : une ligne de base et une analyse » (par exemple, « base_year », « analyse_year ») et réexécutez l’analyse.

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

[23]:

print(datacube.__version__)

1.8.20

Dernier test :

[24]:

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

[24]:

'2025-06-12'