Auteur : Paul Passy

Licence :

Outils de télédétection

Dans cette partie, nous verrons quelques outils utiles à appliquer (si besoin) avant une étude de télédétection. Nous détaillerons par exemple la construction de compositions colorées ou le principe du pansharpening.

Table des matières

Astuce

Concernant les corrections radiométriques et atmosphériques, il est possible de télécharger facilement les images déjà fournies en réflectance de surface (niveau L2) sur les différents poratils de téléchargement (Télécharger des images satellites). Il ne faut pas hésiter à récupérer ces données pour nous éviter ces corrections.

Les compositions colorées en télédétection 

Dans cette section nous allons découvrir les principes théoriques de ce que nous appelons compositions colorées en télédétection. Le principe peut être appliquée à des images optiques ou radars. Les compositions colorées sont une manipulation basiques dans un projet faisant intervenir un traitement d’images satellites. Cette manipulation permet d’avoir un premier aperçu de sa zone d’étude et de ses caractéristiques.

La plupart des logiciels de géomatique permettent de construire concrètement ce genre de compositions. Pour des exemples pratiques, reportez vous à la section dédiée.

Les fondements : la synthèse additive des couleurs 

En ce qui concerne la synthèse des couleurs, il faut différencier la synthèse additive de la synthèse soustractive. La seconde est celle que nous découvrons en premier dans les petites classes. Elle repose sur trois couleurs primaires à base de pigments que nous mélangeons afin d’obtenir toutes les couleurs possibles. Ce n’est pas celle-ci que nous utilisons en télédétection. Nous utilisons la synthèse additive, basée sur les propriétés électromagnétiques (lumineuses) des couleurs. Cette synthèse repose également sur trois couleurs primaires, aussi appelées couleurs fondamentales.

Ces trois couleurs sont le rouge, le vert et le bleu. C’est pourquoi le terme de RVB ou RGB (pour Red, Green, Blue) est couramment utilisé. À partir du mélange des ces trois couleurs, en jouant sur les proportions, il est possible de retrouver toutes les couleurs existantes. Le mélange du bleu et du rouge donne du magenta, celui du bleu et du vert donne du cyan et celui du vert et du rouge donne du jaune. Enfin le mélange du bleu, du rouge et du vert en même temps donne du blanc. Alors qu’une absence totale de ces trois couleurs donne du noir (Fig. 311). Au final, en jouant sur les proportions des mélanges, il est possible d’obtenir toute la gamme de couleurs. Cette synthèse additive est également employée par tous les écrans.

alternate text — Fig. 311 Le principe de la synthèse additive des couleurs.

Note

Vous avez sûrement déjà entendu un ami parfois agaçant vous dire « le blanc et le noir ce ne sont pas des couleurs car c’est un mélange de tout (dans le premier cas) et une absence de tout (dans le second cas) ! » Certes, mais à ce compte, le jaune n’est pas non plus une couleur puisqu’il résulte lui aussi d’un mélange. Au final, seul le bleu, le vert et le rouge seraient des couleurs ? La définition d’une couleur est finalement plus compliquée qu’attendue et nous ne nous y attarderons pas.

Application à la télédétection optique 

En télédétection optique, ce principe est appliqué à la combinaison de bandes spectrales. Chaque bande spectrale est codée sous forme d’une image raster composée de pixels. À chaque pixel est associée une mesure de luminance qui équivaut à l’intensité de l’énergie réfléchie par le pixel sous-jacent dans la bande spectrale donnée. Plus cette intensité est faible, plus la valeur du pixel tend vers 0. Plus l’intensité est élevée, plus la valeur du pixel est élevée. La borne supérieure dépend de l’encodage de l’image. Une fois les luminances converties en réflectance, ces valeurs s’échelonnent entre 0 et 1. Un pixel qui ne réfléchit rien aura une valeur de 0 et un pixel qui réfléchit tout aura une valeur de 1.

D’un point de vue visuel, ces valeurs sont traduites par des teintes de gris. Un pixel qui ne réfléchit rien, donc qui tend vers 0, apparaîtra noir. À l’autre extrémité, un pixel qui réfléchit tout, donc qui tend vers 1, apparaîtra blanc. Pour les valeurs intermédiaires, nous retrouverons toutes les gammes de gris imaginables. C’est pourquoi, lorsque nous chargeons une bande spectrale isolée dans un logiciel SIG, le raster apparaît en teintes de gris (Fig. 312).

Sur la figure (Fig. 312), nous constatons que les pixels correspondant aux bois de Vincennes et de Boulogne apparaissent très clairs. Nous en déduisons que ces pixels réfléchissent beaucoup dans la bande 5 de OLI qui correspond au proche infrarouge. Au contraire, les surfaces bâties du centre de l’agglomération apparaissent très sombres. Les pixels de ces zones réfléchissent donc peu dans ce domaine du proche infrarouge.

L’idée est de combiner trois bandes spectrales en les transformant chacune dans une teinte de rouge, de vert et de bleu. Par exemple, la figure (Fig. 313) présente la même bande spectrale que sur la figure (Fig. 312) mais codée en teintes de rouge plutôt qu’en teintes de gris. Le point délicat est de bien comprendre que les zones qui apparaissent très claires sur la figure (Fig. 312) correspondent à celles qui apparaissent très rouges sur la figure (Fig. 313). Au contraire, les zones très sombres de la figure (Fig. 312) apparaissent très claires sur la figure (Fig. 313). L’explication est simple et logique. Un pixel avec une forte réflectance a une valeur élevée, il va donc être codé avec beaucoup de blanc dans le cas de la figure (Fig. 312) mais avec beaucoup de rouge dans le cas de la figure (Fig. 313). Au contraire, un pixel qui réfléchit peu va être codé avec peu de blanc, donc presque noir, dans le cas de la figure (Fig. 312) et peu de rouge, donc presque blanc, dans le cas de la figure (Fig. 313).

Les trois bandes spectrales choisies pour la composition colorée doivent être « colorées » de la même manière, la première en bleu, la deuxième en vert et la troisième en rouge. Par convention, l’ordre des couleurs et des bandes doit être respectée. La plus « petite » bande sera codée en bleu et la plus « grande » en rouge. Une fois ces colorations effectuées, les trois bandes « colorées » sont « mélangées » afin d’obtenir une vue de la même scène mais en couleurs (figure Fig. 314). Notez bien que la « coloration » en bleu, vert et rouge des bandes ainsi que le « mélange » se fait de façon automatique par le logiciel.

Interprétation d’une composition colorée 

L’interprétation des compositions colorées n’est pas compliquée mais demande d’acquérir une certaine gymnastique. Nous allons nous entraîner en interprétant la composition résultante de la figure Fig. 314, présentée en plus grand format sur la figure Fig. 315.

Sur cette composition colorée, nous voyons tout d’abord des zones qui apparaissent en rouge plus ou moins foncé. Ces zones correspondent donc à des pixels avec beaucoup de rouge et très peu de vert et de bleu. Comme nous avons codé la bande 5 en rouge, la bande 3 en vert et la bande 2 en bleu, ça signifie que ces zones réfléchissent beaucoup dans la bande 5 et très peu dans les autres. Or la bande 5 correspond au proche infrarouge du capteur OLI, la bande trois au vert et la bande 2 au bleu. Ainsi, les zones rouges réfléchissent préférentiellement dans le proche infrarouge et peu dans le vert et le bleu. En s’appuyant sur les signatures spectrales connues, nous pouvons constater que ce comportement est typique des zones en végétation. En effet, nous reconnaissons bien en rouge les bois de Vincennes, de Boulogne ainsi que les forêts de grande couronne francilienne.

En suivant la même logique, la Seine et ses affluents, ainsi que les lacs apparaissent en bleu sombre. Ces zones réfléchissent en effet peu dans les trois bandes spectrales choisies mais tout de même un peu plus dans la bande 2 qui correspond au domaine du bleu. Ce comportement est également en cohérence avec la signature spectrale de l’eau.

Les zones bâties apparaissent plus ou moins cyan. Elles correspondent à un mélange de bleu et de vert mais pas de rouge. Ce sont donc des zones qui réfléchissent préférentiellement dans les domaines du bleu et du vert par rapport à celui du proche infrarouge. Nous retrouvons bien la signature spectrale des surfaces bâties.

Enfin, en zoomant, nous constatons qu’une portion du jardin des Tuileries apparaît très claire, presque blanche. Cette clarté signifie que cette zone réfléchit intensément dans les trois bandes choisies, à savoir le proche infrarouge, le bleu et le vert. Selon les signatures spectrales, ce comportement est typique soit de la neige, ce qui n’est bien sûr pas le cas ici, soit d’un sol nu très clair. Le revêtement des allées de ce jardin est en effet de couleur très clair.

Astuce

Il est tout à fait possible, et parfois très intéressant, de faire des compositions colorées multi-dates. Par exemple, imaginons qu’avec une série d’images Landsat 8 nous calculions un NDVI pour un mois hivernal, un mois printanier et un mois estival. Il est ensuite possible de colorer le NDVI hivernal en rouge, le NDVI printanier en vert et le NDVI estival en bleu. Les zones bleu de la composition résultante correspondra alors aux zones végétalisées en été mais non végétalisées en hiver et au printemps.

Composition colorée avec QGIS 

Version de QGIS : 3.26.1

Nous allons voir dans cette partie comment réaliser une composition colorée avec QGIS. Cela se fait très simplement à partir d’un raster multi-bandes. La construction d’un raster multi-bandes est décrite dans la partie dédiée : Raster multi-bandes. La plupart du temps nous réalisons une composition colorée selon différentes bandes spectrales d’une même image satellite. Mais rien n’empêche de réaliser une composition colorée à partir d’autres rasters comme des NDVI calculés sur une même zone mais à des dates différentes.

Astuce

Il sera plus facile de créer des compositions colorées si le raster multi-bandes a été construit en nommant explicitement les différentes bandes comme décrit ici : Raster multi-bandes en dur (avec noms).

Ici nous nous contenterons de réaliser une composition colorée en vraies couleurs et une en fausses couleurs mettant en avant la végétation à partir d’une image Sentinel-2 de la région parisienne prise le 24 juin 2020. Une fois le raster multi-bandes (qui peut être un raster virtuel) chargé dans QGIS, nous faisons un clic droit sur cette couche dans le panneau des couches puis nous sélectionnons Propriétés ‣ Symbologie. La fenêtre de réglage de la symbologie apparaît (Fig. 316).

Nous commençons par construire une composition colorée en vraies couleurs. Nous colorons donc la bande du domaine du Bleu (bande 2) en bleu, celle du domaine du Vert (bande 3) en vert et celle du domaine du Rouge (bande 4) en rouge. Si besoin, reportez vous à la section présentant les bandes de Sentinel-2 pour revoir les équivalences (Les images Sentinel-2). Ainsi, pour la Bande rouge nous choisissons la Bande 04, pour la Bande verte la Bande 03 et pour la Bande bleue la Bande 02. Il n’y a plus qu’à cliquer sur OK et la composition colorée apparaît (Fig. 317 A).

Avertissement

Il faut être très prudent lorsque nous choisissons nos bandes si les bandes n’ont pas été nommées. Le numéro des bandes qui apparaît dans ce menu n’est que le numéro d’ordre de la bande dans le raster multi-bande, en aucun cas le numéro de bande de l’image de base. Si votre raster multi-bandes contient les bandes 2, 3 et 4 de votre image Sentinel-2, alors la bande 2 de Sentinel sera vue comme étant la bande 1 du raster multi-bandes et ainsi de suite.

Pour réaliser une composition colorée en fausses couleurs mettant en avant la végétation nous procédons de même mais mettons cette fois-ci pour la Bande rouge la bande du proche infrarouge soit la bande 08. Les autres bandes restent inchangées : Bande bleue avec la bande 02 et Bande verte avec la bande 03 (Fig. 317 B).

Astuce

Les réglages des contrastes peuvent se faire au niveau du menu de symbologie de la couche mais le plus simple est d’utiliser à posteriori l’outil d’amélioration des contrastes disponible dans la barre d’outils raster sous l’icône Histogramme cumulatif de l'emprise locale .

Composition colorée avec R 

Version de R : 4.2.1

Il est possible de créer des compositions colorées dans R. Il est nécessaire de charger au préalable un raster multi-bandes ou bien d’en créer un (Création d’un raster multi-bandes avec R). Dans l’exemple suivant nous travaillons sur des objets SpatRaster chargés via le package terra. Nous commençons par charger un raster geoTIFF qui est déjà multi-bandes, composé des quatre premières bandes spectrales d’une image Landsat 8 OLI prise au-dessus des Bouches-du-Rhône. La composition colorée se fait alors via la fonction plotRGB. Le tout est résumé dans les lignes suivantes.

# chargement de la librairie terra
library(terra)

# chargement du raster multi-bandes
stack_Landsat <- terra::rast("stack.tif")

# création de la compostion colorée fausses couleurs
plotRGB(stack_Landsat, r=4, g=2, b=1, stretch='hist')

Dans la fonction plotRGB, nous indiquons les bandes à colorer en Rouge, Vert et Bleu via les clés r=, g= et b= suivies des numéros de bandes à colorer. Ces numéros de bandes correspondent simplement à l’ordre d’empilement dans le raster multi-bandes. Nous indiquons également stretch='hist' afin d’améliorer les contrastes. Plus d’options sont possibles comme indiqué sur la page terra dédiée. Au final, nous obtenons la figure suivante (Fig. 318).

Correction radiométrique 

Lorsque nous récupérons une image satellite dont les pixels représentent des niveaux d’énergie, souvent des images en niveau L1, il est nécessaire d’y appliquer au moins une correction radiométrique afin d’obtenir une image en Réflectance en haut de l’atmosphère (Reflectance Top Of Atmosphere (TOA)). Cette correction est nécessaire pour toutes les images Landsat en niveau L1, se reporter à la section consacrée aux niveaux de traitements pour les détails : Niveaux de traitements. Cette correction peut se faire avec différents outils dont certains sont présentés ici.

Nous appliquerons cette correction sur une image Landsat OLI du 2 avril 2021 prise au-dessus de la commune de Mbandjock, située dans la Haute-Sanaga au centre du Cameroun.

Correction radiométrique dans QGIS avec SCP 

Version de QGIS : 3.20.1

Version de SCP : 7.9.5

Le module SCP de QGIS permet d’appliquer cette correction facilement. Le plus sûr est d’appliquer cette correction sur toutes les bandes et l’intégralité de la scène, puis de supprimer à posteriori les bandes corrigées qui ne nous intéressent pas et de découper ou masquer celles que nous garderons.

Une fois QGIS lancé, nous allons dans le menu SCP ‣ Pré-traitement ‣ Landsat. La fenêtre suivante apparaît (Fig. 319).

Sur la ligne Dossier contenant les bandes Landsat nous pointons vers le répertoire qui contient les bandes Landsat. Ce répertoire doit être exactement celui qui a été téléchargé avec les mêmes fichiers et la même structure. Sur la ligne Sélectionner le fichier MTL nous pointons vers le fichier texte de métadonnées qui se trouve dans le répertoire des images. Ce fichier texte se termine par _MTL.txt. C’est dans ce fichier que se trouve les paramètres nécessaires à la correction radiométrique comme la distance entre la Terre et le Soleil, l’élévation du Soleil par rapport à l’horizon ou l’angle de vue du satellite.

Comme la structure du fichier est normalisée, SCP comprend de lui même quel raster correspond à quelle bande et quelles corrections administrées. Ces informations sont récapitulées dans le tableau sous la ligne Métadonnée. Il suffit ensuite de lancer le processus en cliquant sur Lancer. Il nous est alors demandé de sélectionner un répertoire dans lequel les bandes spectrales corrigées apparaîtront, les rasters résultants s’affichent automatiquement dans QGIS. Les rasters corrigés portent le même nom que les rasters initiaux mais sont précédés du préfixe RT_.

Note

La bande panchromatique (Bande 8 de Landsat 8) ainsi que la bande 9 correspondant aux cirrus ne sont pas corrigées.

Après correction, nous constatons que les pixels ont bien des valeurs comprises entre 0 et 1. Cette valeur est maintenant une réflectance en haut de l’atmosphère (TOA).

Note

Il est rare de travailler sur une scène Landsat entière et sur toutes les bandes spectrales. L’utilisateur doit maintenant supprimer manuellement les bandes qui ne serviront pas et découper, masquer ou mosaïquer celles d’intérêt, voir la section dédiée aux outils rasters : Outils raster.

Correction atmosphérique 

Lorsque nous disposons d’une image en réflectance en haut de l’atmosphère (TOA), il peut être utile de la convertir en réflectance de surface. Cette étape est notamment indispensable dans le cas d’un suivi dans le temps d’une même région. En effet, grâce à cette correction nous serons sûrs que les changements observés dans le temps représentent des changements de surface et non pas de conditions atmosphériques. Cette correction atmosphérique peut se faire via différents outils.

La correction atmosphérique peut être réalisée selon un modèle physique ou empirique. Seule l’utilisation du modèle empirique sera présentée ici. Ce modèle est nommé DOS pour Dark Object Subtraction. Il peut être décliné en plusieurs versions notées DOS1 à DOS6.

Correction atmosphérique dans QGIS avec SCP 

Version de QGIS : 3.20.1

Version de SCP : 7.9.5

Image Landsat 

Nous présenterons ici le cas d’une correction atmosphérique appliquée à une image Landsat. Nous appliquerons plus précisément cette correction atmosphérique à une image Landsat OLI du 2 avril 2021 prise au-dessus de la commune de Mbandjock, située dans la Haute-Sanaga au centre du Cameroun.

Le menu employé est le même que celui dédié à la correction radiométrique, ce qui permet d’enchaîner les deux corrections en une seule manipulation. Une fois QGIS lancé, nous allons dans le menu SCP ‣ Pré-traitement ‣ Landsat. La fenêtre suivante apparaît Fig. 320.

Avertissement

Il semble qu’avec le module SCP, il ne soit pas possible d’effectuer une correction atmosphérique seule. Elle doit obligatoirement se faire dans la foulée d’une correction radiométrique (Correction radiométrique). Ainsi, si par hasard vous récupérez une image Landsat en réflectance TOA, il ne sera pas possible d’y appliquer seulement une correction atmosphérique en utilisant SCP.

Sur la ligne Dossier contenant les bandes Landsat nous pointons vers le répertoire qui contient les bandes Landsat. Ce répertoire doit être exactement celui qui a été téléchargé avec les mêmes fichiers et la même structure. Sur la ligne Sélectionner le fichier MTL nous pointons vers le fichier texte de métadonnées qui se trouve dans le répertoire des images. Ce fichier texte se termine par _MTL.txt. Ce fichier de métadonnées n’est pas utile à la correction atmospéhrique, mais SCP est obligé de partir d’une image en niveau d’énergie pour obtenir une image en réflectance de surface en passant implicitement par une correction radiométrique. Pour indiquer à SCP que nous souhaitons une correction atmosphérique, il ne faut pas oublier de cocher la case Appliquer correction atmosphérique DOS1. Nous cliquons ensuite sur Lancer. Il nous est demandé de renseigner le répertoire qui contiendra les images corrigées.

Les bandes spectrales corrigées portent le même nom que les bandes initiales mais précédées du préfixe RT_. Nous constatons que les pixels ont bien maintenant des valeurs comprises entre 0 et 1 qui correspondent à des réflectances de surface. La bande panchromatique et celle correspondant au cirrus ne sont pas corrigées.

Note

En comparant, pour une bande donnée, les valeurs des pixels en réflectance en haut de l’atmosphère (TOA) et en réflectance de surface, nous constatons logiquement que les réflectances de surface sont moins fortes. En effet, la réflectance atmosphérique a été enlevée.

Image Sentinel-2 

Pour transformer une image Sentinel-2 de niveau L1C (Niveaux de traitements Sentinel-2) en réflectance de surface depuis une image en réflectance au sommet de l’atmosphère avec SCP, la procédure est très semblable. Nous allons dans le menu SCP ‣ Pré-traitement ‣ Sentinel-2. La fenêtre suivante apparaît Fig. 321.

À la ligne Dossier contenant les bandes Sentinel-2, nous pointons vers le répertoire qui contient les rasters des bandes spectrales. Généralement ces rasters se trouvent dans le sous répertoire GRANULE du répertoire se terminant par .SAFE. Dans la ligne Sélectionner un fichier de métadonnées, nous pointons vers le fichier texte de métadonnées au format .XML du répertoire parent. Mais nous pouvons également laisser cette ligne vide, SCP trouvera le fichier par lui-même. Il ne faut bien sûr pas oublier de cocher la case Appliquer correction atmosphérique DOS1. Nous cliquons ensuite sur Lancer. Il nous est demandé de renseigner le répertoire qui contiendra les images corrigées.

Les bandes spectrales corrigées portent le même nom que les bandes initiales mais précédées du préfixe RT_. Nous constatons que les pixels ont bien maintenant des valeurs comprises entre 0 et 1 qui correspondent à des réflectances de surface. Par défaut, les bandes 1, 9 et 10 à 60 m de résolution spatiale ne sont pas corrigées (Les images Sentinel-2). Il est possible de les ajouter en cochant la case Preprocess bands 1, 9, 10 du menu précédent (Fig. 321).

Masquer les nuages 

Il est rare d’avoir une image sans aucun nuage pour une date donnée. Il est fréquent de devoir tolérer un certain pourcentage de nuages sur une scène que nous souhaitons traiter. Par contre, ces nuages ne doivent pas rentrer dans l’analyse de l’image que nous ferons par la suite. En effet, calculer un NDVI sur un nuage n’aura, par exemple, aucun sens. De même, un nuage ne peut pas être une classe d’occupation du sol. Il est donc nécessaire de créer un raster sur lequel ne figureront que les nuages (et éventuellement d’autres objets indésirables comme leurs ombres, les cirrus ou les aérosols). Ce raster nous servira à masquer les images que nous traiterons par la suite.

Il existe diférentes techniques pour extraire les nuages d’une image satellite. Selon les niveaux de traitements, certaines images sont fournies avec un raster représentant une bande de qualité (Qulity band). Sur ce raster, les valeurs des pixels renseignent sur la présence (ou l’absence) de nuages, d’ombres, de cirrus… C’est notamment le cas pour les images Landsat de niveau L2 (et même L1).

Nous verrons ici de façon pratique comment extraire les nuages d’une image selon différentes techniques et pour différentes images.

Nuages et Landsat dans QGIS 

Ici, nous allons extraire les nuages de notre image Landsat OLI du 2 avril 2021 prise au-dessus de la commune de Mbandjock, située dans la Haute-Sanaga au centre du Cameroun. Cette image présente de nombreux nuages comme nous pouvons le voir sur la composition colorée en fausses couleurs présentée ci-dessous (Fig. 322). Le but va être de repérer les nuages de cette image.

Il est tout d’abord nécessaire d’installer l’extension Cloud Masking. Pour des détails sur l’installation d’une extension QGIS, voire la section dédiée : Installation d’extensions dans QGIS. Une fois l’extension installée, il suffit de l’ouvrir via le menu Extensions ‣ Cloud masking for Landsat products ‣ Cloud Masking. Le module s’ouvre en bas à droite de l’interface de QGIS (Fig. 323).

Nous commençons par pointer le fichier texte des métadonnées de l’image à la ligne Landsat Metadata File (MTL) dans l’onglet Open and Load. Ces métadonnées permettent au module de savoir de quel Landsat il s’agît (TM, ETM+ ou OLI) et d’appliquer l’algorithme de masque de nuages correspondant. Une fois ces métadonnées renseignées nous cliquons sur load. Nous allons ensuite dans l’onglet Filters and Mask. Ici, nous pouvons choisir quelle technique appliquer. Trois techniques sont proposées : une basée sur la bande qualité fournie avec l’image (QA Band (C2)), une basée sur la bande bleue (Blue Band) et une basée sur un algorithme plus poussé nommé FMask. La technique basée sur la bande bleue ne donnant pas des résultats très bons, nous ne la présenterons pas ici.

Avec la bande de qualité 

La technique la plus simple est d’utiliser la bande de qualité fournie avec l’image. En fait, le module ne va finalement que recoder la bande de qualité avec des valeurs plus simples à manipuler et lui associer une symbologie. Nous cochons la case QA Band (C2) puis nous spécifions ce que nous souhaitons inclure dans notre masque. Ici nous allons y inclure les pixels correspondant aux Dilated Cloud, Cirrus, Cloud (nuage) et Cloud Shadow (ombre de nuage) (Fig. 324). La bande de qualité contient également un code pour les pixels en neige et en eau dans le cas où nous voudrions aussi exclure ces surfaces ci.

Nous cliquons ensuite sur Generate mask tout en bas de la fenêtre. Le raster des nuages (et des autres éléments sélectionnés précédemment) apparaît. Les nuages sont bien repérés mais les ombres ne sont pas toutes bien identifiées.

Astuce

Nous pouvons restreindre le masque à une zone d’intérêt soit à partir d’une couche de polygone de notre zone d’étude avec l’option In only polygon layer ou en dessinant à la volée une zone d’intérêt avec l’option In only areas of interest.

Avec l’algorithme FMask 

L’algorithme FMask est une méthode de détection des nuages et de leurs ombres parmi les plus robustes. FMask n’est pas basé sur la bande de qualité mais reconnait directement les nuages et leurs ombres à partir de l’analyse de l’image par des calculs relaitvement complexes. Pour utiliser cet algorithme, il suffit de cocher la case FMask de l’onglet Filters to apply. Dans cet algorithme, il est possible de régler 3 paramètres, mais les valeurs par défaut donnent généralement de bons résultats. Nous constatons que cet algorithme est également capable de détecter la neige et l’eau, mais nous n’utiliserons pas cette fonctionnalité ici (Fig. 325).

Nous cliquons ensuite sur Generate mask tout en bas de la fenêtre. Le raster des nuages (et des autres éléments sélectionnés précédemment apparaît). Les nuages sont bien repérés et leurs ombres aussi, ce qui est la vraie plus-value de cette méthode (Fig. 326).

Avertissement

L’algorithme étant assez complexe, il peut prendre un peu de temps à s’exécuter, surtout si nous l’exécutons sur une scène Landsat entière.

Nous constatons que FMask est plus généreux sur l’extension des nuages et qu’il repère bien mieux les ombres qu’en utilisant la bande de qualité seule. Par contre, dans les deux cas les cirrus sont difficilement identifiés.

Astuce

Si nous désirons être sûrs de n’avoir aucun nuage et aucune ombre de nuages, quitte à perdre des portions de l’image qui seraient nettes, nous pouvons augmenter le buffer des nuages et de leurs ombres dans les paramètres de FMask (Fig. 325).

Les indices radiométriques 

Grâce aux images satellites multi-spectrales il est possible de calculer des indices radiométriques en combinant plusieurs bandes spectrales. Un indice radiométrique est une combinaison de bandes spectrales via un calcul raster qui permet de mettre en exergue un état de la surface. Concrètement, à partir de deux, ou plus, bandes spectrales nous construisons un nouveau raster qui sera plus facile à interpréter en fonction de notre objet d’étude. Il existe des indices radiométriques dédiés à l’étude de la végétation, des surfaces en eau, du couvert neigeux, des zones incendiées, du sol nu, de la qualité des eaux … Il en existe énormément, mais seuls les plus fréquemment utilisés seront détaillés ici.

Certains logiciels ou certaines extensions proposent des modules tout prêts pour les calculer, mais la plupart du temps il est tout aussi simple, voire même plus sûr, de les calculer simplement à l’aide de la calculatrice raster de son logiciel de géomatique préféré (Calcul raster mono-bande).

Avertissement

L’ensemble de ces indices doivent se calculer à partir des bandes spectrales exprimées en réflectance de surface. Si les réflectances de surface ne sont vraiment pas disponibles, les réflectances au sommet de l’atmosphère (TOA) peuvent à la rigueur faire l’affaire. Mais il est déconseillé de calculer ces indices sur des bandes spectrales en niveau L1 qui ne seraient pas des réflectances. Voir la fiche Niveaux de traitements.

Indices dédiés à la végétation 

La végétation est sans doute l’état de surface le plus étudié via les images de télédétection. Le but de la plupart de ces indices est de construire un raster sur lequel les pixels aux plus fortes valeurs correspondront aux pixels couverts de végétation luxuriante. Il existe différents indices, mais l’indice roi est le NDVI.

NDVI 

NDVI est l’acronyme anglais de Normalized Difference vegetation Index. La traduction en français n’a pas trop de sens… C’est un indice mis au point dès les années 1970 et qui renseigne sur l’état de la végétation. C’est de loin l’indice le plus utilisé en télédétection de par sa simplicité de calcul, d’interprétation et sa robustesse. Il se calcule simplement comme suit :

\[NDVI = \frac{(PIR - R)}{(PIR + R)}\]

où PIR correspond à la bande spectrale du Proche infrarouge et R correspond à la bande spectrale du Rouge.

Cette formulation est basée sur la signature spectrale typique de la végétation. Pour bien comprendre le principe, nous nous appuierons sur la figure Fig. 327. Nous nous focaliserons sur les signatures spectrales de la végétation (Surface herbacée) et celle de la neige. Sur cette figure, nous avons également ajouté les domaines du Rouge et du Proche infrarouge.

Sur cette figure, nous relevons une réflectance de 0.02 et de 0.5 pour la surface herbacée dans le Rouge et le Proche infrarouge respectivement. Pour la neige nous relevons respectivement 0.98 et 0.92 pour le Rouge et le Proche infrarouge. Connaissant ces valeurs, il est possible de calculer le NDVI pour ces deux états de surface :

\[NDVI végétation = \frac{(0.5 - 0.02)}{(0.5 + 0.02)} = 0.92\]

\[NDVI neige = \frac{(0.92 - 0.98)}{(0.92 + 0.98)} = -0.03\]

Nous constatons que les surfaces en végétation présentent un NDVI élevé et que les surfaces enneigées un NDVI négatif. Sur le même principe, les surfaces en eau libre présentent également un NDVI négatif généralement. D’une manière générale, le raster résultat est borné entre -1 et 1. L’interprétation d’un raster de NDVI est la suivante : plus un pixel tend vers 1, plus ce pixel est couvert de végétation bien développée et en bonne santé. Un pixel qui tend vers 0 est un pixel de sol nu. Enfin, un pixel négatif est généralement un pixel en eau. Mais attention, dans certains cas, si des algues ou une forte turbidité sont présents dans les zones en eau, alors cette eau pourra avoir un NDVI légèrement positif du fait de l’activité chlorophyllienne des algues et du phytoplancton. La figure Fig. 328 présente le NDVI calculé sur les Bouches-du-Rhône le 13 juin 2019 à partir d’une image Landsat 8.

SAVI 

SAVI est l’acronyme anglais de Soil Adjusted Vegetation Index. Tout comme le NDVI, le SAVI donne une indication sur la présence ou l’absence de végétation pour tout pixel. Une forte valeur de SAVI traduit un pixel couvert de végétation en bonne santé et une faible valeur traduit plutôt un sol sans végétation. Le SAVI introduit un facteur correctif qui prend en compte la présence de sol nu au sein de la végétation. Par exemple, dans le cas d’une forêt éparse ou d’une savane sèche, où les arbres se trouvent séparés les uns des autres par du sol nu, le SAVI peut produire des résultats rendant mieux compte du couvert végétal que le NDVI. Le SAVI se calcule comme suit :

\[SAVI = \frac{(1 + L)(PIR - R)}{(PIR + R + L)}\]

où PIR correspond à la bande spectrale du Proche infrarouge, R correspond à la bande spectrale du Rouge et L correspond au facteur correctif de sol nu.

Le facteur L peut varier entre 0 et 1 selon les configurations mais il est admis qu’une valeur de 0.5 donne de bons résultats en règle général. La figure Fig. 329 présente le SAVI calculé sur les Bouches-du-Rhône le 13 juin 2019 à partir d’une image Landsat 8.

Indices dédiés aux sols nus 

Les sols nus sont une occupation du sol qui peut être intéressante à suivre à plusieurs égards. Un sol nu peut être naturel, dans le cas d’une région aride par exemple, mais peut également renseigner sur des dynamiques environnementales. Après un incendie ou une inondation ayant décapé la végétation, un sol nu peut être une trace de l’événement. Dans un contexte agricole, un sol nu peut faire partie d’une rotation culturale en étant le stade post récolte et pré germination. Il existe différents indices pour mettre en avant un sol nu. Quelques exemples sont donnés ici.

BSI 

BSI est l’acronyme anglais de Bare Soil Index (Indice de sol nu). Il se calcule par combinaison des bandes des domaines du bleu, du rouge, du proche infrarouge et du moyen infrarouge. Sa formulation mathématique est présentée ci-dessous.

\[BSI = \frac{(MIR + R)-(PIR + B)}{(MIR + R)+(PIR + B)}\]

où MIR est la bande du Moyen infrarouge (B6, B5 ou B11 pour respectivement les capteurs OLI/OLI2, TM et MSI de Sentinel-2), R la bande du domaine du Rouge, PIR la bande du domaine du Proche infrarouge et B la bande du domaine du Bleu.

Sur le raster résultat, les pixels ayant une valeur positive peuvent être considérés comme étant des pixels de sol nu. Un exemple est montré sur une scène Landsat 8 prise au-dessus des Bouches-du-Rhône le 13 juin 2019 (Fig. 329). Les pixels en rouge ont un BSI positif et sont donc des pixels de sol nu.

En comparant le résultat de cet indice à un indice de végétation comme le NDVI (Fig. 328), il s’agit presque d’une image symétrique, ce qui est tout à fait cohérent.

Indices dédiés aux zones incendiées 

Il existe des indices dédiés au repérage des des zones incendiées et dédiées à l’analyse de la sévérité des incendies.

NBR 

NBR est l’acronyme anglais de Normalized Burn Ratio et est un indice utilisé pour cartographier les zones touchées par un incendie. L’indice repose sur la signature spectrale de la végétation et la signature spectrale des zones incendiées. La végétation en bonne santé présente un pic dans le domaine du proche infrarouge (à peu près 0.75 µm) et un creux dans celui de l’infrarouge moyen (à peu près 2.30 µm). Au contraire, les surfaces récemment incendiées présentent un creux dans le domaine du proche infrarouge et un pic dans celui de l’infrarouge moyen (Fig. 331).

En se basant sur ces signatures spectrales, le NBR se calcule comme suit :

\[NBR = \frac{(PIR - MIR)}{(PIR + MIR)}\]

Avec PIR la bande spectrale du domaine du Proche InfraRouge et MIR celle du domaine du Moyen InfraRouge. Dans le cas de Landsat 8 ou Sentinel-2, la seconde bande du moyen infrarouge, à savoir la MIR2 sera privilégiée.

À la suite de ce calcul, le raster obtenu varie entre -1 et +1. Les pixels négatifs correspondent à des zones récemment incendiées tandis que les pixels positifs correspondent à une végétation en bonne santé.

difference NBR 

L’indice dNBR - diffference NBR - aussi noté ΔNBR, est un indice utilisé pour quantifier la sévérité des incendies. Il se calcule simplement par différence des indices NBR calculés pour une date avant incendie et pour une date après incendie. Concrètement, il s’agît de la soustraction des deux rasters de NBR calculés pour les deux dates.

\[ΔNBR = NBR avant Incendie - NBR après Incendie\]

Le raster résultat varie autour de 0. Les valeurs négatives reflètent généralement une repousse de la végétation et les valeurs positives reflètent la sévérité de l’incendie. Les valeurs sont normalement à interpréter au cas par cas avec une étude de terrain, mais l’USGS propose une échelle de sévérité des incendies qui fonctionne dans la plupart des cas.

Tableau 12 dNBR et sévérité des incendies
Niveau de sévérité	gammes de ΔNBR
Repousse de la végétation (élevée)	-2.000 à -0.250
Repousse de la végétation (modérée)	-0.250 à -0.100
Pas incendié	-0.100 à +0.100
Incendié à faible sévérité	+0.100 à +0.270
Incendié à sévérité modérée	+0.270 à +0.440
Incendié à sévérité élevée	+0.440 à +0.660
Incendié à sévérité très élevée	+0.660 à +2.000

Indices dédiés à l’eau 

Les surfaces en eau sont également un important sujet d’étude de la télédétection. C’est pourquoi il existe plusieurs indices de détection des surfaces en eau, plus ou moins élaborés.

MNDWI 

MNDWI est l’acronyme anglais de Modified Normalized Water Index. C’est un indice largement utilisé pour détecter les surfaces en eau depuis une image satellite multispectrale. Il se calcule comme suit :

\[MNDWI = \frac{(V - MIR)}{(V + MIR)}\]

où V correspond à la bande du domaine du Vert et MIR correspond à la bande du domaine du Moyen Infrarouge. Dans le cas de Landsat 8, il est conseillé de prendre la bande 6 correspondant au MIR 1 pour calculer cet indice.

Avertissement

La formule du MNDWI peut varier quelque peu selon les auteurs.

Au final, le raster obtenu varie entre -1 et 1. Sur un raster de MNDWI, les pixels positifs correspondent aux pixels en eau et les pixels négatifs aux pixels de terre. Il est donc facile de faire la différence entre les zones en eau et les zones en terre. La figure Fig. 332 présente le MNDWI calculé sur les Bouches-du-Rhône le 13 juin 2019 à partir d’une image Landsat 8.

AWEI 

AWEI est l’acronyme anglais de Automated Water Extraction Index. Il a notamment été développé pour être appliqué sur les images Landsat disposant de deux bandes spectrales dans le moyen infrarouge. Cet indice a la particularité d’avoir deux versions. Une version AWEIsh est dédiée à l’extraction des zones en eau dans les régions présentant beaucoup d’ombres sur les images, comme les régions de montagnes. Une seconde version AWEIns est dédiée à l’extraction des zones dans les régions présentant peu d’ombres, comme les régions de plaines et de plateaux.

Avertissement

De par sa formulation, le AWEI ne peut être calculé que sur des bandes spectrales en réflectance exprimées entre 0 et 1.

Cet indice se calcule comme suit pour sa version appliquée aux images présentant des ombres (shadow) :

\[AWEIsh = B + 2.5 × V − 1.5 × (P IR + MIR1) − 0.25 × MIR2\]

Cet indice se calcule comme suit pour sa version appliquée aux images sans ombres (non shadow) :

\[AWEIns = 4 × (V − MIR1) − (0.25 × P IR + 2.75 × MIR2)\]

où :

B correspond à la bande spectrale du domaine du Bleu
V correspond à la bande spectrale du domaine du Vert
PIR correspond à la bande spectrale du domaine du Proche Infrarouge
MIR1 correspond à la bande spectrale du domaine du Moyen Infrarouge 1
MIR2 correspond à la bande spectrale du domaine du Moyen Infrarouge 2

Astuce

La version sans ombres est plus couramment employée que celle avec ombres.

Sur le raster de AWEI résultat, un pixel en eau aura une valeur positive et un pixel de terre aura une valeur négative. Ce seuil est pratique pour la détection automatisée des surfaces en eau. Il présente généralement des résultats un peu plus précis que les autres indices d’extraction des surfaces en eau. La figure Fig. 333 présente le AWEI (sans ombres) calculé sur les Bouches-du-Rhône le 13 juin 2019 à partir d’une image Landsat 8.

Indices dédiés à la neige 

Des indices dédiés à l’extraction des zones enneigées depuis les images de télédétection ont été développées depuis longtemps. Nous présenterons ici l’indice le plus utilisé.

NDSI 

NDSI est l’acronyme anglais de Normalized Difference Snow Index. Il s’agît d’un indice simple à calculer en utilisant l’équation suivante.

\[NDSI = \frac{(V - MIR)}{(V + MIR)}\]

où V correspond à la bande du domaine du Vert et MIR correspond à la bande du domaine du Moyen Infrarouge. Dans le cas de Landsat 8, il est conseillé de prendre la bande 6 correspondant au MIR 1 pour calculer cet indice.

Note

Nous pouvons noter que le NDSI se formule de la même manière que le MNDWI dédié à l’eau. La distinction de la neige de l’eau libre se fera en trouvant un seuil de NDSI au-dessus duquel nous serons quasiment certains de n’extraire que les surfaces en eau.

Sur la figure précédente nous constatons que les surfaces de terre ont des valeurs de NDSI négatives et que les surfaces en eau libre ou en neige présentent des valeurs de NDSI positives. Cependant nous constatons que les surfaces en neige présentent un NDSI bien plus élevé (au-dessus de 0.5 généralement) que l’eau (qui reste au-dessous de 0.1 dans la plupart des cas).

Outils dédiés 

Même si il est très simple d’utiliser le calcul raster (Calcul raster mono-bande), via par exemple la calculatrice raster de QGIS (Calcul raster avec QGIS) ou de R (Calcul raster avec R), pour calculer ces différents indices, il existe des outils dédiés. Ces outils peuvent être contenus dans des logiciels ou peuvent faire partie de librairies de langages de programmation comme R ou Python.

Indices avec OTB 

Version de Orfeo ToolBox : 7.2.0

OTB (Orfeo ToolBox) propose un menu dédié au calcul de différents indices. Ce menu se nomme RadiometricIndices et s’utilise soit via l’interface Monteverdi soit via la ligne de commande otbcli.

Via Monteverdi

Le menu se trouve dans le panneau Navigateur d’OTB-Applications ‣ Feature extraction ‣ RadiometricIndices. Le menu suivant apparaît (Fig. 335).

Dabs le champ Input image nous indiquons le raster à utiliser pour le calcul. Il doit nécessairement s’agir d’un raster multi-bandes (Raster multi-bandes) (en dur car les rasters virtuels ne sont pas acceptés) contenant au moins les bandes spectrales nécessaires au calcul de l’indice en question. Ici, nous travaillons avec une scène Landsat 8 constituée d’un raster multi-bandes contenant les bandes 2 à 6 sur les Bouches-du-Rhône. Dans le champ Output image, nous spécifions un chemin et un nom pour le raster qui sera créé, ici nous allons calculer un NDVI. Dans le panneau Channels selection, nous indiquons quelles sont les bandes qui correspondent aux différents domaines radiométriques.

Note

Il n’est pas indispensable de renseigner les 5 canaux. Si le raster multi-bandes ne dispose que d’une bande rouge (1ère bande) et proche-infrarouge (2nde bande), nous indiquons pour le Red Channel 1 et pour le NIR Channel 2, puis des valeurs fictives, 3 par exemple, pour les autres canaux.

Enfin, dans le panneau Available Radiometric Indices, nous sélectionnons le ou les indices à calculer. Puis nous cliquons sur Execute pour lancer le calcul. Si nous avons sélectionné plusieurs indices, le raster résultat sera un raster multi-bandes dans lequel la première bande correspondra au premier indice et ainsi de suite.

Via otbcli

Le menu RadiometricIndices est également utilisable via la ligne de commande otbcli, ce qui peut s’avérer pratique dans une optique d’automatisation de tâches. Pour l’exemple, nous allons calculer un NDVI et un (MNDWI) à partir d’un raster multi-bandes Landsat 8 contenant les bandes 2 à 6. La commande à employer est la suivante :

otbcli_RadiometricIndices -in stack_LC08_L2SP_196030_20190613.TIF -channels.blue 1 -channels.green 2 -channels.red 3 -channels.nir 4 -channels.mir 5 -list Vegetation:NDVI Water:MNDWI -out RadiometricIndicesImage.tif

Où :

otbcli_RadiometricIndices : la commande à utiliser
-in : le mot clé pour spécifier le raster multi-bandes à utiliser
-channels.blue 1 : la précision du numéro de bande à utiliser pour le canal du Bleu, et ainsi de suite pour tous les canaux
-list Vegetation:NDVI Water:MNDWI : la liste des indices à calculer, se reporter à la documentation pour la liste complète
-out : le mot clé pour spécifier le chemin et le nom du raster qui sera créé