Animations en vraies couleurs avec des précipitations

Mots-clés : données utilisées ; sentinel-2, données utilisées ; GeoMAD roulant, ensembles de données ; CHIRPS, climat, précipitations, données utilisées ; géomédiane sentinel-2, données utilisées ; géomédiane landsat 8 et 9, données utilisées ; géomédiane landsat 8, données utilisées ; géomédiane landsat 5 et 7

Aperçu

Il peut être instructif de visualiser des images en vraies couleurs d’une zone au fil du temps, parallèlement à une variable qui contribue à l’état de la surface terrestre, comme les précipitations. L’inspection d’une série chronologique animée d’une zone donnée au fil du temps, parallèlement aux précipitations cumulées, peut montrer comment les caractéristiques du paysage (végétation, plans d’eau, etc.) réagissent aux épisodes de pluie ou à leur absence.

Ce cahier montre comment les données de précipitations et d’imagerie en couleurs réelles peuvent être assimilées pour produire des animations informatives.

Description

Nous pouvons utiliser la réflectivité de surface et les données de précipitations CHIRPS pour produire ces animations. Grâce au DE Africa Sandbox, nous pouvons intégrer ces ensembles de données sur une période donnée et les aligner temporellement et spatialement pour les visualiser de manière synchrone.

Les animations sont réalisées dans ce notebook en suivant les étapes ci-dessous :

  1. Sélectionnez une zone d’intérêt.

  2. Charger des images de réflectance de surface (vraies couleurs).

  3. Rééchantillonnez et interpolez les images pour générer une série chronologique régulière.

  4. Charger les données de précipitations.

  5. Rééchantillonner les données de précipitations pour générer une série chronologique régulière.

  6. Définir les axes et la fonction de mise à jour pour l’animation.

  7. Exécutez et interprétez l’animation.

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

[1]:

%matplotlib inline

import datacube
import datetime as dt
import numpy as np
import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import matplotlib.dates as mdates

from IPython.display import HTML
from odc.geo.geom import Geometry
from deafrica_tools.plotting import display_map, rgb
from deafrica_tools.areaofinterest import define_area
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.datahandling import load_ard

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='animation-rainfall')

[3]:

create_local_dask_cluster()

/opt/venv/lib/python3.12/site-packages/distributed/node.py:187: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 40959 instead
  warnings.warn(

Client

Client-73720b06-d4a1-11ef-8a95-1ede66a139be

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /user/victoria@kartoza.com/proxy/40959/status

Cluster Info

Paramètres d’analyse

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

  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 39c304c1e28d438d821bf2fa35d74214 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.

La zone par défaut est Mosiotunya - la fumée qui tonne - (chutes Victoria) sur le fleuve Zambèze à la frontière du Zimbabwe et de la Zambie.

[4]:

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-17.92, lon=25.85, 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)

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

Charger des données satellite à partir de Datacube

Dans un premier temps, nous chargerons les bandes géomédianes en couleurs réelles (rouge, vert, bleu) pour visualiser les changements dans le paysage au cours d’une année récente. Nous chargeons un mois avant et après la période d’intérêt, car cela permet de générer une série chronologique régulière au cours de l’année d’intérêt.

[5]:

# Create a reusable query
query = {
    'x': lon_range,
    'y': lat_range,
    'time': ('2022-01-01', '2022-12-31'),
    'resolution': (-20, 20)
}

# Load available data
ds = dc.load(product=['gm_s2_rolling'],
             measurements=['red', 'green', 'blue'],
             group_by='solar_day',
             output_crs='epsg:6933',
             **query)

# Print output data
ds

[5]:

<xarray.Dataset> Size: 25MB
Dimensions:      (time: 14, y: 609, x: 483)
Coordinates:
  * time         (time) datetime64[ns] 112B 2021-12-16T23:59:59.999999 ... 20...
  * y            (y) float64 5kB -2.244e+06 -2.244e+06 ... -2.256e+06 -2.256e+06
  * x            (x) float64 4kB 2.489e+06 2.489e+06 ... 2.499e+06 2.499e+06
    spatial_ref  int32 4B 6933
Data variables:
    red          (time, y, x) uint16 8MB 1257 1299 1223 1177 ... 411 370 357 368
    green        (time, y, x) uint16 8MB 1015 1021 964 911 ... 665 638 609 611
    blue         (time, y, x) uint16 8MB 727 755 684 660 736 ... 384 364 359 362
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

Interpolation de séries chronologiques

Nous souhaitons générer une animation pour une seule année civile. Nous avons donc besoin d’une série chronologique régulière avec des intervalles de temps correspondant aux données de précipitations. Cette fonction n’est généralement pas disponible à partir d’un ensemble de données prêt à être analysé, car certaines images seront omises en raison de la couverture nuageuse, en particulier pendant les saisons les plus humides.

Pour réguler les séries temporelles, nous sous-échantillonnons la résolution temporelle à l’aide d’une interpolation. Cela signifie que nous générons de nouvelles données, ou que nous comblons les lacunes au fil du temps, sur la base des données existantes. Il est important que cela soit communiqué aux téléspectateurs et aux utilisateurs, car les données interpolées ne sont pas observationnelles et sont donc sujettes à incertitude.

Il existe plusieurs méthodes d’interpolation : « linéaire », « la plus proche », « zéro », « s-linéaire », « quadratique », « cubique » et « polynomiale ». La méthode détermine la flexibilité ou « ondulation » autorisée dans le remplissage temporel des espaces vides. Vous pouvez essayer quelques méthodes différentes ci-dessous. Vous trouverez de plus amples informations dans la documentation scipy.

Dans ce cas, nous générons une série temporelle avec une image pour chaque deux jours de l’année, ce qui donne un ensemble de données avec 182 pas de temps (environ = 365/2). Cela donne une animation fluide et agréable sans imposer une demande excessive en termes de calcul et de mémoire.

[6]:

ds = ds.resample(time="2D").interpolate("linear").sel(
    time=slice('2022-01-01', '2022-12-31'))

ds

[6]:

<xarray.Dataset> Size: 1GB
Dimensions:      (y: 609, x: 483, time: 183)
Coordinates:
  * y            (y) float64 5kB -2.244e+06 -2.244e+06 ... -2.256e+06 -2.256e+06
  * x            (x) float64 4kB 2.489e+06 2.489e+06 ... 2.499e+06 2.499e+06
    spatial_ref  int32 4B 6933
  * time         (time) datetime64[ns] 1kB 2022-01-01 2022-01-03 ... 2022-12-31
Data variables:
    red          (time, y, x) float64 431MB 1.082e+03 1.139e+03 ... 443.4 443.0
    green        (time, y, x) float64 431MB 967.4 978.1 937.1 ... 682.4 676.7
    blue         (time, y, x) float64 431MB 667.5 705.3 644.7 ... 433.0 429.8
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

Charger les données de précipitations

Nous allons charger les données de précipitations quotidiennes à partir de CHIRPS, puis les rééchantillonner pour additionner les précipitations tous les deux jours, de sorte que la résolution temporelle corresponde à celle de l’imagerie en vraies couleurs. Nous allons également accumuler les précipitations tout au long de l’année à l’aide de la fonction « cumsum() ». Cela nous permet de visualiser les précipitations totales, ce qui peut être plus informatif que la visualisation d’événements pluvieux sporadiques. Cependant, il existe de nombreuses autres façons de visualiser les précipitations au fil du temps que vous souhaiterez peut-être explorer.

Un graphique statique des précipitations annuelles cumulées est généré ci-dessous. Pour l’emplacement par défaut de Mosiotunya (chutes Victoria), une saison sèche distincte au milieu de l’année est évidente.

[7]:

ds_rf_daily = dc.load(product='rainfall_chirps_daily',
                time='2022',
                y = lat_range,
                x = lon_range,
                resolution=(-5000, 5000),
                output_crs='epsg:6933')

ds_rf_daily['rainfall'].resample(time='2D').sum().mean(['y', 'x']).cumsum().plot(figsize=(16, 4))
plt.xlabel('Day')
plt.ylabel('%s (%s)' % ('Total Precipitation', ds_rf_daily['rainfall'].attrs['units']));
../../../_images/sandbox_notebooks_Real_world_examples_Animation_with_rainfall_18_0.png

Définir et exécuter l’animation

La cellule ci-dessous définit et exécute l’animation à travers les étapes suivantes :

  1. Création de deux sous-lots.

  2. Définition des paramètres des axes.

  3. Définition d’une fonction de mise à jour.

  4. Exécution et traçage de l’animation.

Vous pouvez éventuellement enregistrer l’animation dans un « .gif » en supprimant le # avant la commande « ani.save() ». Notez que cela ferme le tracé, vous devrez donc rajouter le # et exécuter à nouveau la cellule pour visualiser l’animation sous forme dynamique.

Pour la région par défaut, nous pouvons voir que la végétation est verte et que la cascade coule dans les premiers mois. Au début de la période sèche, la végétation vieillit, bien que le débit de la cascade mette plus de temps à ralentir.

[8]:

# create a figure and axes
fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(122)
ax2 = plt.subplot(121)

ax1.set_title("Cumulative Rainfall")
ax1.set_xlabel("Date")
ax1.set_ylabel("Total Precipitation (mm)")

bands=['red', 'green', 'blue']


cax = (ds[bands].isel(time=0).to_array().transpose('y','x','variable').squeeze()/np.max(
    ds[bands].isel(time=0).to_array(
    ).transpose('y','x','variable').squeeze())).plot.imshow(rgb='variable', animated = True, robust=True)

def update(num, x, y, line):
    cax.set_array((ds[bands].isel(time=num).to_array().transpose('y','x','variable').squeeze().clip(0,2500)/np.max(
        ds[bands].isel(time=num).to_array().transpose('y','x','variable').squeeze().clip(0,2500))))
    ax2.set_title("Time = " + str(ds[bands].coords['time'].values[(num)])[:12])
    line.set_data(x[:num], y[:num])
    return line,

x = ds_rf_daily['rainfall'].resample(time='2D').sum().mean(['y', 'x']).cumsum().time
y = ds_rf_daily['rainfall'].resample(time='2D').sum().mean(['y', 'x']).cumsum()

line, = ax1.plot(x, y)

plt.tight_layout()

ani = animation.FuncAnimation(fig, update, len(x), fargs=[x, y, line],
                              interval=200, blit=True)
#ani.save('rgb_rain_anim.gif')

plt.close()
HTML(ani.to_html5_video())

[8]:

Précipitations annuelles à long terme

Nous pouvons également inspecter les modèles de précipitations à plus long terme parallèlement aux images géomédianes annuelles pour obtenir des informations sur les modèles interannuels des précipitations et des paysages.

Charger les géomédianes

Les produits géomédians dérivés de Landsat nous permettent d’inspecter de longues périodes de temps. Dans ce cas, nous chargeons des images annuelles de Landsat 5, 7, 8 et 9 couvrant la période 1990-2021.

[9]:

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,
             dask_chunks={'time': 1, 'x': 3000, 'y': 3000},
             time=("1990", "2021")
             ).compute()
ds

[9]:

<xarray.Dataset> Size: 25MB
Dimensions:      (time: 32, y: 406, x: 322)
Coordinates:
  * time         (time) datetime64[ns] 256B 1990-07-02T11:59:59.999999 ... 20...
  * y            (y) float64 3kB -2.244e+06 -2.244e+06 ... -2.256e+06 -2.256e+06
  * x            (x) float64 3kB 2.489e+06 2.489e+06 ... 2.499e+06 2.499e+06
    spatial_ref  int32 4B 6933
Data variables:
    red          (time, y, x) uint16 8MB 1008 994 997 1127 ... 1044 1022 974
    green        (time, y, x) uint16 8MB 814 798 791 864 988 ... 790 789 769 750
    blue         (time, y, x) uint16 8MB 562 551 539 580 645 ... 505 502 485 472
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

Charger les précipitations mensuelles

Dans cette animation, nos pas de temps sont annuels, nous allons donc générer une série chronologique des précipitations totales annuelles. Pour ce faire, nous chargeons les précipitations annuelles et additionnons les valeurs par année.

Nous calculerons également l’anomalie des précipitations, car cela peut être un moyen utile de visualiser la relativité des précipitations annuelles. Pour cette raison, nous calculons et stockons la moyenne et l’écart type des précipitations annuelles sur toute la période.

[10]:

ds_rf_year = dc.load(product='rainfall_chirps_monthly',
                time=("1990", "2021"),
                y = lat_range,
                x = lon_range,
                resolution=(-5000, 5000),
                dask_chunks={'time': 1, 'x': 3000, 'y': 3000},
                output_crs='EPSG:6933').resample(time='1YE').sum()

# mean across years
rf_mean = ds_rf_year.mean('time').compute()

#calculate annual std dev
rf_std = ds_rf_year.std('time').compute()

Calculer les anomalies standardisées

La fonction ci-dessous calcule les anomalies standardisées en utilisant les valeurs moyennes et d’écart type stockées ci-dessus.

[11]:

stand_anomalies = xr.apply_ufunc(
    lambda x, m, s: (x - m) / s,
    ds_rf_year,
    rf_mean,
    rf_std,
    output_dtypes=[rf_mean.rainfall.dtype],
    dask="allowed"
).compute().mean(['y', 'x']).to_dataframe().drop(
    ['spatial_ref'], axis=1).fillna(0)

Préparer les données pluviométriques

Enfin, nous réduisons les dimensions spatiales des données de précipitations pour générer un cadre de données contenant des colonnes pour le temps et le volume annuel de précipitations (mm).

[12]:

ds_rf_year_1D = ds_rf_year.mean(['y', 'x']).to_dataframe().drop(
    ['spatial_ref'], axis=1).fillna(0)

Définir et exécuter l’animation

La cellule ci-dessous définit et exécute l’animation à travers les étapes suivantes :

  1. Création de deux sous-lots.

  2. Définition des paramètres des axes.

  3. Définition d’une fonction de mise à jour.

  4. Exécution et traçage de l’animation.

Vous pouvez éventuellement enregistrer l’animation dans un « .gif » en supprimant le # avant la commande « ani.save() ». Notez que cela ferme le tracé, vous devrez donc rajouter le # et exécuter à nouveau la cellule pour visualiser l’animation sous forme dynamique.

[13]:

fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(10)

ax1 = plt.subplot2grid(shape=(2, 2), loc=(0, 1), colspan=2)
ax2 = plt.subplot2grid(shape=(2, 2), loc=(0, 0), rowspan=2)
ax3 = plt.subplot2grid(shape=(2, 2), loc=(1, 1), colspan=2)

bands=['red', 'green', 'blue']

ax1.xaxis_date()
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
ax1.set_xlim((min(mdates.date2num(ds_rf_year_1D.index))), (max(mdates.date2num(ds_rf_year_1D.index))))
ax1.set_ylim(0, max(ds_rf_year_1D.rainfall)+100)
ax1.set_title("Annual Rainfall")
ax1.set_xlabel("Year")
ax1.set_ylabel("Rainfall (mm)")

ax3.set_ylim([-4, 4])
ax3.xaxis_date()
ax3.set_xlim((min(mdates.date2num(ds_rf_year_1D.index))), (max(mdates.date2num(ds_rf_year_1D.index))))
ax3.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
ax3.set_title("Rainfall Anomaly")
ax3.set_xlabel("Year")
ax3.set_ylabel("Std Deviations from annual mean")
ax3.axhline(0, color='black', linestyle='--')

x = ds_rf_year_1D.index.strftime("%Y")
y = [0] * ds_rf_year_1D.rainfall

x2 = stand_anomalies.index.strftime("%Y")
y2 = [0] * stand_anomalies.rainfall
bars = ax3.bar(x2, y2, align='center', width=80,color=(stand_anomalies.rainfall > 0).map({True: 'g', False: 'brown'}))

bars2 = ax1.bar(x, y, align='center', width=80, color='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)

def update(num, x, y, x2, y2, bars, bars2):
    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))
    y2.iloc[num] = stand_anomalies.rainfall.iloc[num]
    bars[num].set_height(y2.iloc[num])
    ax2.set_title("Year = " + str(ds[bands].coords['time'].values[(num)])[:4])
    y.iloc[num] = ds_rf_year_1D.rainfall.iloc[num]
    bars2[num].set_height(y.iloc[num])
    return bars2

plt.tight_layout()
ani = animation.FuncAnimation(fig, update, len(x), fargs=[x, y, x2, y2, bars, bars2],
                              interval=400, blit=False, repeat=True)
#ani.save('geom_rain_lt_anim.gif')

plt.close()

HTML(ani.to_html5_video())

[13]:

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

[14]:

print(datacube.__version__)

1.8.20

Dernier test :

[15]:

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

[15]:

'2025-01-17'