Turbidité dans les zones humides

Produits utilisés : ls8_sr, ls9_sr, wofs_ls

**Keywords**: :index:`data used; landsat 9`, :index:`data used; landsat 8`, :index:`data used; WOfS`, :index:`NDTI`, :index:`wetlands`

Aperçu

Dans ce cahier, nous calculons la turbidité dans les zones humides grâce à l’indice de turbidité par différence normalisée (NDTI). La turbidité fait référence à la clarté optique de l’eau. Bien qu’elle ne soit pas strictement corrélée aux matières en suspension, à la teneur en chlorophylle ou à d’autres paramètres, elle est souvent utilisée conjointement pour fournir une analyse plus globale de la qualité de l’eau.

Le NDTI est un indice de bande défini comme :

\(\text{NDTI} = \frac{\text{Rouge} - \text{Vert}}{\text{Rouge} + \text{Vert}}\)

Lacaux et al, 2007

Ici, nous étudions le NDTI sur les zones humides dans une série chronologique de plusieurs années, à la fois par le biais d’une animation saisonnière et d’un graphique de séries chronologiques.

Description

Charger les packages et se connecter au datacube
Définir la zone d’intérêt et la plage horaire
Masquer les zones avec de l’eau en utilisant WOfS
Tracer et animer
Étudier l’étendue des zones humides au fil du temps

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 numpy as np
import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt
import rioxarray
from IPython.display import Image
from geopandas import GeoDataFrame
from odc.geo.geom import Geometry

from deafrica_tools.plotting import display_map, rgb, xr_animation, map_shapefile
from deafrica_tools.spatial import xr_vectorize, xr_rasterize
from deafrica_tools.datahandling import load_ard, mostcommon_crs, wofs_fuser
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.areaofinterest import define_area

Se connecter au datacube

[2]:

dc = datacube.Datacube(app="Wetland_turbidity")

Définir la zone d’intérêt

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 permettra de démontrer le fonctionnement de l’analyse et de fournir des résultats significatifs.

[3]:

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=13.83408, lon=-16.56464 , lat_buffer =0.2025, lon_buffer = 0.1025)

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

# Time period
time_range = ("2017-01-01", "2021-12-31")

display_map(x=lon_range, y=lat_range)

[3]:

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

[4]:

# # Define the area of interest
# lat = 13.83408
# lon = -16.56464
# lat_buffer = 0.2025
# lon_buffer = 0.1025

# # Time period
# time_range = ("2017-01-01", "2021-12-31")

# # Combine central lat,lon with buffer to get area of interest
# lat_range = (lat - lat_buffer, lat + lat_buffer)
# lon_range = (lon - lon_buffer, lon + lon_buffer)

# display_map(x=lon_range, y=lat_range)

Charger les données satellite

La première étape de l’analyse consiste à charger les données Landsat 8 et 9 pour la zone d’intérêt spécifiée. Cela utilise la fonction utilitaire prédéfinie load_ard. La fonction load_ard est utilisée ici pour charger un ensemble de données prêt à l’analyse, exempt d’ombres et de données manquantes.

Notez que la spécification du paramètre « min_gooddata=0.9 » réduit les données du milieu de l’année, lorsque cette zone connaît un temps nuageux maximal.

[5]:

query = {
    "x": lon_range,
    "y": lat_range,
    "time": time_range,
    "resolution": (-30, 30),
    "align": (15, 15),
}

output_crs = mostcommon_crs(dc=dc, product='ls8_sr', query=query)

ds = load_ard(
    dc=dc,
    products=['ls8_sr', 'ls9_sr'],
    measurements=["red", "green", "blue", "nir"],
    group_by="solar_day",
    dask_chunks={"time": 1, "x": 2000, "y": 2000},
    output_crs=output_crs,
    min_gooddata=0.9,
    **query
)

/opt/venv/lib/python3.12/site-packages/deafrica_tools/datahandling.py:245: 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 USGS Collection 2
Finding datasets
    ls8_sr
    ls9_sr
Counting good quality pixels for each time step
Filtering to 39 out of 114 time steps with at least 90.0% good quality pixels
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Returning 39 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])

Charger WOfS

Les valeurs de NDTI sont relatives et l’indice ne doit être appliqué qu’à l’eau. L’inclusion de terres dans la zone d’intérêt faussera la carte des couleurs. Le WOfS temporel est utilisé ici pour définir un seuil pour les pixels des zones susceptibles d’être de l’eau.

[6]:

water = dc.load(
            product="wofs_ls_summary_alltime",
            like=ds.geobox,
            dask_chunks={"time": 1, "x": 2000, "y": 2000}
            )
#extract from mask the areas classified as water
water_extent = (water.frequency > 0.3).squeeze()

print(water_extent)

<xarray.DataArray 'frequency' (y: 1499, x: 750)> Size: 1MB
dask.array<getitem, shape=(1499, 750), dtype=bool, chunksize=(1499, 750), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 12kB 1.552e+06 1.552e+06 ... 1.508e+06 1.507e+06
  * x            (x) float64 6kB 3.196e+05 3.197e+05 ... 3.421e+05 3.421e+05
    time         datetime64[ns] 8B 2004-12-31T11:59:59.999999
    spatial_ref  int32 4B 32628
Attributes:
    units:         1
    nodata:        nan
    crs:           PROJCRS["WGS 84 / UTM zone 28N",BASEGEOGCRS["WGS 84",ENSEM...
    grid_mapping:  spatial_ref

/tmp/ipykernel_2589/3674913840.py:3: DeprecationWarning: Geobox extraction logic has moved to odc-geo and the .geobox property is now deprecated. Please access via .odc.geobox instead.
  like=ds.geobox,

Inspecter l’étendue de l’eau

[ ]:

# Plot the geomedian composite and water extent
fig, ax = plt.subplots(1, 2, figsize=(12, 6))

#plot the true colour image
rgb(ds.isel(time=5), ax=ax[0])

#plot the water extent from WOfS
water_extent.plot.imshow(ax=ax[1], cmap="Blues", add_colorbar=False)

# Titles
ax[0].set_title(""), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
ax[1].set_title("Water Extent"), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False);

Parcelle NDTI

Ici, nous calculons le NDTI et l’ajoutons à l’ensemble de données chargé. Il génère également une animation qui est exportée vers le sandbox et peut être enregistrée en cliquant avec le bouton droit sur le « .gif » et en sélectionnant « Télécharger ».

[ ]:

ds = calculate_indices(ds.where(water_extent), index="NDTI", satellite_mission="ls")

[ ]:

fig, ax = plt.subplots(1, 2, figsize=(18, 12))
cmap='BrBG_r'
vmin=-0.4
vmax=0.4

ds.NDTI.isel(time=0).plot(ax=ax[0], cmap=cmap, add_colorbar=False, vmin=vmin, vmax=vmax)
ds.NDTI.isel(time=-1).plot(ax=ax[1], cmap=cmap, add_colorbar=True, vmin=vmin, vmax=vmax);

[ ]:

# Produce time series animation of NDTI
xr_animation(ds=ds,
             output_path='NDTI_'+time_range[0]+'_'+time_range[1]+'.gif',
             bands=['NDTI'],
             show_text='NDTI',
             interval=250,
             imshow_kwargs={'cmap': 'BrBG_r', 'vmin': vmin, 'vmax': vmax},
             width_pixels=300)

# Plot animated gif
plt.close()
Image(filename='NDTI_'+time_range[0]+'_'+time_range[1]+'.gif')

Étudier la taille des plans d’eau

Les séries chronologiques animées et les graphiques ci-dessus montrent une nette différence dans le NDTI entre les grands cours d’eau, qui apparaissent généralement en bleu-blanc, et les plans d’eau plus petits, qui apparaissent plutôt orange à marron dans la palette de couleurs.

Nous pouvons émettre l’hypothèse que les plans d’eau plus petits ont généralement des valeurs NDTI plus élevées que les plus grands, ce qui signifie qu’ils sont plus troubles. L’histogramme ci-dessous suggère que cela pourrait être le cas car il y a deux pics évidents. Nous pouvons tester cela plus en détail en divisant les plans d’eau par taille ci-dessous.

[ ]:

plt.hist(ds.NDTI.to_numpy().flatten(), color = 'blue', edgecolor = 'black',
         bins = 8000)
plt.xlim((-1,1))

# Add labels
plt.xlabel('NDTI')
plt.ylabel('Count of pixels')

Vectoriser les plans d’eau

Nous allons diviser les plans d’eau en deux classes de taille, petite et grande, à un seuil de 50 000 mètres carrés. Vous pouvez ajuster le seuil et le tester dans les graphiques ci-dessous.

[ ]:

polygons = xr_vectorize(water_extent, crs='epsg:32628', mask=water_extent==True)
polygons['area'] = polygons.area
polygons['size'] = np.where(polygons['area']<50000, 1, 2) # se

Le graphique interactif ci-dessous montre les grands plans d’eau en jaune et les plus petits en violet.

[ ]:

map_shapefile(polygons, attribute='size')

Ici, nous rastérisons les polygones et masquons le NDTI pour les petits et les grands plans d’eau, puis calculons la moyenne et l’écart type des séries chronologiques. Cela nous permet de tracer la série chronologique dans la cellule ci-dessous. Le tracé confirme que les petits plans d’eau ont des valeurs NDTI plus élevées que les plus grands, et que les deux classes de taille suivent le même schéma au fil du temps.

[ ]:

water_size = xr_rasterize(gdf=polygons, da=water, attribute_col='size')

ds_large = ds.where(water_size==2)
ds_small = ds.where(water_size==1)

NDTI_ts_large = ds_large.NDTI.mean(['x','y'])
NDTI_ts_small = ds_small.NDTI.mean(['x','y'])
NDTI_ts_large_std = ds_large.NDTI.std(['x','y'])
NDTI_ts_small_std = ds_small.NDTI.std(['x','y'])

[ ]:

fig, ax = plt.subplots(1, 1, figsize=(11, 5))
NDTI_ts_large.plot(ax=ax, label='Large waterbodies', linestyle='dashed', marker='o')
NDTI_ts_small.plot(ax=ax, label='Small waterbodies', linestyle='dashed', marker='o')
ax.fill_between(
    NDTI_ts_large.time,
    NDTI_ts_large-NDTI_ts_large_std,
    NDTI_ts_large+NDTI_ts_large_std,
    alpha=0.2,
)
ax.fill_between(
    NDTI_ts_small.time,
    NDTI_ts_small-NDTI_ts_small_std,
    NDTI_ts_small+NDTI_ts_small_std,
    alpha=0.2,
)
plt.legend(loc="upper left")
plt.title("")
plt.xlabel("");

Enfin, nous pourrions émettre l’hypothèse que la variation saisonnière de la turbidité observée ci-dessus est liée aux précipitations. Le graphique des précipitations sur la même période, ci-dessous, suggère que c’est le cas. Cependant, nous devrions effectuer une analyse plus approfondie pour déterminer les conditions qui varient en fonction de la turbidité.

[ ]:

ds_rf_month = dc.load(product='rainfall_chirps_monthly',
                resolution=(-5000, 5000),
                output_crs='epsg:6933',
                x = lon_range,
                y = lat_range,
                time = time_range)

ds_rf_month.resample(time='1M').mean().rainfall.mean(['x','y']).plot(figsize=(11, 5))
plt.title("");

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

[ ]:

print(datacube.__version__)

Dernier test :

[ ]:

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