Exercice: déterminer l’étendue des plans d’eau

Aperçu

In this exercise, we will create a new notebook to determine the extent of water bodies using the Water Observation from Space (WOfS) product. The WOfS product uses an automated water mapping algorithm to identify water in Landsat 8 images.

Le cahier comprendra les étapes suivantes:

• Charger le produit de couche d’entités WOfS et les données Landsat 8

• Identifier les pixels d’eau de WOfS

• Tracez une image en couleurs vraies à l’aide des données Landsat

• Tracez la superficie du plan d’eau pour la même zone à l’aide des données WOfS

• Personnalisez les tracés

À l’issue de cet exercice, vous serez en mesure de déterminer l’étendue des plans d’eau à l’aide du produit WOfS.

Configurer le notebook

Dans votre dossier Training, créez un nouveau notebook Python 3. Nommez-le Water_extent_exercise.ipynb. Pour plus d’instructions sur la création d’un nouveau notebook, consultez les instructions de la session 2.

Charger des packages et des fonctions

Dans la première cellule, tapez le code suivant, puis exécutez la cellule pour importer les dépendances Python nécessaires.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import datacube

from deafrica_tools.datahandling import load_ard, wofs_fuser
from deafrica_tools.plotting import display_map, rgb

Dans cet exercice, nous importons une nouvelle fonction, wofs_fuser`'. ``wofs_fuser assure que les données WOfS entre les scènes sont combinées correctement.

Connectez-vous au datacube

Entrez le code suivant et exécutez la cellule pour créer notre objet dc, qui donne accès au datacube.

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

Sélectionnez la zone d’intérêt

Dans la cellule suivante, saisissez le code suivant, puis exécutez-le pour sélectionner une zone et une heure. Dans cet exercice, nous utilisons un point central et un tampon pour définir notre zone d’intérêt, de manière similaire à ce que nous avons fait dans l’exercice Session 5.

La seule différence ici est que nous fournissons un tampon de latitude et un tampon de longitude. Dans cet exemple, ils ont la même valeur. Cependant, le choix de différentes valeurs de tampon vous permet de sélectionner des zones d’intérêt rectangulaires plutôt que des carrés.

# Define the area of interest
lat = -6.0873
lon = 35.1817

lat_buffer = 0.2
lon_buffer = 0.2

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

# Define the year
time = '2018'

Note

Rappelons que les lignes de code commençant par # sont des commentaires. Ils aident à expliquer le code et peuvent être supprimés ou ajoutés sans avoir d’impact sur les scripts Python.

Dans la cellule suivante, entrez le code suivant, puis exécutez-le pour afficher la zone sur une carte. Puisque nous avons défini notre zone en utilisant les variables lon_range et lat_range, nous pouvons les utiliser au lieu de retaper (lat - lat_buffer, lat + lat_buffer) et (lon - lon_buffer, lon + lon_buffer).

display_map(x=lon_range, y=lat_range)

Créer un objet de requête

Remarquez que lat_range, lon_range et time ont tous été définis dans la cellule précédente, donc nous pouvons les utiliser comme variables ici. Nous allons les utiliser pour créer une requête.

La variable ``query`”” ci-dessous est un dictionnaire Python. Elle peut être utilisée pour stocker des paramètres. La création d’une variable objet telle que ``query`” permet de réutiliser des paramètres dans différentes fonctions qui acceptent les mêmes paramètres d’entrée.

Cela nous est utile car nous pouvons l’utiliser pour charger les données Landsat 8, puis l’utiliser à nouveau pour charger les données WOfS.

Dans la cellule suivante, entrez le code suivant, puis exécutez-le.

query = {
    'x': lon_range,
    'y': lat_range,
    'resolution': (-30, 30),
    'output_crs':'EPSG:6933',
    'group_by': 'solar_day',
    'time': (time),
}

Note

Remarquez que la structure du dictionnaire query est légèrement différente de dc.load ou load_ard. Chaque nom de paramètre est entre guillemets '' et est suivi de deux points :.

Charger des données

Dans la cellule suivante, nous chargeons les jeux de données Landsat et WOfS, en les nommant respectivement ds_landsat et ds_wofs.

Dans les fonctions ci-dessous, nous pouvons passer directement l’objet query en utilisant **query - cela donnera tous les paramètres définis dans query à la fonction.

Le principal avantage est que nous pouvons utiliser la même requête query pour Landsat 8 et WOfS, ce qui nous évite de la retaper et de faire des erreurs.

Charger Landsat 8

Pour Landsat 8, nous pouvons utiliser la fonction load_ard.

ds_landsat = load_ard(dc=dc,
                     products=['ls8_sr'],
                     measurements=['red', 'green', 'blue'],
                     **query)

ds_landsat

Note some users will find ds_landsat loads with 43 timesteps. That will not affect the next steps.

Charger WOfS

Pour WOfS, nous devons utiliser la fonction dc.load.

ds_wofs = dc.load(product="wofs_ls",
         fuse_func=wofs_fuser,
         **query
        )

ds_wofs

Note

As of September 2021, WOfS can best be accessed through the product name wofs_ls. This replaces the deprecated ga_ls8c_wofs_2 product name. For more information on wofs_ls, visit the WOfS technical specifications.

Calcul de l’étendue de l’eau

Comprendre les couches d’entités WOfS

Les couches d’entités WOfS sont stockées sous forme de nombre binaire, où chaque chiffre du nombre est défini indépendamment ou non en fonction de la présence (1) ou de l’absence (0) d’une entité particulière. Vous trouverez ci-dessous une ventilation de la valeur décimale représentant les fonctionnalités.

Attribut	Valeur décimale
Pas de données	1
Non contigu	2
Mer	4
Terrain ou faible angle solaire	8
Pente élevée	16
Ombre de nuage	32
Nuage	64
L’eau	128

Par exemple, une valeur de 128 indique que de l’eau a été observée pour le pixel, tandis qu’une valeur de 32 indiquerait une ombre de nuage.

Dans la cellule suivante, nous allons extraire uniquement les caractéristiques de l’eau à partir des données WOfS. Cela se fait en trouvant les valeurs où la mesure water est égale à 128. En Python, nous pouvons trouver quels pixels ont une valeur de 128 en utilisant l’expression == :

Extraire les pixels d’eau

ds_valid_water = ds_wofs.water == 128

Le tableau ds_valid_water ne contient pas les valeurs décimales des couches de caractéristiques WOfS. Au lieu de cela, il a la valeur False si le pixel n’est pas de l’eau, et True s’il est de l’eau. Vous pouvez le vérifier en visualisant le DataArray ds_valid_water.

ds_valid_water

Calculez la surface par pixel

Le nombre de pixels peut être utilisé pour la superficie de la masse d’eau si la superficie du pixel est connue. Nous pouvons extraire la taille d’un pixel à partir du paramètre resolution de notre requête query, puis diviser la surface d’un seul pixel (en mètres carrés) par le nombre de mètres carrés dans un kilomètre carré.

En Python, number**2 renvoie la valeur au carré de number.

pixel_length = query["resolution"][1]  # in metres
m_per_km = 1000  # conversion from metres to kilometres
area_per_pixel = pixel_length**2 / m_per_km**2

Calculer la surface des pixels d’eau

Maintenant que nous savons combien de surface est couverte par un pixel, nous pouvons compter le nombre de pixels d’eau et le multiplier par cette valeur pour obtenir la surface totale couverte par l’eau.

Comme nous l’avons vu ci-dessus, le tableau ds_valid_water contient des valeurs True pour les pixels d’eau, et False sinon. Lorsque nous utilisons la fonction .sum, elle compte les valeurs True comme 1, et False comme 0. Par conséquent, la somme sera égale au nombre total de pixels d’eau pour ce pas de temps.

Ci-dessous, nous définissons les dimensions x et y pour être sûrs de faire la somme de tous les pixels sur les dimensions spatiales. Cela signifie que nous obtenons une somme de pixels pour chaque pas de temps. Cela nous permettra de suivre l’évolution de la surface de l’eau dans le temps.

ds_valid_water_pixel_sum = ds_valid_water.sum(dim=['x', 'y'])
ds_valid_water_area = ds_valid_water_pixel_sum * area_per_pixel

Tracer des séries chronologiques

Maintenant que nous avons la superficie de l’eau dans chaque observation, nous pouvons tracer une série chronologique pour nous aider à identifier les dates où il y avait plus ou moins d’eau dans la zone d’intérêt.

Ci-dessous, il y a un code pour configurer, afficher et personnaliser le tracé. Les paramètres sont les suivants:

• plt.figure(figsize=(18, 4)) : crée un objet figure pour contenir le graphique, et le rend long de 18 pouces et haut de 4 pouces.

• ds_valid_water_area.plot(marker='o', color='#9467bd') : trace les données de la zone d’eau avec des marqueurs circulaires en violet (code couleur HEX #9467bd)

• plt.title('Analyse des séries temporelles de la surface observée en eau') : Donnez un titre au graphique

• plt.xlabel('Dates') : Étiquette l’axe des x

• plt.ylabel('Waterbody area (km$^2$)') : Étiquette l’axe des ordonnées. Les symboles $ permettent l’utilisation de LaTeX, un langage de composition mathématique

• Plt.tight_layout()`` : Formate l’image de façon à ce que tous les axes soient clairement visibles.

plt.figure(figsize=(18, 4))
ds_valid_water_area.plot(marker='o', color='#9467bd')
plt.title('Time Series Analysis of water observed area')
plt.xlabel('Dates')
plt.ylabel('Waterbody area (km$^2$)')
plt.tight_layout()

Affichage de la couverture d’eau pour un pas de temps sélectionné

À partir du graphique ci-dessus, vous pouvez choisir n’importe quel pas de temps (entre 0 et 45) pour afficher le résultat sur WOfS et Landsat 8.

Par exemple, regardons le cinquième pas de temps, timestep = 4.

timestep = 4

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

#plot the true colour image
ds_nearest_landsat = ds_landsat.sel(time=ds_wofs.time.isel(time=timestep), method='nearest')
rgb(ds_nearest_landsat, ax=ax[0])

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

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

plt.show()

Ce code utilise des paramètres supplémentaires pour personnaliser le tracé, notamment en permettant d’avoir deux tracés ensemble. Si vous souhaitez en savoir plus sur la réalisation de ce type d’intrigue, veuillez demander aux instructeurs lors d’une session en direct.

Essayez différentes valeurs de timestep - pouvez-vous trouver une image où le lac est asséché ?

Conclusion

Toutes nos félicitations! Vous avez créé votre propre carnet d’étendue d’eau. Il est comparable au «cahier d’étendue d’eau Sandbox existant <https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/master/Real_world_examples/Water_extent.ipynb>« __.

Vous avez maintenant construit votre deuxième étude de cas! Vous avez envie de réfléchir à ce qui était similaire et différent entre les deux. Y a-t-il des morceaux de code que vous pourriez réutiliser pour une nouvelle analyse? Comment pourriez-vous modifier vos études de cas pour faire des analyses plus complexes?

Si vous souhaitez expérimenter davantage, essayez d’exécuter le code avec différentes zones. Avez-vous appris quelque chose d’intéressant à partager avec nous?

Activité optionnelle

Si vous êtes curieux de savoir comment fonctionne l’étude de cas existante, vous pouvez l’ouvrir et l’exécuter dans le bac à sable:

1. From the main Sandbox folder, open the Real_world_examples folder

2. Double-click the Water_extent.ipynb notebook to open it

Le cahier a déjà été exécuté, vous pouvez donc le lire étape par étape. Cependant, vous pouvez trouver utile d’effacer les sorties et d’exécuter chaque cellule étape par étape pour voir comment elle fonctionne. Vous pouvez le faire en cliquant sur Kernel -> Restart Kernel and Clear All Outputs. Lorsqu’on vous demande si vous voulez redémarrer le noyau, cliquez sur Restart.

Note

Si vous voulez le modifier de manière significative, nous vous recommandons d’en faire une copie, comme vous l’avez fait dans Session 1.

Il existe de nombreuses similitudes entre le bloc-notes que vous avez créé dans cette session et le bloc-notes Sandbox existant. Notez peut-être ce qui est similaire et ce qui est différent. Si vous avez des questions sur le fonctionnement du bloc-notes existant, veuillez demander aux instructeurs lors d’une session en direct.