Modélisation de l’élévation intertidale à l’aide de données de marée
Products used: s2_l2a_c1
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> » __
Keywords: data used; Sentinel-2 c1, water; tide modelling, water; waterline extraction, band index; NDWI, data methods; compositing
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 :
Load in a cloud-free Sentinel-2 c1 time series.
Calculer un indice hydrique (NDWI).
Étiquetez et triez les images satellites par hauteur de marée.
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).
Extraire et visualiser la topographie de la zone intertidale sous forme de courbes de profondeur.
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 eo_tides import tag_tides
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-6dcbefde-3320-11f1-8eb3-66c27c7d89dc
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/42279/status |
Cluster Info
LocalCluster
50410277
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/42279/status | Workers: 1 |
| Total threads: 4 | Total memory: 26.21 GiB |
| Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-8b984144-2c4f-4643-a7aa-e6a95f7f8955
| Comm: tcp://127.0.0.1:43311 | Workers: 0 |
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/42279/status | Total threads: 0 |
| Started: Just now | Total memory: 0 B |
Workers
Worker: 0
| Comm: tcp://127.0.0.1:39415 | Total threads: 4 |
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/33577/status | Memory: 26.21 GiB |
| Nanny: tcp://127.0.0.1:32783 | |
| Local directory: /tmp/dask-scratch-space/worker-4ad2_w5y | |
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 :
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.
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
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', '2021')
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]:
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_c1'],
output_crs=output_crs,
min_gooddata=0.5,
group_by='solar_day',
dask_chunks={},
**query)
Using pixel quality parameters for Sentinel 2
Finding datasets
s2_l2a_c1
Counting good quality pixels for each time step
/opt/venv/lib/python3.12/site-packages/rasterio/warp.py:385: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
dest = _reproject(
Filtering to 354 out of 435 time steps with at least 50.0% good quality pixels
Applying pixel quality/cloud mask
Re-scaling Sentinel-2 C1 data
Returning 354 time steps as a dask array
/opt/venv/lib/python3.12/site-packages/deafrica_tools/datahandling.py:565: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
ds = xr.merge([ds_data, ds_masks])
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: 4GB
Dimensions: (time: 354, y: 889, x: 845)
Coordinates:
* time (time) datetime64[ns] 3kB 2017-11-28T07:24:55.563000 ... 202...
* 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 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
green (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
blue (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
nir (time, y, x) float32 1GB 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)
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 :
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()
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.
The tag_tides function below uses the EOT20 model (by default) to calculate the height of the tide at the exact moment each satellite image in our dataset was taken, and adds this as a new tide_m attribute in our dataset (for more information about this function, refer to the Tidal modelling notebook).
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
model_dir="/var/share/tide_models"
s2_ds_tidal = tag_tides(data=s2_ds,directory=model_dir, tidepost_lat=None, tidepost_lon=None)
# Print the output dataset with new `tide_m` variable
print(s2_ds_tidal)
Setting tide modelling location from dataset centroid: 44.32, -19.26
Modelling tides with EOT20
<xarray.DataArray 'tide_height' (time: 354)> Size: 1kB
array([ 0.39306295, -0.7169483 , 0.8547459 , -1.0570892 , -1.3540006 ,
0.7487494 , 0.4216002 , -0.92028236, 0.27469084, 0.40118676,
0.7122888 , 0.5972058 , -0.32074735, -1.3419393 , -0.57915896,
0.3328396 , -1.4734659 , -0.06816646, 0.84975755, 0.33673579,
-0.531434 , -1.4372492 , -1.3834155 , -0.5045415 , 0.60837865,
-0.8928216 , -1.4254504 , -1.097326 , 0.15669289, 0.45031512,
-0.10479415, -1.1515855 , -0.2946027 , 0.44667545, 0.08142319,
-1.4411024 , -1.9138846 , -0.9104993 , -0.2276676 , -1.2434506 ,
-1.6270472 , -1.1103843 , -0.05712806, 0.51531523, -0.18608418,
-1.6180627 , -1.7449201 , 0.23211198, 0.3429487 , -1.2608172 ,
-1.6091816 , -0.75132567, 0.4098541 , 0.6055012 , -0.3702263 ,
-1.5586641 , -1.4210099 , 0.38662583, 0.47273168, -0.12416134,
-1.3095279 , -1.5145022 , -0.2319423 , 0.8729848 , 0.52599424,
-0.5585511 , -1.4597821 , -1.1677879 , -0.18750156, 0.516716 ,
0.5346426 , -0.1981567 , -1.4022771 , 0.21935251, 1.0394722 ,
0.20821114, -0.80758435, -1.462378 , -1.0778812 , -0.12395494,
0.5635668 , 0.41503513, -0.5412106 , -1.8211501 , -1.264274 ,
0.4359955 , 0.85512155, -1.0660336 , -1.0875367 , -0.09157696,
0.55928326, -1.0708017 , -2.0053718 , -1.0361816 , 0.46058995,
0.55588067, -0.4739317 , -1.6123146 , -1.0606712 , 0.06124351,
...
-1.2659749 , -0.3809658 , 0.27727115, 0.2327453 , -0.40068245,
-1.6474174 , -1.7319325 , -0.4543758 , 0.4933956 , -0.03723864,
-0.94873583, -1.6224958 , -1.2880163 , -0.43266362, 0.20569777,
-0.6394708 , -1.7155238 , -0.18717477, 0.4636609 , -0.27835113,
-1.0982447 , -1.5771451 , -1.2107817 , -0.29526484, 0.35682547,
0.20724222, -0.7757151 , -1.8972222 , -1.3317802 , 0.22903933,
0.55058885, -1.0461912 , -1.4509429 , -1.0236506 , 0.03174575,
0.6888924 , 0.2818686 , -0.8768377 , -1.6536078 , -0.7667877 ,
0.5607486 , 0.6308623 , -0.25051507, -0.9947342 , -0.8143968 ,
0.41217443, 0.9662385 , -1.0852911 , -1.3899453 , -0.34714535,
0.643293 , 0.5916012 , -1.094879 , -1.4335526 , -0.6639648 ,
0.6587496 , 0.9075786 , -0.4168858 , -1.3995167 , -1.2412237 ,
-0.20715041, 0.52902275, 0.4360091 , -0.55083996, -1.3706405 ,
-1.5852262 , -0.5402477 , 0.6848246 , 0.48836878, -1.0109856 ,
-1.6433502 , -1.1732332 , -0.22480284, 0.39703104, 0.2746532 ,
-0.8251906 , -1.7078724 , -1.5996971 , -0.30079836, 0.5967017 ,
0.00732199, -1.3494126 , -1.6799028 , -0.18972693, -1.0455463 ,
-1.8898833 , -1.247205 , 0.5515938 , -0.25369942, -1.5344704 ,
-0.8549546 , 0.29055828, -1.1675752 , -1.7623956 , -0.5961458 ,
-1.2759078 , -1.325917 , -0.5647778 , 0.2719121 ], dtype=float32)
Coordinates:
* time (time) datetime64[ns] 3kB 2017-11-28T07:24:55.563000 ... 2021...
tide_model <U5 20B 'EOT20'
[11]:
# Plot the resulting tide heights for each Sentinel-2 image:
s2_ds["tide_height"] = s2_ds_tidal
s2_ds.tide_height.plot()
[11]:
[<matplotlib.lines.Line2D at 0x7f454e146420>]
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_height')
# Bin tide heights into 9 tidal intervals from low (1) to high tide (9)
binInterval = np.linspace(s2_ds.tide_height.min().values,
s2_ds.tide_height.max().values,
num=10)
tide_intervals = pd.cut(s2_ds.tide_height,
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: 5GB
Dimensions: (time: 354, y: 889, x: 845)
Coordinates:
* time (time) datetime64[ns] 3kB 2019-09-29T07:25:07.662000 ... 2...
* 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
tide_model <U5 20B 'EOT20'
Data variables:
red (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
green (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
blue (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
nir (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
NDWI (time, y, x) float32 1GB dask.array<chunksize=(1, 889, 845), meta=np.ndarray>
tide_height (time) float32 1kB -2.005 -1.919 -1.914 ... 1.039 1.06 1.067
tide_interval (time) int64 3kB 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_height.plot()
for i in binInterval: plt.axhline(i, c='black', alpha=0.5)
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_height']]
.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()
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_height': s2_intervals.tide_height.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_height', cmap='YlOrRd', legend=True)
plt.show()
/tmp/ipykernel_3763/3950071209.py:7: DeprecationWarning: Geobox extraction logic has moved to odc-geo and the .geobox property is now deprecated. Please access via .odc.geobox instead.
crs=s2_ds.geobox.crs,
Operating in single z-value, multiple arrays mode
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_height', 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_height')
# 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()
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()
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)
[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.9.13
[21]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[21]:
'2026-04-08'