Cartographie de l’étendue des eaux et des précipitations à l’aide de Sentinel-1 et CHIRPS#

Produits utilisés : s1_rtc, rainfall_chirps_monthly

**Keywords**: :index:`data used; sentinel-1`, :index:`data used; CHIRPS`, :index:`water; extent`, :index:`analysis; time series`

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. The spatial extent of water-related ecosystems
ii. The quantity of water contained within these ecosystems
iii. The quality of water within these ecosystems
iv. The health or state of these ecosystems

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

Description#

Le bloc-notes charge les données Sentinel-1 et, à l’aide de l’indice Sentinel-Water, cartographie l’étendue spatiale des plans d’eau. Il charge et trace également les précipitations totales mensuelles à partir de CHIRPS. La dernière section comparera l’étendue de l’eau entre deux périodes pour permettre de visualiser où les changements se produisent.

Charger des paquets#

Importez les packages Python utilisés pour l’analyse.

[1]:

%matplotlib inline

import datacube
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

from deafrica_tools.datahandling import load_ard
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.datahandling import wofs_fuser

from long_term_water_extent import (
    load_vector_file,
    get_resampled_labels,
    resample_water_observations,
    resample_rainfall_observations,
    calculate_change_in_extent,
    compare_extent_and_rainfall,
)

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.

[ ]:

create_local_dask_cluster()

Se connecter à Data Cube#

[3]:

dc = datacube.Datacube(app="long_term_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.

Téléchargez un fichier vectoriel pour votre étendue d’eau et votre bassin versant dans le dossier « données ».
Définissez la plage horaire que vous souhaitez utiliser.
Définissez la stratégie de rééchantillonnage. Les options possibles sont les suivantes :
- « 1Y » - Rééchantillonnage annuel, utilisez cette option pour une surveillance à plus long terme
- « QS-DEC » - Rééchantillonnage trimestriel à partir de décembre
- « 3M » - Rééchantillonnage trimestriel
- « 1M » - Rééchantillonnage mensuel

Pour plus de détails sur le rééchantillonnage des périodes de temps, consultez la documentation xarray et pandas.

[4]:

water_extent_vector_file = "data/lake_baringo_extent.geojson"

water_catchment_vector_file = "data/lake_baringo_catchment.geojson"

time_range = ("2018-07", "2021") #earliest date of S1 is July 1st 2018

resample_strategy = "Q-DEC"

dask_chunks = dict(x=1000, y=1000)

Obtenir les géométries des plans d’eau et des bassins versants#

La cellule suivante extraira les géométries du plan d’eau et du bassin versant à partir des fichiers vectoriels fournis, qui seront utilisés pour charger les observations de l’eau depuis l’espace et les produits de précipitations CHIRPS.

[5]:

extent, extent_geometry = load_vector_file(water_extent_vector_file)
catchment, catchment_geometry = load_vector_file(water_catchment_vector_file)

Charger Sentinel-1 pour le plan d’eau#

La première étape consiste à charger les données Sentinel-1 à l’aide de la géométrie d’étendue.

[6]:

extent_query = {
    "time": time_range,
    "resolution": (-20, 20),
    "output_crs": "EPSG:6933",
    "geopolygon": extent_geometry,
    "group_by": "solar_day",
    "dask_chunks":dask_chunks
}

ds = load_ard(dc=dc, products=["s1_rtc"], measurements=["vv", "vh"], **extent_query)

Using pixel quality parameters for Sentinel 1
Finding datasets
    s1_rtc
Applying pixel quality/cloud mask
Returning 197 time steps as a dask array

Convertir les nombres numériques Sentinel-1 en dB#

Bien que la rétrodiffusion Sentinel-1 soit fournie sous forme d’intensité linéaire, il est souvent utile de convertir la rétrodiffusion en décibels (dB) pour l’analyse. La rétrodiffusion en dB présente un profil de bruit plus symétrique et une distribution de valeurs moins biaisée pour une évaluation statistique plus facile.

Les données de rétrodiffusion Sentinel-1 sont converties du nombre numérique (DN) en rétrodiffusion en décibels (dB) à l’aide de la fonction :

\begin{equation} 10 * \log_{10}(\text{DN}) \end{equation}

[7]:

# Convert DN to db values.
ds["vv_db"] = 10 * np.log10(ds.vv)
ds["vh_db"] = 10 * np.log10(ds.vh)

Calculer l’indice d’eau SWI#

L’indice d’eau Sentinel-1A (SWI) est calculé comme suit :

\begin{equation} \text{SWI} = 0,1747 * \beta _{vv} + 0,0082 * \beta _{vh} * \beta _{vv} + 0,0023 * \beta _{vv}^{2} - 0,0015 * \beta _{vh}^{2} + 0,1904 \end{equation}

où βvh et βvv représentent respectivement le coefficient de rétrodiffusion en polarisation VH et en polarisation VV (Tian et al., 2017).

[8]:

# Calculate the Sentinel-1A Water Index (SWI).
ds['swi'] = (
    (0.1747 * ds.vv_db)
    + (0.0082 * ds.vh_db * ds.vv_db)
    + (0.0023 * ds.vv_db ** 2)
    - (0.0015 * ds.vh_db ** 2)
    + 0.1904
)

swi = ds[['swi']]
print(swi)

<xarray.Dataset>
Dimensions:      (time: 197, y: 1818, x: 852)
Coordinates:
  * time         (time) datetime64[ns] 2018-07-01T15:56:48.569171 ... 2021-12...
  * y            (y) float64 9.433e+04 9.431e+04 ... 5.801e+04 5.799e+04
  * x            (x) float64 3.473e+06 3.473e+06 3.473e+06 ... 3.49e+06 3.49e+06
    spatial_ref  int32 6933
Data variables:
    swi          (time, y, x) float32 dask.array<chunksize=(1, 1000, 852), meta=np.ndarray>
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Identifier l’eau à chaque période de rééchantillonnage#

La deuxième étape consiste à rééchantillonner les observations pour obtenir une mesure cohérente de la masse d’eau, puis à calculer la quantité d’eau classée pour chaque période.

[9]:

resampled_water_ds, resampled_water_area_ds = resample_water_observations(
    swi, resample_strategy, radar=True
)

date_range_labels = get_resampled_labels(ds, resample_strategy)

/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in add
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in subtract
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in add
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in subtract
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in add
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in subtract
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in log10
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in subtract
  return func(*(_execute_task(a, cache) for a in args))
/usr/local/lib/python3.10/dist-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(

Automatic SWI threshold: 0.59

Tracer l’évolution de la surface de l’eau au fil du temps#

[10]:

fig, ax = plt.subplots(figsize=(15, 5))

ax.plot(
    date_range_labels,
    resampled_water_area_ds.values,
    color="red",
    marker="^",
    markersize=4,
    linewidth=1,
)
plt.xticks(date_range_labels, rotation=65)
plt.title(f"Observed Area of Water from {time_range[0]} to {time_range[1]}")
plt.ylabel("Waterbody area (km$^2$)")
plt.tight_layout()

../../../../_images/sandbox_notebooks_Use_cases_Monitoring_water_extent_Monitoring_water_extent_radar_25_0.png

Charger les précipitations mensuelles CHIRPS#

[11]:

catchment_query = {
    "time": time_range,
    "resolution": (-5000, 5000),
    "output_crs": "EPSG:6933",
    "geopolygon": catchment_geometry,
    "group_by": "solar_day",
    "dask_chunks":dask_chunks
}

rainfall_ds = dc.load(product="rainfall_chirps_monthly", **catchment_query)

Rééchantillonner pour estimer les précipitations pour chaque période#

Cela se fait en calculant la pluviométrie moyenne sur l’étendue du bassin versant, puis en additionnant ces moyennes sur la période de rééchantillonnage pour estimer la pluviométrie totale pour le bassin versant.

[12]:

catchment_rainfall_resampled_ds = resample_rainfall_observations(
    rainfall_ds, resample_strategy, catchment
)

Comparer la superficie du plan d’eau aux précipitations du bassin versant#

Cette étape trace la somme des précipitations moyennes pour le bassin versant sur chaque période sous forme d’histogramme, superposée à la superficie du plan d’eau calculée précédemment.

[13]:

figure = compare_extent_and_rainfall(
    resampled_water_area_ds, catchment_rainfall_resampled_ds, "mm", date_range_labels
)

../../../../_images/sandbox_notebooks_Use_cases_Monitoring_water_extent_Monitoring_water_extent_radar_31_0.png

Sauver la figure#

[14]:

figure.savefig("waterarea_and_rainfall_radar.png", bbox_inches="tight")

Comparer l’étendue d’eau pour deux périodes différentes#

Pour l’étape suivante, entrez une date de référence et une date d’analyse pour construire un graphique montrant où l’eau est apparue et a disparu, en comparant les deux dates.

[15]:

baseline_time = "2018-07-01"
analysis_time = "2021-10-01"

[16]:

figure = calculate_change_in_extent(baseline_time, analysis_time, resampled_water_ds, radar=True)

../../../../_images/sandbox_notebooks_Use_cases_Monitoring_water_extent_Monitoring_water_extent_radar_36_0.png

Enregistrer la figure#

[17]:

figure.savefig("waterarea_change_radar.png", bbox_inches="tight")

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 :

[18]:

print(datacube.__version__)

1.8.15

Dernier test :

[19]:

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

[19]:

'2023-08-21'