- Cours de télédétection avancée pour Master Sigma.
- Introduction
- Gestion des tableaux
- D'une image au tableau
- Écriture d'une image géoréférencée
- Filtre spatial (tenant compte des voisins)
- Museo ToolBox
- Apprentissage automatique sans images
- Apprentissage automatique à partir d'images
Ce cours a pour but de vous rendre le plus autonome possible en télédétection à partir de python. À la fin de ce cours vous serez capable d'interagir avec des images directement depuis python et ce de manière optimisée et efficiente (lecture par bloc par exemple). Vous allez aussi découvrir comment manier les algorithmes d'apprentissage automatique afin d'apprendre des modèles pour cartographier par l'exemple l'occupation du sol.
Les données utilisées pour ce cours sont disponibles ici en ligne.
Les bibliothèques python utilisées dans ce cours sont :
- gdal (raster)
- osr (projection)
- ogr (vecteur)
- numpy (tableau)
- scipy (calcul scientifique)
- scikit-learn (apprentissage automatique)
- museotoolbox (facilite la lecture/écriture des rasters)
Si une des bibliothèques venait à manquer, vous pouvez l'installer en tapant la commande suivante dans le terminal :
python3 -m pip install scikit-learn --userToutes les classes, fonctions et variables que nous allons créer pendant ce cours doivent suivre le convention PEP8.
| Type | Écriture |
|---|---|
| Class | MaSuperClasse |
| Function | ma_super_fonction |
| Variable | ma_variable |
Les codes doivent être bien documentés, en suivant la docstring numpy.
/!\ Attention, un canard en plastique pourra venir relire et valider votre code.
| Ordre | Date | Horaire | Thème |
|---|---|---|---|
| 1 | Lun. 02/12/2019 | 8h/12h | Manipulation des tableaux/matrices. |
| 2 | Lun. 09/12/2019 | 8h/12h | Du tableau au raster. |
| 3 | Lun. 16/12/2019 | 13h30/17h30 | Filtre spatial (tenant compte des voisins) |
| 4 | Ven. 18/12/2019 | 8h/12h | Apprentissage automatique avec scikit-learn |
| 5 | Mer. 08/01/2020 | 8h/12h | Découverte de Museo ToolBox |
| 6 | Ven. 08/01/2020 | 13h30/17h30 | PARTIEL |
Pour gérer un tableau ou une matrice (c'est-à-dire une image sans géoréférencement)
import numpy as np
X = np.load('sentinel2_3a_20180815.npy')
X.shape # la taille de l'image (lignes, colonnes, bandes...)
>>> (200, 200, 4)
X.ndim # le nombre de dimensions
>>> 3
X.size # le nombre de cellules
>>> 160000- Comment pouvez-vous obtenir uniquement le nombre de lignes ?
- Comment pouvez-vous obtenir uniquement le nombre de bandes ?
- De combien de pixels est composée chaque bande ?
Les : permettent de sélectionner toutes les cellules de la dimension en question.
Les ... permettent de sélectionner toutes les dimensions suivantes ou précédentes.
X[0,:,0] # première ligne, toutes les colonnes, première bande
X[0,...] # première ligne, et toutes les autres dimensions
X[...,0] # toutes les lignes et colonnes de la bande 0from matplotlib import pyplot as plt
plt.imshow((X[...,0])/(np.amax(X[...,0])),cmap='gray') # Si vous voulez afficher une bande- Accéder au pixel/ à la cellule de la première ligne et première colonne
- Accéder au pixel / à la cellule de la dernière ligne et dernière colonne
- Parcourir la matrice cellule par cellule (boucle for)
- Calculer le ndvi :
Sachant que l'infra-rouge est la dernière bande (la numéro 4, donc en partant de 0 la numéro 3, et le rouge est la bande numéro 3).
Pour ceux qui ont terminé, vous pouvez chronométrer et essayer d'optimiser le temps de traitement pour calculer le NDVI.
import gdal
data_src = gdal.Open('sentinel2_3a_20180815.tif') # objet gdal
data_src.RasterCount # nombre de bandes
data_src.RasterXSize # nombre de colonnes
data_src.RasterYSize # nombre de lignes
data_src.GetRasterBand(1) # renvoie le tableau de la bande 1
data_src.GetRasterBand(2) # renvoie le tableau de la bande 2
data_src.ReadAsArray() # renvoie toute l'image en tableau- De combien de bandes, de lignes, et de colonnes est composée l'image ?
- Ouvrir l'image en objet gdal
- Calculer le NDVI pour toute l'image.
- Optimiser le calcul du NDVI pour ne plus le faire cellule par cellule (éviter la boucle for).
Maintenant que l'on sait ouvrir, lire, et faire un traitement sur une image, il ne reste plus qu'à sauvegarder le résultat dans une nouvelle image.
Pour cela quelques lignes pour définir nos besoins sont nécessaires :
driver = gdal.GetDriverByName("GTiff") # on choisi de créer un GeoTIFF
out_data = driver.Create('/tmp/mon_ndvi.tif', data_src.RasterXSize, data_src.RasterYSize, 1, gdal.GDT_Float32) # 1 pour une bande
out_data.SetGeoTransform(data_src.GetGeoTransform()) # même géotransformation que l'image d'origine
out_data.SetProjection(data_src.GetProjection()) # même projection que l'image d'origineLa variable out_data est désormais une instance de votre fichier de sortie. Vous pouvez désormais intéragir avec cette instance et sauvegarder le résultat du ndvi dans votre nouvelle image :
out_data.GetRasterBand(1).WriteArray(ndvi) # j'écris mon NDVI dans la bande 1
out_data.FlushCache() # écrit sur le disque
out_data = None # précaution, permet de bien spécifier que le fichier est ferméCréer pour chacune consignes suivantes une fonction qui vous permet :
- de lire une image (fonction retournant l'objet gdal)
- de calculer le NDVI (à partir de l'objet gdal, on demande à l'utilisateur la position de la bande rouge et infrarouge)
- d'écrire une image (on donne un objet gdal, un tableau et le nom de fichier que l'on veut écrire)
Une fois terminée, vous pourrez à partir de 3 lignes :
- ouvrir une image
- calculer le ndvi (en donnant un tableau par bande à la fonction, afin de ne presque pas modifier votre fonction précédente)
- écrire le ndvi
Comment identifier les forêts des cartes de l'État-Major ?
Des traits plus ou moins fins et plus ou moins serrés représentent les pentes dans les cartes. Nous avons donc besoin de supprimer ces rayures afind d'avoir une couleur verte homogène pour identifier la forêt.
Pour cela plusieurs filtres doivent être utilisés :
- closing filter (filtre de fermeture)
- médian (récupérer les contours)
La carte de l'État-major est à télécharger ici.
- Quelle est la valeur des pixels blancs ?
- Quelle est la valeur des pixels noirs ?
- Que fait et à quoi sert le closing filter ?
Créer une fonction qui parcourt l'image selon un nombre de voisins défini (1 = les premiers voisins (8 donc), 2 = les voisins jusqu'à 2 pixels de distance (24)...) Puis en tenant compte des voisins, appliquez :
- un filtre max
- un filtre min
- un filtre median dont on peut définir le nombre d'itération
L'idée derrière Museo ToolBox est de créer une bibliothèque python assez généraliste pour la télédétection. L'un des points les plus cruciaux et la lecture et l'écriture par bloc des images. Pour faciliter cette étape, la classe RasterMath a été développée afin que l'utilisateur n'ait pas à se soucier du traitement du raster (projection, nombre de bandes, lecture par bloc, compression de l'image, gestion des no-data et du type de données). L'utilisateur a juste besoin de créer une fonction qui traite un tableau numpy, et RasterMath se charge de tout le reste.
Pour installer Museo ToolBox, ouvrez le terminal linux (ctrl + alt + T) et coller la ligne suivante :
python3 -m pip install museotoolbox --user -ULe --user indique que vous l'installez que pour votre compte utilisateur (vous n'avez pas besoin d'avoir les droits administrateurs), le -U indique que si museotoolbox est déjà installé, alors on fera la mise à jour.
Si vous souhaitez avoir la toute dernière version en développement, coller alors la ligne suivante :
python3 -m pip install https://github.com/nkarasiak/museotoolbox/archive/develop.zip --user -U
RasterMath est un outil du module processing du Museo ToolBox. Cette classe permet d'interagir facilement avec une ou plusieurs images et d'appliquer une fonction de votre choix. La lecture et l'écriture se fera bloc par bloc (cf graphique ci-dessus), et les paramètres de création du raster seront optimisés (compression par défaut, tuilage de l'image, choix du type de données selon le type du tableau numpy...)
import museotoolbox as mtbPuis créez une instance de RasterMath avec l'image Sentinel-2, et donner un bloc de l'image choisi aléatoirement.
my_image = mtb.processing.RasterMath('sentinel2_3a_20180815.tif')
X = my_image.get_random_block()
print(X.shape)
(65536, 10)RasterMath permet d'interagir facilement avec votre image. Par défaut, cette classe vous renvoie un pixel par ligne (ici 65536), et en colonne chaque bande du raster (ici 10).
def ndvi(X,r,ir):
return np.divide(X[...,ir]-X[...,r],X[...,ir]+X[...,r])Ensuite tester votre fonction tout simplement avec le bloc que vous avez eu en échantillon.
ndvi(X,2,3)Si tout se déroule comme prévu, vous pouvez désormais demander à RasterMath d'ajouter la fonction qui calcule le NDVI. Pour ajouter des arguments qu'utilise votre fonction, il suffit de les renseigner par leur nom (ici : r=2,ir=3).
my_image.add_function(ndvi,r=2,ir=3,out_image='/tmp/ndvi.tif')Il est à noter que l'on peut donner autant de fonctions que l'on souhaite à RasterMath.
Une fois que vos fonctions sont validées, il ne vous reste plus qu'à demander l'exécution du calcul.
my_image.run()Maintenant lire et écrire sur une raster ne prend plus que 3 lignes, il vous suffit juste d'avoir une fonction qui traite un tableau 2d ou 3d, et le tour est joué !
Scikit-Learn est une bibliothèque python d'apprentissage automatique. Elle regroupe de nombreux pré-traitements et algorithmes.
Pour installer et/ou mettre à jour scikit-learn (et joblib qui gère les calculs sur plusieurs cœurs), il vous suffit d'ajouter l'option -U à la fin du code :
python3 -m pip install scikit-learn joblib --user -UPour apprendre un modèle il faut deux éléments :
- X : un tableau qui contient pour chaque échantillon (un échantillon par ligne) les variables en colonnes
- y : le label de chaque échantillon (un label par ligne).
from sklearn.datasets import load_iris
X,y = load_iris(return_X_y=True)Une fois qu'on a notre jeu de données, il faut choisir un algorithme.
Scikit-learn propose de très nombreux algorithmes. Je vous en ai sélectionné 3 :
Prenons le cas de l'arbre de décision (Decision Tree) :
from sklearn.tree import DecisionTreeClassifier
model = DecisionTreeClassifier()
model.fit(X,y) # on apprend le modèle avec nos échantillons
model.predict(X) # on prédit le modèle sur nos échantillonsDans le cas précédent, nous avons entraîné un modèle avec l'ensemble de données à notre disposition, nous ne pouvons donc pas le valider.
Pour valider le modèle, nous devons mettre de côté quelques échantillons afin de voir si l'algorithme peut retrouver le label sans jamais les avoir vu.
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)En utilisant la fonction train_test_split, nous gardons 33% des données pour valider notre modèle et ainsi estimer sa qualité, et 66% des données pour l’entraîner.
from sklearn.tree import DecisionTreeClassifier
model = DecisionTreeClassifier()
model.fit(X_train,y_train) # on apprend le modèle avec nos échantillons
y_pred = model.predict(X_test) # on prédit le modèle sur nos échantillonsLe plus y_pred sera semblable au y_test, le plus performant sera le modèle.
Tout d'abord on peut calculer la matrice de confusion pour regarder la qualité générale du modèle et voir où le modèle s'est trompé et où il a bien réussi.
from sklearn.metrics import confusion_matrix,accuracy_score,f1_score
confusion_matrix = confusion_matrix(y,y_pred)
Si vous avez oublié comment lire une matrice de confusion, je vous conseille de lire la page wikipédia qui y est dédiée.
Accord global (ou Overall accuracy), il s'agit du total de la diagonale sur le total d'échantillons prédits:
oa = np.sum(np.diag(confusion_matrix))/np.sum(confusion_matrix)On peut aussi calculer l'OA avec scikit-learn
oa = accuracy_score(y,y_pred)Le F1 est le nombre de commissions et d'omissions par classe :
f1 = f1_score(y,y_pred)- Quelle est la qualité globale de votre modèle ? Calculer le f1 moyen et l'accord global
- Quelle classe a été la mieux prédite ?
- Quelle classe a été la moins bien prédite ?
Pour extraire les données d'une image et le label d'un polygone, nous utilisons la fonction extract_ROI de Museo ToolBox.
in_field correspond à la colonne qui contient pour chaque point/polygone le label de la classe.
X,y = mtb.processing.extract_ROI('sentinel2_3a_20180815.tif','ROI.gpkg','class')
Nous pouvons demander autant de colonnes de type numérique que l'on souhaite à extract_ROI. Par exemple pour obtenir le groupe, il suffit de rajouter un argument :
X,y,g = mtb.processing.extract_ROI('sentinel2_3a_20180815.tif','ROI.gpkg','class','group')
Ensuite, vous pouvez appliquer la même méthode d'apprentissage automatique que vous avez vu précédemment, ou alors partir sur une validation croisée pour fixer les paramètres de l'algorithme.
- Apprenez un modèle avec les données (X, y et g) extraite de l'image Sentinel-2
- Utiliser une grille de recherche pour apprendre à partir d'une validation croisée (de type Leave-One-SubGroup-Out)
- Fixer les hyper-paramètres de votre algorithme. Pour RandomForestClassifier, vous pouvez faire varier le n_estimators (nombre d'arbres), ou le min_sample_split par exemple.
- Quelles sont les paramètres qui permettent d'avoir le meilleur résultat ? (nombre d'arbres et min_sample_split)
- Prédire maintenant le modèle sur l'image Sentinel-2 en utilisant RasterMath.