Modélisation de l’élévation intertidale à l’aide de données de marée

  • Produits utilisés : s2_l2a

  • Prérequis : Pour plus d’informations sur l’étape de modélisation des marées de cette analyse, reportez-vous au « Cahier de modélisation des marées <../Frequently_used_code/Tidal_modelling.ipynb> » __

Mots-clés : données utilisées ; Sentinel-2, eau ; modélisation des marées, eau ; extraction de la ligne de flottaison, indice de bande ; NDWI, méthodes de données ; composition

Aperçu

Les zones intertidales abritent des habitats écologiques importants (p. ex. plages et rivages sablonneux, vasières, rivages rocheux et récifs) et offrent de nombreux avantages précieux tels que la protection contre les ondes de tempête, le stockage du carbone et les ressources naturelles à usage récréatif et commercial. Cependant, les zones intertidales sont confrontées à des menaces croissantes liées à l’érosion côtière, à la remise en valeur des terres (p. ex. construction de ports) et à l’élévation du niveau de la mer. Des données d’élévation précises décrivant la hauteur et la forme du littoral sont nécessaires pour aider à prédire quand et où ces menaces auront le plus d’impact. Cependant, ces données sont coûteuses et difficiles à cartographier sur l’ensemble de la zone intertidale d’un continent de la taille de l’Afrique.

Cas d’utilisation de Digital Earth Africa

La montée et la descente de la marée peuvent être utilisées pour révéler la forme tridimensionnelle du littoral en cartographiant la frontière entre l’eau et la terre sur une gamme de marées connues (par exemple de la marée basse à la marée haute). En supposant que la frontière terre-eau est une ligne de hauteur constante par rapport au niveau moyen de la mer (MSL), les élévations peuvent être modélisées pour la zone du littoral située entre la marée la plus basse et la marée la plus haute observée.

Les images provenant de satellites tels que le programme Copernicus Sentinel-2 sont disponibles gratuitement sur l’ensemble de la planète, ce qui fait de l’imagerie satellite un outil puissant et rentable pour modéliser la forme et la structure 3D de la zone intertidale à l’échelle régionale ou nationale.

Référence : Bishop-Taylor et al. 2019, le premier modèle 3D de l’ensemble du littoral australien : le National Intertidal Digital Elevation Model.

Description

Dans cet exemple, nous démontrons une version simplifiée de la méthode NIDEM (National Intertidal Digital Elevation Model) qui combine les données du satellite Sentinel-2 avec la modélisation des marées, la composition d’images et les techniques d’interpolation spatiale. Nous cartographions d’abord la frontière entre la terre et l’eau de la marée basse à la marée haute, et utilisons ces informations pour générer des cartes d’élévation 3D continues et fluides de la zone intertidale. Les données obtenues peuvent aider à cartographier les habitats des espèces côtières menacées, à identifier les zones d’érosion côtière, à planifier les événements extrêmes tels que les ondes de tempête et les inondations, et à améliorer les modèles de la façon dont l’élévation du niveau de la mer affectera le littoral. Cet exemple pratique guide les utilisateurs à travers le code requis pour :

  1. Chargez dans une série chronologique Sentinel-2 sans nuage.

  2. Calculer un indice hydrique (NDWI).

  3. Étiquetez et triez les images satellites par hauteur de marée.

  4. Créez des images « récapitulatives » ou composites qui montrent la répartition des terres et de l’eau à des intervalles discrets de l’amplitude des marées (par exemple, à marée basse, à marée haute).

  5. Extraire et visualiser la topographie de la zone intertidale sous forme de courbes de profondeur.

  6. Interpoler les contours de profondeur dans un modèle numérique d’élévation (DEM) lisse et continu de la zone intertidale.

Commencer

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Une fois l’analyse terminée, revenez à la cellule « Paramètres d’analyse », modifiez certaines valeurs (par exemple, choisissez un autre lieu ou une autre période à analyser) et relancez l’analyse. Vous trouverez des instructions supplémentaires sur la modification du bloc-notes à la fin.

Charger des paquets

Chargez les principaux packages Python et les fonctions de support pour l’analyse.

[1]:

%matplotlib inline

import os
import datacube
import xarray as xr
import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from odc.geo.cog import write_cog
from odc.geo.geom import Geometry

from deafrica_tools.datahandling import load_ard, mostcommon_crs
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.coastal import tidal_tag
from deafrica_tools.spatial import subpixel_contours, interpolate_2d, contours_to_arrays
from deafrica_tools.plotting import display_map, map_shapefile, rgb
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.areaofinterest import define_area

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 utiliser Dask, configurez le cluster de calcul local à l’aide de la cellule ci-dessous.

[2]:

create_local_dask_cluster()

Client

Client-546f4068-d3dd-11ef-9a93-1e01ca941b00

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

Cluster Info

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.

[3]:

dc = datacube.Datacube(app='Intertidal_elevation_Sentinel-2')

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 l’analyse. Les paramètres sont

  • « lat » : la latitude centrale à analyser (par exemple « -11,384 »).

  • « lon » : la longitude centrale à analyser (par exemple « 43,287 »).

  • « buffer » : le nombre de degrés carrés à charger autour de la latitude et de la longitude centrales. Pour des temps de chargement raisonnables, définissez cette valeur sur « 0,1 » ou moins.

  • time_range : la plage de dates à analyser (par exemple ('2017', '2022'))

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 a5f2ea26fbec4604a71ed52f83c54b29 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 génère un modèle d’élévation intertidale pour une zone au large des côtes de la Guinée-Bissau.

Pour garantir que la partie modélisation des marées de cette analyse fonctionne correctement, assurez-vous que le centre de la zone d’étude est situé au-dessus de l’eau lors du réglage de « lat » et « lon ».

[4]:

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-19.25612, lon=44.32470, buffer=0.04)

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

# Set the range of dates for the analysis
time_range = ('2017', '2022')

Afficher l’emplacement sélectionné

La cellule suivante affichera la zone sélectionnée sur une carte interactive. N’hésitez pas à zoomer et dézoomer pour mieux comprendre la zone que vous allez analyser. En cliquant sur n’importe quel point de la carte, vous découvrirez les coordonnées de latitude et de longitude de ce point.

[5]:

display_map(x=lon_range, y=lat_range)

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

Charger les données Sentinel-2 masquées par le cloud

La première étape de cette analyse consiste à charger les données Sentinel-2 pour les « lat_range », « lon_range » et « time_range » que nous avons fournies ci-dessus. Le code ci-dessous utilise la fonction « load_ard » pour charger les données prêtes pour l’analyse Sentinel-2 pour la zone et l’heure spécifiées. Pour plus d’informations, consultez le bloc-notes « Using load_ard <../Frequently_used_code/Using_load_ard.ipynb> ». La fonction masquera également automatiquement les nuages de l’ensemble de données, ce qui nous permettra de nous concentrer sur les pixels qui contiennent des données utiles :

[6]:

# Create the 'query' dictionary object, which contains the longitudes,
# latitudes and time provided above
query = {
    'x': lon_range,
    'y': lat_range,
    'time': time_range,
    'measurements': ['red', 'green', 'blue', 'nir'],
    'resolution': (-10, 10),
}

# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='s2_l2a', query=query)

# Load available data
s2_ds = load_ard(dc=dc,
                      products=['s2_l2a'],
                      output_crs=output_crs,
                      min_gooddata=0.5,
                      group_by='solar_day',
                      dask_chunks={},
                      **query)

/home/jovyan/dev/deafrica-sandbox-notebooks/Tools/deafrica_tools/datahandling.py:244: UserWarning: Setting 'min_gooddata' percentage to > 0.0 will cause dask arrays to compute when loading pixel-quality data to calculate 'good pixel' percentage. This can slow the return of your dataset.
  warnings.warn(

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Counting good quality pixels for each time step
Filtering to 642 out of 821 time steps with at least 50.0% good quality pixels
Applying pixel quality/cloud mask
Returning 642 time steps as a dask array

Une fois le chargement terminé, examinez les données en les imprimant dans la cellule suivante. L’argument « Dimensions » révèle le nombre d’intervalles de temps dans l’ensemble de données, ainsi que le nombre de pixels dans les dimensions « x » (longitude) et « y » (latitude).

[7]:

print(s2_ds)

<xarray.Dataset> Size: 8GB
Dimensions:      (time: 642, y: 889, x: 845)
Coordinates:
  * time         (time) datetime64[ns] 5kB 2017-01-07T07:23:27 ... 2022-12-27...
  * y            (y) float64 7kB 7.875e+06 7.875e+06 ... 7.866e+06 7.866e+06
  * x            (x) float64 7kB 4.248e+05 4.248e+05 ... 4.332e+05 4.333e+05
    spatial_ref  int32 4B 32738
Data variables:
    red          (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    green        (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    blue         (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    nir          (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
Attributes:
    crs:           epsg:32738
    grid_mapping:  spatial_ref

Exemple de tracé de pas de temps en vraies couleurs

Pour visualiser les données, utilisez la fonction utilitaire « rgb » préchargée pour tracer une image en vraies couleurs pour un intervalle de temps donné. Les zones blanches indiquent les endroits où les nuages ou autres pixels non valides de l’image ont été masqués.

Modifiez la valeur de « timestep » et réexécutez la cellule pour tracer un pas de temps différent (les pas de temps sont numérotés de « 0 » à « n_time - 1 » où « n_time » est le nombre total de pas de temps ; voir la liste « time » sous la catégorie « Dimensions » dans l’impression du jeu de données ci-dessus).

[8]:

# Set the timesteps to visualise
timestep = 12

# Generate RGB plots at each timestep
rgb(s2_ds, index=timestep)
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_19_0.png

Calculer l’indice d’eau par différence normalisée

Pour extraire les contours de profondeur intertidale, nous devons être capables de séparer l’eau de la terre dans notre zone d’étude. Pour ce faire, nous pouvons utiliser nos données Sentinel-2 pour calculer un indice d’eau appelé « Normalised Difference Water Index », ou NDWI. Cet indice utilise le rapport entre le rayonnement vert et le rayonnement proche infrarouge pour identifier la présence d’eau. La formule est :

\[\begin{aligned} \text{NDWI} &= \frac{(\text{Vert} - \text{NIR})}{(\text{Vert} + \text{NIR})} \end{aligned}\]

où « Vert » est la bande verte et « NIR » est la bande proche infrarouge.

Pour interpréter l’indice, les valeurs élevées (supérieures à 0, couleurs bleues) représentent généralement des pixels d’eau, tandis que les valeurs faibles (inférieures à 0, couleurs rouges) représentent la terre. Vous pouvez utiliser la cellule ci-dessous pour calculer et tracer l’une des images après avoir calculé l’indice.

[9]:

# Calculate the water index
s2_ds = calculate_indices(s2_ds, index='NDWI',
                               satellite_mission ='s2')

# Plot the resulting image for the same timestep selected above
s2_ds.NDWI.isel(time=timestep).plot(cmap='RdBu',
                                          size=6,
                                          vmin=-0.8,
                                          vmax=0.8)
plt.show()
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_21_0.png

Comment le tracé de l’indice se compare-t-il à l’image optique précédente ? Y avait-il de l’eau ou de la terre à un endroit auquel vous ne vous attendiez pas ?

Modèles de hauteurs de marée

L’emplacement du littoral peut varier considérablement entre la marée basse et la marée haute. Dans le code ci-dessous, nous cherchons à calculer la hauteur de la marée au moment exact où chaque image Sentinel-2 a été acquise. Cela nous permettra de construire une série chronologique triée d’images prises à marée basse et à marée haute, que nous utiliserons pour générer le modèle d’élévation intertidale.

La fonction « tidal_tag » ci-dessous utilise le modèle de marée OTPS TPXO8 <https://www.tpxo.net/global/tpxo8-atlas>`__ pour calculer la hauteur de la marée au moment exact où chaque image satellite de notre ensemble de données a été prise, et l’ajoute comme un nouvel attribut « tide_m » dans notre ensemble de données (pour plus d’informations sur cette fonction, reportez-vous au « cahier de modélisation des marées <../Frequently_used_code/Tidal_modelling.ipynb>`__).

Remarque : cette fonction ne peut modéliser correctement les marées que si le centre de votre zone d’étude est situé au-dessus de l’eau. Si ce n’est pas le cas, vous pouvez spécifier un emplacement de modélisation de marée personnalisé en transmettant une coordonnée à « tidepost_lat » et « tidepost_lon » (par exemple, « tidepost_lat=11.228, tidepost_lon=-15.860 »).

[10]:

# Calculate tides for each timestep in the satellite dataset
s2_ds = tidal_tag(ds=s2_ds, tidepost_lat=None, tidepost_lon=None)

# Print the output dataset with new `tide_m` variable
print(s2_ds)

Setting tide modelling location from dataset centroid: 44.32, -19.26
Modelling tides using FES2014 tidal model
<xarray.Dataset> Size: 10GB
Dimensions:      (time: 642, y: 889, x: 845)
Coordinates:
  * time         (time) datetime64[ns] 5kB 2017-01-07T07:23:27 ... 2022-12-27...
  * y            (y) float64 7kB 7.875e+06 7.875e+06 ... 7.866e+06 7.866e+06
  * x            (x) float64 7kB 4.248e+05 4.248e+05 ... 4.332e+05 4.333e+05
    spatial_ref  int32 4B 32738
Data variables:
    red          (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    green        (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    blue         (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    nir          (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    NDWI         (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    tide_m       (time) float64 5kB 0.6238 -0.9891 -0.7224 ... -1.41 0.1797
Attributes:
    crs:           epsg:32738
    grid_mapping:  spatial_ref

[11]:

# Plot the resulting tide heights for each Sentinel-2 image:
s2_ds.tide_m.plot()
plt.show()
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_25_0.png

Créer des images récapitulatives de l’indice hydrique de la marée basse à la marée haute

En utilisant ces hauteurs de marée, nous pouvons trier notre ensemble de données Sentinel-2 par hauteur de marée pour révéler quelles parties du paysage sont inondées ou exposées de la marée basse à la marée haute.

Les images individuelles de télédétection peuvent être affectées par le bruit, notamment les nuages, les reflets du soleil et les mauvaises conditions de qualité de l’eau (par exemple, les sédiments). Pour produire des images plus nettes qui peuvent être comparées plus facilement entre les stades de marée, nous pouvons créer des images « récapitulatives » ou composites qui combinent plusieurs images en une seule image pour révéler l’apparence « typique » ou médiane du paysage à différents stades de marée. Dans ce cas, nous utilisons la médiane comme statistique récapitulative, car elle empêche les valeurs aberrantes importantes (comme les nuages errants) de fausser les données, ce qui ne serait pas le cas si nous devions utiliser la moyenne.

Dans le code ci-dessous, nous prenons la série chronologique d’images, trions par marée et classons chaque image en 9 intervalles de marée discrets, allant des marées les plus basses (intervalle de marée 1) aux marées les plus hautes observées par Sentinel-2 (intervalle de marée 9). Pour plus d’informations sur cette méthode, reportez-vous à Sagar et al. 2018.

[12]:

# Sort every image by tide height
s2_ds = s2_ds.sortby('tide_m')

# Bin tide heights into 9 tidal intervals from low (1) to high tide (9)
binInterval = np.linspace(s2_ds.tide_m.min(),
                          s2_ds.tide_m.max(),
                          num=10)
tide_intervals = pd.cut(s2_ds.tide_m,
                        bins=binInterval,
                        labels=range(1, 10),
                        include_lowest=True)

# Add interval to dataset
s2_ds['tide_interval'] = xr.DataArray(tide_intervals.astype(int),
                                           [('time', s2_ds.time.data)])

print(s2_ds)

<xarray.Dataset> Size: 10GB
Dimensions:        (time: 642, y: 889, x: 845)
Coordinates:
  * time           (time) datetime64[ns] 5kB 2019-09-29T07:25:07 ... 2020-07-...
  * y              (y) float64 7kB 7.875e+06 7.875e+06 ... 7.866e+06 7.866e+06
  * x              (x) float64 7kB 4.248e+05 4.248e+05 ... 4.332e+05 4.333e+05
    spatial_ref    int32 4B 32738
Data variables:
    red            (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    green          (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    blue           (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    nir            (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    NDWI           (time, y, x) float32 2GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
    tide_m         (time) float64 5kB -1.995 -1.911 -1.901 ... 1.001 1.013
    tide_interval  (time) int64 5kB 1 1 1 1 1 1 1 1 1 1 ... 9 9 9 9 9 9 9 9 9 9
Attributes:
    crs:           epsg:32738
    grid_mapping:  spatial_ref

[13]:

s2_ds.sortby('time').tide_m.plot()
for i in binInterval: plt.axhline(i, c='black', alpha=0.5)
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_28_0.png

Maintenant que nous disposons d’un ensemble de données dans lequel chaque image est classée dans une plage distincte de la marée, nous pouvons combiner nos images en un ensemble de neuf images individuelles qui montrent où se trouvent la terre et l’eau de la marée basse à la marée haute. Cette étape peut prendre plusieurs minutes à traiter.

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.

[14]:

# For each interval, compute the median water index and tide height value
s2_intervals = (s2_ds[['tide_interval', 'NDWI', 'tide_m']]
                     .compute()
                     .groupby('tide_interval')
                     .median(dim='time'))

# Plot the resulting set of tidal intervals
s2_intervals.NDWI.plot(col='tide_interval', col_wrap=5, cmap='RdBu')
plt.show()

/opt/venv/lib/python3.12/site-packages/xarray/core/concat.py:540: UserWarning: No index created for dimension tide_interval because variable tide_interval is not a coordinate. To create an index for tide_interval, please first call `.set_coords('tide_interval')` on this object.
  ds.expand_dims(dim_name, create_index_for_new_dim=create_index_for_new_dim)
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_30_1.png

Le graphique ci-dessus devrait montrer clairement comment la forme et la structure du littoral changent considérablement entre la marée basse et la marée haute, car les vasières basses sont rapidement inondées par l’augmentation du niveau de l’eau.

Extraire les contours de profondeur à partir d’images

Nous souhaitons maintenant extraire une limite précise entre la terre et l’eau pour chacun des intervalles de marée ci-dessus. Le code ci-dessous identifie les contours de profondeur en fonction de la limite entre la terre et l’eau en traçant une ligne le long des pixels avec une valeur d’indice d’eau de « 0 » (la limite entre la terre et l’eau). Il renvoie un « geopandas.GeoDataFrame » avec un contour de profondeur pour chaque intervalle de marée étiqueté avec les hauteurs de marée en mètres par rapport au niveau moyen de la mer.

[15]:

# Set up attributes to assign to each waterline
attribute_df = pd.DataFrame({'tide_m': s2_intervals.tide_m.values})

# Extract waterlines
contours_gdf = subpixel_contours(da=s2_intervals.NDWI,
                                 z_values=0,
                                 crs=s2_ds.geobox.crs,
                                 attribute_df=attribute_df,
                                 min_vertices=20,
                                 dim='tide_interval')

# Plot output shapefile over the top of the first tidal interval water index
fig, ax = plt.subplots(1, 1, figsize=(15, 10))
s2_intervals.NDWI.sel(tide_interval=1).plot(ax=ax,
                                                 cmap='Greys',
                                                 add_colorbar=False)
contours_gdf.plot(ax=ax, column='tide_m', cmap='YlOrRd', legend=True)
plt.show()

Operating in single z-value, multiple arrays mode
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_33_1.png

Le graphique ci-dessus est une visualisation basique des contours de profondeur renvoyés par la fonction « subpixel_contours ». Les contours plus profonds (en m par rapport au niveau moyen de la mer) sont colorés en jaune ; les contours moins profonds sont colorés en rouge. Maintenant que nous avons le fichier de formes, nous pouvons utiliser une fonction plus complexe pour créer un graphique interactif permettant de visualiser la topographie de la zone intertidale.

Tracer une carte interactive des courbes de profondeur colorées en fonction du temps

La cellule suivante fournit une carte interactive avec une superposition des courbes de profondeur identifiées dans la cellule précédente. Exécutez-la pour afficher la carte.

Zoomez sur la carte ci-dessous pour explorer l’ensemble des courbes de profondeur obtenues. Grâce à ces données, nous pouvons facilement identifier les zones du littoral qui ne sont exposées qu’à marée basse ou d’autres zones qui ne sont recouvertes d’eau qu’à marée haute.

[16]:

map_shapefile(gdf=contours_gdf, attribute='tide_m', continuous=True)

Bien que les contours ci-dessus fournissent des informations précieuses sur la topographie de la zone intertidale, nous pouvons extraire des informations supplémentaires sur la structure 3D du littoral en les convertissant en un raster d’élévation (c’est-à-dire un modèle numérique d’élévation ou DEM).

Dans la cellule ci-dessous, nous convertissons le fichier de formes ci-dessus en un tableau de points avec des coordonnées X, Y et Z, où la coordonnée Z est l’élévation du point par rapport au niveau moyen de la mer. Nous utilisons ensuite ces points XYZ pour interpoler des élévations lisses et continues sur la zone intertidale à l’aide d’une interpolation linéaire.

[17]:

# First convert our contours shapefile into an array of XYZ points
xyz_array = contours_to_arrays(contours_gdf, 'tide_m')

# Interpolate these XYZ points over the spatial extent of the Sentinel-2 dataset
intertidal_dem = interpolate_2d(ds=s2_intervals,
                                x_coords=xyz_array[:, 0],
                                y_coords=xyz_array[:, 1],
                                z_coords=xyz_array[:, 2])

# Plot the output
intertidal_dem.plot(cmap='viridis', size=8, robust=True)
plt.show()
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_38_0.png

Vous pouvez voir dans le résultat ci-dessus que nos résultats d’interpolation sont très confus. Cela est dû au fait que l’interpolation s’étend à des zones de notre zone d’étude qui ne sont pas affectées par les marées (par exemple, des zones d’eau situées au-delà de la marée la plus basse observée et sur terre). Pour nettoyer les données, nous pouvons restreindre le MNT à la seule zone comprise entre les marées les plus basses et les plus hautes observées :

[18]:

# Identify areas that are always wet (e.g. below low tide), or always dry
above_lowest = s2_intervals.isel(tide_interval=0).NDWI < 0
below_highest = s2_intervals.isel(tide_interval=-1).NDWI > 0

# Keep only pixels between high and low tide
intertidal_dem_clean = intertidal_dem.where(above_lowest & below_highest)

# Plot the cleaned dataset
intertidal_dem_clean.plot(cmap='viridis', size=8, robust=True)
plt.show()
../../../_images/sandbox_notebooks_Real_world_examples_Intertidal_elevation_S2_40_0.png

Exporter le DEM intertidal au format GeoTIFF

En guise d’étape finale, nous pouvons prendre le DEM intertidal que nous avons créé et l’exporter sous forme de GeoTIFF qui peut être chargé dans un logiciel SIG comme QGIS ou ArcMap (pour télécharger l’ensemble de données depuis le Sandbox, localisez-le dans le navigateur de fichiers à gauche, faites un clic droit sur le fichier et sélectionnez « Télécharger »).

[19]:

# Export as a GeoTIFF
write_cog(intertidal_dem_clean,
          fname='intertidal_dem_s2.tif',
          overwrite=True)

WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option WIDTH
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option HEIGHT
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option COUNT
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option DTYPE
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option CRS
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option TRANSFORM
WARNING:rasterio._env:CPLE_NotSupported in driver GTiff does not support creation option NODATA

[19]:

PosixPath('intertidal_dem_s2.tif')

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

[21]:

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

[21]:

'2025-01-16'