Cartographie des changements à long terme dans l’étendue des eaux avec WOfS

Mots-clés : données utilisées ; WOfS, eau ; étendue, analyse ; séries chronologiques, 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 bloc-notes montre comment charger, visualiser et analyser le produit « résumé annuel WOfS <https://docs.digitalearthafrica.org/en/latest/data_specs/Landsat_WOfS_specs.html> » pour recueillir des informations sur l’étendue à long terme des masses d’eau. Il complète le bloc-notes « Water_extent_sentinel_2 <Water_extent_sentinel_2.ipynb> » qui se concentre sur les étendues d’eau plus récentes à des intervalles de temps saisonniers.


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

# Force GeoPandas to use Shapely instead of PyGEOS
# In a future release, GeoPandas will switch to using Shapely by default.
import os
os.environ['USE_PYGEOS'] = '0'

import datacube
import numpy as np
import xarray as xr
import seaborn as sns
import geopandas as gpd
import matplotlib.pyplot as plt
from IPython.display import Image
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
from datacube.utils.geometry import Geometry

from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.plotting import display_map, xr_animation
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')

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 ('1990', '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 334a9d1eb26a490a841b12da9c9892bb 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.

[3]:
# 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 = '1990'
end_year = '2021'

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.

[4]:
display_map(lon_range, lat_range)
[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Télécharger les résumés annuels du WOfS

[5]:
#Create a query object
query = {
    'x': lon_range,
    'y': lat_range,
    'resolution': (-30, 30),
    'output_crs':'EPSG:6933',
    'time': (start_year, end_year),
}

#load wofs
ds = dc.load(product="wofs_ls_summary_annual",
             **query)

print(ds)
<xarray.Dataset> Size: 13MB
Dimensions:      (time: 32, y: 255, x: 194)
Coordinates:
  * time         (time) datetime64[ns] 256B 1990-07-02T11:59:59.999999 ... 20...
  * y            (y) float64 2kB -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:
    count_wet    (time, y, x) int16 3MB 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
    count_clear  (time, y, x) int16 3MB 3 3 3 3 3 3 3 3 ... 21 21 18 19 19 19 20
    frequency    (time, y, x) float32 6MB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
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.

[6]:
#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)

Graphique à facettes, un sous-ensemble des résumés annuels du WOfS

[7]:
ds.isel(time=[5,10,15,20,25]).frequency.plot(col='time', col_wrap=5, cmap=sns.color_palette("mako_r", as_cmap=True));
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_WOfS_18_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.

[8]:
out_path = 'annual_water_frequency.gif'

xr_animation(ds=ds,
             output_path=out_path,
             interval=400,
             bands=['frequency'],
             show_text='WOfS Annual Summary',
             show_date = '%Y',
             width_pixels=300,
             annotation_kwargs={'fontsize': 15},
             imshow_kwargs={'cmap': sns.color_palette("mako_r", as_cmap=True), 'vmin': 0.0, 'vmax': 0.9},
             colorbar_kwargs={'colors': 'black'},
             show_colorbar=False)

# Plot animated gif
plt.close()
Image(filename=out_path)
Exporting animation to annual_water_frequency.gif
[8]:
<IPython.core.display.Image object>

Calculer la superficie annuelle de l’étendue d’eau

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.

[9]:
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

Seuil de fréquence annuelle WOfS pour classer eau/non-eau

Calcule la surface des pixels classés comme eau (si « ds.frequency » est > 0,20, alors le pixel sera considéré comme une eau libre régulière pendant l’année)

[10]:
water_threshold = 0.20
[11]:
#threshold
water_extent = ds.frequency > water_threshold

#calculate area
ds_valid_water_area = water_extent.sum(dim=['x', 'y']) * area_per_pixel

Tracer la superficie annuelle des eaux libres

[12]:
plt.figure(figsize=(18, 4))
ds_valid_water_area.plot(marker='o', color='#9467bd')
plt.title(f'Observed Annual 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_WOfS_27_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 ».

[13]:
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(['1992-07-01T23:59:59.999999000', '2021-07-02T11:59:59.999999000'],
      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.

[14]:
water_extent.sel(time=time_xr).plot.imshow(col="time", col_wrap=2, figsize=(14, 6));
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_WOfS_31_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

[15]:
baseline_time = '2019'
analysis_time = '2021'

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

Traçage

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

[16]:
compare = water_extent.sel(time=[baseline_ds.time.values, analysis_ds.time.values])

compare.plot(col="time",col_wrap=2,figsize=(10, 5), cmap='viridis', add_colorbar=False);
../../../_images/sandbox_notebooks_Real_world_examples_Water_extent_WOfS_35_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

[17]:
analyse_total_value = compare.isel(time=1).astype(int)
change = analyse_total_value - compare.isel(time=0).astype(int)

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

[18]:
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

[19]:
water_appeared_color = "Green"
water_disappeared_color = "Yellow"
stable_color = "Blue"
land_color = "Brown"

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

compare[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_WOfS_41_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

[20]:
print(datacube.__version__)
1.8.19

Dernier test :

[21]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[21]:
'2024-11-05'