Animations à long terme sur les plans d’eau
Produits utilisés: `"Waterbodies" de DE Africa<https://docs.digitalearthafrica.org/en/latest/data_specs/Waterbodies_specs.html>`__; DE Africa WOfS; GeoMADs; CHIRPS
Mots clés: données utilisées; Waterbodies (pans d’eau), données utilisées; wofs_ls_summary_annual, données utilisées; CHIRPS, données utilisées; géomédiane sentinel-2, données utilisées; géomédiane landsat 5 et 7, données utilisées; géomédiane landsat 8 et 9
Aperçu
Le service de surveillance des plans d’eau continentaux de Digital Earth Africa identifie plus de 700 000 plans d’eau à partir de plus de trois décennies d’observations par satellite. Ce service cartographie les plans d’eau persistants et saisonniers et l’évolution de leur surface d’eau au fil du temps. Les plans d’eau cartographiés peuvent inclure, sans s’y limiter, des lacs, des étangs, des réservoirs artificiels, des zones humides et des segments de certains systèmes fluviaux. Pour plus d’informations sur le service de surveillance des plans d’eau, reportez-vous au « carnet de données » <../Datasets/Waterbodies.ipynb>`__.
Il est souvent judicieux de visualiser les changements dans les masses d’eau au fil du temps à l’aide d’animations. Les animations sont une forme d’analyse exploratoire qui peut aider à montrer s’il existe une tendance générale dans l’étendue de l’eau ou illustrer des modèles saisonniers répétitifs. Ce bloc-notes montre comment générer des animations pour la composante de séries chronologiques des masses d’eau, ainsi que des données récapitulatives provenant des observations de l’eau depuis l’espace, des données sur les précipitations et des images en vraies couleurs.
Dans un premier temps, une animation montrant les changements annuels de l’étendue des eaux est produite pour un lac au Kenya, dont l’étendue des eaux tend à augmenter. Ensuite, une animation est générée pour le lac Tchad au Nigéria, dont l’étendue des eaux tend à diminuer.
Avertissement : DE Africa Waterbodies Surface Area Change mesure la surface humide des plans d’eau telle qu’estimée à partir des satellites. Ce produit ne mesure pas la profondeur, le volume, la fonction du plan d’eau, ni la source de l’eau.
Description
Ce notebook démontre la génération de deux animations permettant d’inspecter les changements dans l’étendue de l’eau sur plusieurs années, en s’appuyant sur le service « DE Africa Waterbodies » <https://docs.digitalearthafrica.org/en/latest/data_specs/Waterbodies_specs.html>.
Les mesures prises sont les suivantes :
Chargement et préparation de séries chronologiques annuelles, de WOfS et de données en couleurs réelles.
Création d’une animation sur les changements annuels de l’étendue d’eau.
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]:
import matplotlib.pyplot as plt
import datacube
import matplotlib.dates as mdates
import numpy as np
import matplotlib.animation as animation
from deafrica_tools.plotting import display_map
from IPython.display import HTML
from deafrica_tools.spatial import xr_rasterize
from deafrica_tools.waterbodies import (
get_geohashes,
get_waterbodies,
get_waterbody,
get_time_series,
display_time_series,
)
[2]:
dc = datacube.Datacube(app="Waterbody-anim")
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 (en degrés)
La zone par défaut est le lac Baringo au Kenya qui a connu « une augmentation de son étendue au cours des dernières années <https://www.theguardian.com/world/2022/mar/17/kenya-quiet-slide-underwater-great-rift-valley-lakes-east-africa-flooding> »__ et a « déjà été étudié avec Digital Earth Africa <https://www.digitalearthafrica.org/media-center/blog/rising-lakes-rift-valley-kenya> »__.
[3]:
# Set the central latitude and longitude
lat = 0.62
lon = 36.08
# Set the buffer to load around the central coordinates
buffer = 0.3
# Compute the bounding box coordinates
xlim = (lon - buffer, lon + buffer)
ylim = (lat + buffer, lat - buffer)
# Preview area on a map
display_map(xlim, ylim)
[3]:
Obtenir des données
Le module deafrica_waterbodies vous permet d’interroger les plans d’eau par localisation ou par géohash. Pour plus d’informations, consultez la documentation sur le service « waterbodies »<https://docs.digitalearthafrica.org/en/latest/data_specs/Waterbodies_specs.html>`__ et le notebook des jeux de données sur les plans d’eau <../Datasets/Waterbodies.ipynb>`__.
Répertorier les polygones de plans d’eau dans une zone
Nous pouvons obtenir une liste de polygones de plans d’eau à l’intérieur d’une boîte englobante de coordonnées en utilisant « get_waterbodies ».
[4]:
# Create a bounding box from study area coordinates
bbox = (xlim[0], ylim[1], xlim[1], ylim[0])
# Select all water bodies located within the bounding box
polygons = get_waterbodies(bbox, crs="EPSG:4326")
Le cadre de données géographiques renvoyé inclut tous les plans d’eau qui se trouvent dans la zone de délimitation. Cet ensemble de données contient des métadonnées pour chaque plan d’eau de l’ensemble de données, y compris l’ID, l’UID, le WB_UID, la superficie, le périmètre et la série chronologique. Consultez la documentation sur l’étendue historique des plans d’eau <https://docs.digitalearthafrica.org/en/latest/data_specs/Waterbodies_specs.html#Waterbodies-Historical-Extent>`__ pour obtenir une description de chaque attribut.
Explorez et sélectionnez le polygone
Une fois les polygones du plan d’eau en mémoire, vous pouvez les tracer directement ou les explorer dans une fenêtre interactive.
[5]:
# Explore the waterbody polygons located within the bounding box
polygons.explore()
[5]:
Sélectionnez le géohash
En explorant les polygones ci-dessus, nous pouvons voir que notre plan d’eau d’intérêt, le lac Baringo, est représenté par un seul grand polygone. Il existe de nombreux polygones plus petits à proximité ou immédiatement autour du lac. Le géohachage du polygone principal (affiché comme « uid » lorsque vous passez la souris dessus) est « sb1ekyxw88 ».
Obtenir des données de séries chronologiques
[6]:
selected_waterbody_geohash = "sb1ekyxw88"
selected_waterbody = get_waterbody(selected_waterbody_geohash)
Obtenir la série chronologique de la surface mouillée pour le plan d’eau sélectionné
Pour tout géohachage ou polygone donné, nous pouvons également utiliser la fonction get_time_series() pour obtenir diverses mesures de la surface de la masse d’eau au fil du temps. Consultez la documentation Waterbodies Historical Extent pour obtenir une description des différentes mesures de surface.
La fonction calcule également une médiane mobile du pourcentage de surface mouillée de la masse d’eau. Elle permet de visualiser la tendance générale du pourcentage de surface mouillée. La médiane mobile utilise les trois dernières observations pour déterminer la médiane à une date donnée.
[7]:
# Get time series for the selected water body
selected_waterbody_timeseries = get_time_series(waterbody=selected_waterbody)
selected_waterbody_timeseries.head()
[7]:
| area_wet_m2 | percent_wet | area_dry_m2 | percent_dry | area_invalid_m2 | percent_invalid | area_observed_m2 | percent_observed | percent_wet_rolling_median | |
|---|---|---|---|---|---|---|---|---|---|
| date | |||||||||
| 1984-05-30 | 142439400.0 | 62.88 | 79190100.0 | 34.96 | 4910400.0 | 2.17 | 226539900.0 | 100.0 | NaN |
| 1984-06-15 | 147272400.0 | 65.01 | 79223400.0 | 34.97 | 44100.0 | 0.02 | 226539900.0 | 100.0 | NaN |
| 1984-07-01 | 147078900.0 | 64.92 | 79350300.0 | 35.03 | 110700.0 | 0.05 | 226539900.0 | 100.0 | 64.92 |
| 1984-07-17 | 147500100.0 | 65.11 | 78994800.0 | 34.87 | 45000.0 | 0.02 | 226539900.0 | 100.0 | 65.01 |
| 1984-12-08 | 132145200.0 | 58.33 | 83509200.0 | 36.86 | 10885500.0 | 4.81 | 226539900.0 | 100.0 | 64.92 |
Afficher la série chronologique de la surface mouillée pour le plan d’eau sélectionné
Après avoir chargé la série chronologique du plan d’eau, nous pouvons utiliser la fonction « display_time_series() » pour créer une visualisation interactive de la série chronologique.
La visualisation montre le pourcentage invalide et le pourcentage humide.
La disponibilité limitée des données avant 2009 est évidente, en raison de la combinaison de la présence de capteurs (Landsat-8 a été lancé en 2013 et Landsat-9 en 2021) et de la qualité des données.
L’augmentation de l’étendue d’eau de 2010 à aujourd’hui est clairement observable dans les séries chronologiques.
[8]:
display_time_series(selected_waterbody_timeseries)
Charger les GeoMADs en vraies couleurs
Les images composites annuelles géomédianes sans nuages sont une excellente source de données en vraies couleurs à présenter parallèlement aux changements annuels de l’étendue des eaux. La cellule ci-dessous charge les bandes de vraies couleurs (rouge, vert et bleu) de trois produits annuels GeoMAD : Landsat-5 et Landsat-7 combinés, Landsat-8 et Landsat-8 et Landsat-9 combinés. Cela produit une série chronologique soignée d’images annuelles pour notre période d’intérêt, commençant en 2009 en raison de la disponibilité des données.
[9]:
lat_range = (selected_waterbody.total_bounds[1], selected_waterbody.total_bounds[3])
lon_range = (selected_waterbody.total_bounds[0], selected_waterbody.total_bounds[2])
ds = dc.load(
product=["gm_ls5_ls7_annual", "gm_ls8_annual", "gm_ls8_ls9_annual"],
measurements=["red", "green", "blue"],
x=lon_range,
y=lat_range,
resolution=(-30, 30),
dask_chunks={"time": 1, "x": 3000, "y": 3000},
time=("2009-01-01", "2022-12-31"),
).compute()
Charger le résumé annuel des observations de l’eau depuis l’espace (WOfS)
Notre première animation superposera WOfS sur les images en vraies couleurs. Comme nous produisons une animation annuelle, le résumé annuel WOfS est un produit logique à utiliser.
[10]:
query = {
"x": lon_range,
"y": lat_range,
"time": ("2009-01-01", "2022-12-31"),
"resolution": (-30, 30),
}
wofs = dc.load(product="wofs_ls_summary_annual", output_crs="epsg:6933", **query)
Masque WOfS pour plan d’eau
Les données WOfS doivent être masquées sur toute la surface du plan d’eau pour produire une superposition nette. Le polygone du plan d’eau est d’abord rastérisé puis appliqué comme masque. La couche WOfS résultante est ensuite tracée pour vérifier que le masquage a fonctionné comme prévu.
[11]:
fig, ax = plt.subplots()
ax.set_aspect('equal')
mask = xr_rasterize(selected_waterbody, wofs)
clipped_wofs = wofs.where(mask)
clipped_wofs.frequency.isel(time=0).plot()
[11]:
<matplotlib.collections.QuadMesh at 0x7f3249232ce0>
Découper et rééchantillonner des séries chronologiques de plans d’eau
La série temporelle doit avoir la même longueur que les jeux de données raster pour produire une animation intuitive. Cela signifie que nous devons rééchantillonner la série temporelle en pas de temps annuels et la limiter à la période comprise entre 2009 et 2022.
[12]:
tsw = (
selected_waterbody_timeseries.percent_wet.loc["2009-01-01":"2022-12-31"]
.resample("1Y")
.mean()
.interpolate()
)
tsw
[12]:
date
2009-12-31 58.107500
2010-12-31 58.945000
2011-12-31 61.466250
2012-12-31 70.880000
2013-12-31 77.166000
2014-12-31 82.610000
2015-12-31 82.553125
2016-12-31 82.785882
2017-12-31 80.502353
2018-12-31 83.196000
2019-12-31 83.286923
2020-12-31 90.783636
2021-12-31 95.685000
2022-12-31 96.086667
Freq: A-DEC, Name: percent_wet, dtype: float64
L’impression de la longueur de chaque série temporelle à inclure dans l’animation garantit qu’elles ont des longueurs égales.
[13]:
print(len(ds.time))
print(len(clipped_wofs.time))
print(len(tsw))
14
14
14
Générer une animation annuelle
Le code de la cellule ci-dessous produit une animation à des pas de temps annuels. Les premières sections de la cellule définissent la structure du tracé, puis les données de chaque axe sont définies. Une fonction « update » détermine ensuite les règles de mise à jour de chaque série temporelle à des pas de temps annuels.
La sortie peut être enregistrée sous forme de fichier « .gif » en décommentant (en supprimant le « # ») la ligne « ani.save() ».
L’animation montre une augmentation marquée de l’étendue d’eau entre 2010 et 2015, puis une stabilisation jusqu’en 2020 avant une forte augmentation. Ces deux fortes augmentations de l’étendue d’eau sont associées à un verdissement du paysage en arrière-plan.
[ ]:
# create a figure and axes
fig = plt.figure(figsize=(10, 6))
ax1 = plt.subplot(122)
ax2 = plt.subplot(121)
ax1.set_title("Waterbodies Timeseries")
ax1.set_xlabel("Date")
ax1.set_ylabel("Surface water extent (%)")
years = mdates.YearLocator(2)
ax1.xaxis.set_major_locator(years)
ax1.tick_params(axis="x", labelrotation=45)
bands = ["red", "green", "blue"]
cax = (
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
/ np.max(
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
)
).plot.imshow(rgb="variable", animated=True, robust=True, ax=ax2)
ax2.set_aspect('equal')
dax = clipped_wofs.frequency[0, :, :].plot(cmap="Blues", ax=ax2)
selected_waterbody.plot(ax=ax2, edgecolor="black", color="none")
def update(num, x, y, line):
dax.set_array(clipped_wofs.frequency[num, :, :])
cax.set_array(
(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
/ np.max(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
)
ax2.set_title("Time = " + str(clipped_wofs.frequency.coords["time"].values[(num)])[:12])
line.set_data(x[:num], y[:num])
return (line,)
x = tsw.index
y = tsw.values
(line,) = ax1.plot(x, y)
plt.tight_layout()
ani = animation.FuncAnimation(
fig, update, len(ds.time), fargs=[x, y, line], interval=800, blit=True
)
# ani.save('LakeBaringo_Kenya_ts.gif')
plt.close()
HTML(ani.to_html5_video())
Obtenir des données sur le lac Tchad
Prenons un autre exemple pour examiner une tendance à l’assèchement. Le lac Tchad au Nigéria est considéré comme « asséché de 90 % depuis les années 1960 » <https://www.un.org/africarenewal/magazine/december-2019-march-2020/drying-lake-chad-basin-gives-rise-crisis>. Bien sûr, le service WOfS (dérivé de Landsat) ne remonte pas aussi loin, mais nous pouvons examiner les tendances récentes.
Les cellules ci-dessous suivent le même processus que celui présenté pour Baringo, avec moins de sorties pour plus de concision.
Tout d’abord, nous sélectionnons le polygone. La nature de la zone signifie qu’il existe de nombreux polygones de plans d’eau qui composent la plus grande superficie du lac Tchad. Nous sélectionnons le plus grand polygone qui couvre le plan d’eau le plus grand et le plus permanent, à savoir « s66615ywe3 ».
[ ]:
# Set the central latitude and longitude
lat = 13.63
lon = 13.77
# Set the buffer to load around the central coordinates
buffer = 0.9
# Compute the bounding box coordinates
xlim = (lon - buffer, lon + buffer)
ylim = (lat + buffer, lat - buffer)
# Create a bounding box from study area coordinates
bbox = (xlim[0], ylim[1], xlim[1], ylim[0])
# Select all water bodies located within the bounding box
polygons = get_waterbodies(bbox, crs="EPSG:4326")
# Explore the waterbody polygons located within the bounding box
polygons.explore()
Inspecter les séries chronologiques du lac Tchad
Nous pouvons constater que des données cohérentes sont disponibles depuis fin 2013. Nous examinerons la période 2013-2023 pour visualiser les tendances.
[ ]:
selected_waterbody_geohash = "s66615ywe3"
selected_waterbody = get_waterbody(selected_waterbody_geohash)
# Get time series for the selected water body
selected_waterbody_timeseries = get_time_series(waterbody=selected_waterbody)
display_time_series(selected_waterbody_timeseries)
Maintenant que nous avons la période de visualisation, nous allons charger et traiter WOfS et les séries chronologiques des plans d’eau. L’impression de cette cellule garantit que chaque objet a la même longueur en années.
[ ]:
lat_range = (selected_waterbody.total_bounds[1], selected_waterbody.total_bounds[3])
lon_range = (selected_waterbody.total_bounds[0], selected_waterbody.total_bounds[2])
ds = dc.load(
product=["gm_ls5_ls7_annual", "gm_ls8_annual", "gm_ls8_ls9_annual"],
measurements=["red", "green", "blue"],
x=lon_range,
y=lat_range,
resolution=(-30, 30),
dask_chunks={"time": 1, "x": 3000, "y": 3000},
time=("2013-01-01", "2023-12-31"),
).compute()
query = {
"x": lon_range,
"y": lat_range,
"time": ("2013-01-01", "2023-12-31"),
"resolution": (-30, 30),
}
wofs = dc.load(product="wofs_ls_summary_annual", output_crs="epsg:6933", **query)
plt.figure(figsize=(6,6))
mask = xr_rasterize(selected_waterbody, wofs)
clipped_wofs = wofs.where(mask)
tsw = (
selected_waterbody_timeseries.percent_wet.loc["2009-01-01":"2023-12-31"]
.resample("1Y")
.mean()
.interpolate()
)
print(len(ds.time))
print(len(clipped_wofs.time))
print(len(tsw))
Générer une animation sur le lac Tchad
Enfin, nous pouvons exécuter le code d’animation pour le lac Tchad. Le résultat montre une tendance à la baisse, même s’il convient de noter que l’axe se situe entre 81 % et 88 %.
La sortie peut être sauvegardée en supprimant le # de la commande ani.save('LakeChad_Nigeria_ts.gif').
[ ]:
# create a figure and axes
fig = plt.figure(figsize=(10, 6))
ax1 = plt.subplot(122)
ax2 = plt.subplot(121)
ax1.set_title("Waterbodies Timeseries")
ax1.set_xlabel("Date")
ax1.set_ylabel("Surface water extent (%)")
years = mdates.YearLocator(2)
ax1.xaxis.set_major_locator(years)
ax1.tick_params(axis="x", labelrotation=45)
bands = ["red", "green", "blue"]
cax = (
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
/ np.max(
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
)
).plot.imshow(rgb="variable", animated=True, robust=True, ax=ax2)
ax2.set_aspect('equal')
dax = clipped_wofs.frequency[0, :, :].plot(cmap="Blues", ax=ax2)
selected_waterbody.plot(ax=ax2, edgecolor="black", color="none")
def update(num, x, y, line):
dax.set_array(clipped_wofs.frequency[num, :, :])
cax.set_array(
(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
/ np.max(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
)
ax2.set_title("Time = " + str(clipped_wofs.frequency.coords["time"].values[(num)])[:12])
line.set_data(x[:num], y[:num])
return (line,)
x = tsw.index
y = tsw.values
(line,) = ax1.plot(x, y)
plt.tight_layout()
ani = animation.FuncAnimation(
fig, update, len(ds.time), fargs=[x, y, line], interval=800, blit=True
)
# ani.save('LakeChad_Nigeria_ts.gif')
plt.close()
HTML(ani.to_html5_video())
Conclusion
Ce notebook a présenté deux animations permettant de visualiser les masses d’eau sur plusieurs années afin de montrer les tendances à moyen et long terme. Les techniques et le code présentés ici peuvent être adaptés aux objectifs définis par l’utilisateur. Ces visualisations peuvent être utiles pour une analyse exploratoire des tendances de l’étendue des eaux sur plusieurs années. Une analyse plus formelle des tendances et des modèles devrait être entreprise pour mieux comprendre ces phénomènes, en particulier pour comprendre les effets de facteurs tels que les précipitations et les changements de couverture terrestre.
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 :
[ ]:
import datacube
print(datacube.__version__)
Dernier test :
[ ]:
from datetime import datetime
datetime.today().strftime("%Y-%m-%d")