(Un)fold blocks

Objectifs d'apprentissage

À la fin de ce TD, tu seras en mesure de classifier de manière "experte" une image de télédétection multispectrale à partir d'une analyse de la dynamique spectrale de l'image. De façon plus détaillée, tu seras en mesure de :

• identifier une classe spectrale :
  • à l'aide d'une histogramme d'une bande ;
  • à l'aide d'une histogramme 2D ;
• caractériser une classe spectrale par :
  • des seuil sur un histogramme 1D ;
  • des droites sur un histogramme 2D ;
• produire une carte de classes spectrales préalablement identifiées
  • en reclassant des plages de valeurs continues en labels discrets ;
  • à l'aide d'outils type calculatrice raster ;
• caractériser le comportement spectrale d'un occupation du sol à l'aide d'une histogramme 2D.

Pour atteindres ces objectifs d'apprentissages, tu seras amené à utiliser les outils suivants :

• Histogramme ;
• RasterDataPlotting ;
• Calculatrice raster et/ou BandMath ;
• RadiometricIndices.

Indice de végétation

Parmi les indices radiométriques existants, seul le NDVI sera utilisé dans ce TD car il est l'indice de végétation le plus utilisé en télédétection.

Pour rappel, la formule du NDVI est la suivante :

$$NDVI = \frac{IR - R}{IR + R}$$

Il est possible de calculer le NDVI de plusieurs manières. Soit :

• avec la calculatrice raster ;
• avec l'application OTB BandMath, dans QGIS, avec Mapla ou en ligne de commande ;
• avec l'application OTB RadiometricIndices, dans QGIS, avec Mapla ou en ligne de commande.

Les deux premiers outils permettent de faire des opérations entre les bandes des différentes images que tu as chargées dans QGIS. On peut les utiliser pour calculer le NDVI mais également pour d'autres calculs. En revanche, les deux outils suivent une syntaxe qui n'est pas la même (ce serait trop facile). Pour la calculatrice raster tu peux te référer à cette doc et pour la syntaxe de BandMath à cette doc.

La Calculatrice Raster est disponible depuis le menu Raster ► Calculatrice Raster.

Pour la calculatrice raster il faudra rentrer la formule suivante pour le calcul du NDVI (dans la partie expression de la calculatrice raster sur la figure ci-dessus):

("fabas_octobre_2013@4"-"fabas_octobre_2013@1")/("fabas_octobre_2013@4"+"fabas_octobre_2013@1")

L'application BandMath est quant à elle disponible depuis la Boîte à outils de traitements ► OTB ► Image Manipulation ► BandMath.

L'expression à rentrer pour le calcul du NDVI (dans la partie expression sur la figure ci-dessus) est :

(im1b4-im1b1)/(im1b4+im1b1)

Voici également un exemple de syntaxe de pour l'application OTB BandMath en ligne de commande :

otbcli_BandMath \
    -il fabas_octobre_2013.tif \
    -out ndvi_fabas.tif \
    -exp "(im1b4-im1b1)/(im1b4+im1b1)"

Enfin, voici un exemple de syntaxe de pour l'application OTB RadiometricIndices en ligne de commande :

otbcli_RadiometricIndices \
    -in fabas_12_10_2013.tif \
    -list Vegetation:NDVI \
    -channels.red 1 \
    -channels.nir 4  \
    -out ndvi_fabas.tif

Tu vas te tromper, et c'est normal ;)

Tu vas faire des erreurs, c'est normal. Ce serait de ne pas en faire qui n'est pas normal. Certaines peuvent être dues à une mauvaise utilisation des paramètres d'une application, il faut dans ce cas là se référer à l'aide des applications. D'autres erreurs sont dues à de l'inattention, parce que tu es allé·e trop vite, emporté·e par ton enthousiame.

Dans tous les cas, prends le temps de regarder les messages d'erreur. Ils ne sont pas toujours faciles à comprendre mais ils t'aideront souvent (pas toujours malheureusement) à trouver la source de ton erreur (voir exemple sur la figure ci-dessous).

Seuillage d'image (classification experte)

Dans ce TD, tu vas réaliser ta première classification : tu vas attribuer à chaque pixel d'une image un label qui correspondra à une occupation du sol. Aujourd'hui, nous nous focaliserons sur une approche non supervisée et manuelle : tu vas déterminer toi même les règles permettant de discrimer les labels de chaque pixel.

En fonction de ton interlocuteur, le terme utilisé pourra varier. Certains parleront de segmentation (simple) d'images, d'autres de classification experte (au sens d'expertisée par une personne) ou encore de seuillage.

De quelle occupation du sol parle-t-on ?

Dans ce TD nous ne connaissons pas a priori les occupations du sol qu'il est possible de discriminer (--> c'est pour ça qu'on parle d'approche non supervisée). Tu vas analyser les dynamiques spectrales des images, aussi bien des bandes natives (Rouge, Vert, Bleu et Proche Infrarouge) que de l'indice NDVI. Tu verras dans un premier temps comment analyser la dynamique spectrale d'une bande (ou d'un indice) puis dans un second temps comment analyser la dynamique de deux bandes en même temps.

Analyse spectrale 1D

Dans le meilleur des mondes, les occupations du sol principales présentes dans une image peuvent être identifées sur l'histogramme d'une bande. Dans l'histogramme théorique ci dessous, les classes de sol nu, Herbe, Ligneux bas et Ligneux Hauts se distinguent par des modes différents (ou des pics). Dans ce cas là on peut imaginer différencier ces classes grâce à des seuils entre les modes. La "hauteur" des modes nous renseigne quant à elle à quel point c'est classe est présente dans la zone d'étude.

Aujourd'hui, les classes d'occupation du sol que tu vas classifier ne correspondront pas forcément à du sol nu, de l'herbe, des ligneux bas ou des ligneux hauts. À toi de le découvrir !

Tu peux visualiser les histogrammes d'une bande dans QGIS en : Clique Droit sur la couche votre choix ► Propriétés ► Histogramme :

Remarque : les histogrammes peuvent parfois être plus difficiles à lire selon les machines avec lesquelles tu vas travailler.

Seuillage 1D

Ici, l'objectif est de regrouper les pixels ayant des caractéristiques spectrales similaires. L'analyse de l'histogramme que tu as faite dans la partie précédente t'a permis d'identifier plusieurs classes spectrales. Tu vas maintenant attribuer un label à chacun des pixels grâce à l'aide de l'application BandMath ou de la calculatrice raster.

Le principe est d'attribuer un label (par exemple 2) à des pixels ayant une gamme de valeurs que tu auras déterminée (par exemple les pixels compris entre -0.2 et 0). Selon l'application choisie la syntaxe n'est pas la même.

La syntaxe entre les calculatrices est légèrement différente mais le principe reste le même.

Syntaxe et fonctionnement des calculatrices raster

Pour les explications, prenons l'exemple d'un raster ImgA qui contient des entiers entere 1 et 7. Pour sélectionner les pixels compris entre 4 et 5 la syntaxe à utiliser est la suivante pour la calculatrice raster est :

("ImgA@1" > 3 ) AND ("ImgA@1" < 6)

ou la suivante pour l'application BandMath :

(im1b1 > 3 ) and (im1b1 < 6)

Si cette formule est entrée telle quelle, le raster en sortie aura pour valeur 1 là où la condition est respectée et 0 ailleurs. Pour attribuer la valeur de 3 à la plage de valeur sélectionnée, on multiplie tout simplement le résultat par 3, les autres zones auront toujours pour valeur 0. Ces deux étapes sont résumées dans la figure ci dessous.

Il est possible d'obtenir l'image C avec un seul traitement en entrant une formule légèrement plus complexe :

Si tu veux attribuer un autre label (par ex. 2) à une autre plage de valeurs (par ex. les valeurs inférieures à 3), tu peux additionner des conditions ensemble :

Erreur courante : conditions qui ne sont pas mutuellement exlusives

La combinaison de plusieurs conditions n'est possible que si les plages de valeurs sélectionnées sont mutuellement exclusives. C'était le cas dans les exemples ci dessus. Imaginons maintenant que la condition pour obtenir l'image D soit légèrement différente et est "ImgA@1" < 5 au lieu de "ImgA@1" < 4 (voir figure ci-dessous).

Dans ce cas l'image E en sortie comportera 4 valeurs au lieu de 3 :

• 0 pour les pixels qui ne respectent aucune condition ;
• 3 pour les pixels qui respectent uniquement la première condition ;
• 2 pour les pixels qui respectent uniquement la deuxième condition ;
• 5 (2 + 3) pour les pixels qui respectent les deux conditions --> résultat non souhaité (dans la majorité des cas).

Si tu souhaites simplement masquer une partie de l'image, c'est possible aussi. Dans ce cas là, selectionne la plage de valeurs à garder et multiplie la par l'image d'origine :

(("fabas_decembre_2013_ndvi@1" > -0.1) AND ("fabas_decembre_2013_ndvi@1" < 0.1)) * 2 +
((("fabas_decembre_2013_ndvi@1" > 0.2) AND ("fabas_decembre_2013_ndvi@1" < 0.35))) * 5

("fabas_decembre_2013_ndvi@1" < -0.1) +
(("fabas_decembre_2013_ndvi@1" > -0.1) AND ("fabas_decembre_2013_ndvi@1" < 0.1)) * 2 +
((("fabas_decembre_2013_ndvi@1" > 0.2) AND ("fabas_decembre_2013_ndvi@1" < 0.35))) * 5 +
("fabas_decembre_2013_ndvi@1" > 0.35) * 7

(("fabas_decembre_2013_ndvi@1" > -0.1) AND ("fabas_decembre_2013_ndvi@1" < 0.1)) * 2 +
((("fabas_decembre_2013_ndvi@1" > 0.2) AND ("fabas_decembre_2013_ndvi@1" < 0.35))) * 5 +
((("fabas_decembre_2013_ndvi@1" > 0.05) AND ("fabas_decembre_2013_ndvi@1" < 0.25))) * 3

Comme tu as pu le voir dans ces différents exemples, les syntaxes de BandMath et de la calculatrice raster de QGIS ne sont pas très différentes. Les principales différences sont résumées dans le tableau ci-dessous :

Rappel : l'avantage de BandMath par rapport à la calculatrice raster est que tu peux définir le type de l'image en sortie (par exemple en entier 8 bits).

Mise en pratique

Analyse spectrale 2D

Pour cette partie, il faut installer le plugin RasterDataPlotting (installer un plugin).

L'objectif est de repérer des classes spectrales à partir des bandes Rouge et Infrarouge. On parle de Scatter Plot 2D. Cette représentation est très utile pour distinguer les différents types de végétation, de sol et d’enneigement (voir figure ci-dessus).

Tu peux obtenir cette représentation en configurant le pluggin RasterDataPlotting comme suit :

Comment faire pour faire le lien entre le scatter plot et les pixels de l'image ?

Deux fonctionnalités peuvent t'aider à faire le lien entre les classes spectrales que tu as distinguées sur le scatter plot et les occupations du sol. La première est le fait que seul les pixels contenus dans l'emprise de ton Canevas QGIS sont affichés dans dans le scatter plot : pour voir le comportement spectrale d'une occupation du sol, zoome sur celle-ci comme sur l'exemple ci-dessous.

De manière inverse tu peux mettre en évidence sur le Canevas les pixels correspondant à une classe spectrale en utilisant le bouton Spectral ROI comme sur l'exemple ci-dessous. Attention toutefois, l'utilisation des Spectral ROI demande un peu de doigté et sera certainement source de frustration au bout de quelques essais infructueux. Prends patience ! cet outil en vaut la peine.

Astuce : équation d'une droite

À partir de deux points (identifiés visuellement, en zoomant sur le scatter plot) vous pouvez déterminer l'équation d'une droite. Ça remonte à loin tout ça, je sais je sais ...

Trop loin ?

Seuillage du Scatter splot

Ici, l'objectif est encore et toujours de regrouper les pixels ayant des caractéristiques spectrales similaires. L'analyse du scatter plot que tu as faite dans la partie précédente t'a permis d'identifier plusieurs classes spectrales. Tu vas maintenant attribuer un label à chacun des pixels.

La distinction de deux classes (2 et 0 par exemple) séparées par une droite dessinée dans le plan spectral constitué par la bande Infra rouge (bande 4) et la bande Rouge (bande 1) comme sur la figure ci-dessus, peut se faire avec BandMath comme suit :

(im1b4 > (a * im1b1 + b)) * 2

: Seuillage d'une image en utilisant deux bandes.

1. À l'aide de l'application BandMath ou la calculatrice raster, seuille ton image et enregistre ton seuillage finale (par ex: fabas_octobre_2013_seuil.tif). Ne supprime pas l'image à la fin du TD, on s'en reservira plus tard (quel suspense ! 😱)
2. Quel type de d'image en sortie faut il renseigner ? et surtout, pourquoi ?
3. Change éventuellement le mode de représentation de l'image : Clique droit ► propriétés ► syle de couche ► Palette /Uniques (Doc) pour que le résultat de ton seuillage soit plus lisible :

4. Les classes issues ton seuillage sont-elles conformes à ce que tu avais repéré précédemment ?

5. Marque dans l'exercice comme terminé dans Moodle.