Vincent GODARD - V1- 07/04/2022
Inspiré de :
Sources :
Dossier compressé à télécharger => ici.
Le jeu de données a été téléchargé lors d'un précédent TD (tdpy21tel)¶
Pour ce TD, une fenêtre a été découpée autour de l'île d'Yeu. Les données sont du 23 septembre 2021 (T30TWS_20210923T110659_TCI_b?.jp2).
Pour rendre ce TD indépendant, ces deux premiers chapitres sont identiques à ceux du TD "tdpy22tel".
# Installation de la bibliothèque rasterio qui n'est généralement pas installée en native.
!pip install rasterio
## Import des modules nécessaires dans ce TD.
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from pathlib import Path
Charger préalablement ces images dans la zone dédiée à cet effet si vous êtes sur une plateforme comme Google Colab (qui ne semble accepter que du jpg en direct ou du JPEG2000 en passant par Google Drive), Jupyter Lab, Kaggle... sinon adaptez le chemin d'accès si vous êtes sur votre disque dur.
Pour importer les fichiers jp2 (JPEG2000), charger les fichiers dans Google Drive, puis les installer dans la zone des données ("Fichiers").
from google.colab import drive
drive.mount('/content/drive')
Utiliser la commande "Copier le chemin d'accès" pour assurer un accès "sans stress" aux images avec Rasterio.
# Ouverture des 4 fichiers comme on l'a déjà fait
b2 = rasterio.open('/content/drive/MyDrive/yeu_b2.jp2')
b3 = rasterio.open('/content/drive/MyDrive/yeu_b3.jp2')
b4 = rasterio.open('/content/drive/MyDrive/yeu_b4.jp2')
b8 = rasterio.open('/content/drive/MyDrive/yeu_b8.jp2')
Nous allons voir au chapitre suivant, la description d'une méthode automatisée d'ouverture !
Une fonction est un morceau de code qui reçoit des arguments en entrée, effectue un calcul, puis génère une sortie. Elle peut être utilisée plusieurs fois dans un même script ou appelée dans d'autre scripts.
# Quelques remarques :
# Création d'une fonction qui contient un boucle de lecture de certains fichiers (qui contiennent yeu_ et finissent par jp2 !)
# Comme les images Sentinel contiennent "yeu_.jp2", ils seront seuls sélectionnés !
# La lettre "f" au début de la chaîne "next(path.glob(f'yeu_{band}.jp2'))" indique qu'il s'agit d'une chaîne formatée
# ceci permet d'inclure entre parenthèses des commandes Python pour : afficher un résultat.
def load_sentinel_image(img_folder, bands): # nom de la fonction
image = {}
path = Path(img_folder)
for band in bands: # début de la boucle
file = next(path.glob(f'yeu_{band}.jp2')) # ici la chaîne formatée
print(f'Opening file {file}')
ds = rasterio.open(file) # création d'un dataset (jeu de données) avec rasterio.open()
image.update({band: ds.read(1)})
return image
# Un dictionnaire Python est un objet qui peut stocker des paires de clés (keys) et de valeurs (values).
# Les clés peuvent correspondre à n'importe quel type de données Python "interne" (c'est-à-dire une chaîne, des nombres ou des tuples)
# Les valeurs peuvent être de n'importe quel type/objet arbitraire (ici une image).
# Utilisation de la fonction précédente pour compléter le dictionnaire des "clées" avec les bandes chargées de notre image
img = load_sentinel_image('/content/drive/MyDrive', ['b2', 'b3', 'b4', 'b8'])
img.keys()
# Interrogation de la taille lignes x colonnes de nos images
img['b2'].shape, img['b3'].shape, img['b4'].shape, img['b8'].shape
# Description de nos variables
type(img), type(img['b8'])
# Contrôle par affichage de nos variables
fig, ax = plt.subplots(1, 4, figsize=(20, 5))
ax[0].imshow(img['b2'], cmap='Blues_r')
ax[1].imshow(img['b3'], cmap='Greens_r')
ax[2].imshow(img['b4'], cmap='Reds_r')
ax[3].imshow(img['b8'], cmap='gray_r')
Ici, l'objectif est de trouver un indice pour séparer les terres immergées des terres émergées !
# Comme il y a un risque de division par "0" si le pixel du dénominateur comporte cette valeurs :
# Les deux premières lignes de code vérifient cette condition !
# "==" vérifie si la valeur du pixel est égale à "0"
# "&" signifie "AND" (ET), qui doit être strictement équivalent
# De même, si un même pixel est égal à "0" dans les deux bandes alors ce n'est pas un nombre (NaN, not a number).
def normalized_difference(img, B1, B2, eps=0.0001):
Band1 = np.where((img[B1]==0) & (img[B1]==0), np.nan, img[B1])
Band2 = np.where((img[B1]==0) & (img[B2]==0), np.nan, img[B2])
return (Band1 - Band2) / (Band1 + Band2)
# Calcul de l'IVN
ndvi = normalized_difference(img, 'b8', 'b4')
# Calcul du MNDVI (Modified Normalized Difference Water Index)
mndwi = normalized_difference(img, 'b3', 'b8')
# Contrôle des images
fig, ax = plt.subplots(1, 2, figsize=(18, 9))
ax[0].imshow(ndvi, cmap='Greens')
ax[1].imshow(mndwi, cmap='Blues')
# Création du masque des valeurs de MNDWI supérieures à -0.0
water_mask = mndwi > -0.0
# Affichage d'un tableau indiquant si la condition est remplie
water_mask
# Affichage de l'image du masque
plt.figure(figsize=(18,9))
plt.imshow(water_mask)
# Vous remarquerez que quelques pixels dans les terres sont en eau !
Pourriez-vous faire la même manipulation avec le NDVI ?
Que pensez-vous du résultat ?