Tutoriel classification d’images avec QGis: 2.1 Correction atmosphérique des images Landsat 8

Nous avons vu précédemment comment corriger les images Landsat pour avoir la réflectance TOA (top-of-atmosphere). Voilà un schéma qui va nous aider à comprendre la suite:

La réflectance que nous avons calculé est donc le % de lumière réfléchie par rapport à la totalité de la lumière visible incidente. Mais comme on le voit dans le schéma, le capteur du satellite mesure deux choses en même temps: la lumière réfléchie par nos cibles à la surface de la terre, plus la lumière diffusée par les particules en suspension dans l’air.

On peut pousser la correction atmosphérique des images satellites pour enlever la partie de lumière due à la diffusion.

Comme nous l’avons déjà dit, ceci n’a d’intérêt que si vous travaillez sur des images séparées dans le temps. Le pourcentage de diffusion étant le même pour l’ensemble d’une image, c’est du temps perdu de faire cette correction.



Le mode d’emploi pour effectuer cette correction sur les images Landsat est disponible à l’adresse http://www.gisagmaps.com/landsat-8-atco-guide/ 

Mais si vous êtes débutant, vous aurez du mal à comprendre la marche à suivre. Nous allons la voir ici, en utilisant les outils de QGis.

Methode DOS (
Dark Object Subtraction)

La méthode que nous allons voir est la plus couramment employée.
La soustraction d’objets noirs est une méthode de correction atmosphérique empirique simple pour l’imagerie satellitaire, qui suppose que la réflectance des objets sombres inclut une composante importante de la diffusion atmosphérique.

Le concept est simple: si on a un objet complètement noir, et donc qui ne réfléchit pas du tout de lumière, cet objet qui devrait avoir une valeur de 0 réflectance va avoir une valeur non nulle au niveau du capteur satellite. Cette valeur sera la réflectance due aux particules présentes dans l’atmosphère.

La soustraction d’objets sombres recherche dans chaque bande la valeur de pixel la plus sombre. La diffusion est supprimée en soustrayant cette valeur de chaque pixel de la bande.

Première chose: la correction que nous allons faire ne s’applique qu’à la lumière visible et le proche infra-rouge. Nous travaillerons donc, pour les images Landsat 8, sur les bandes 2 (Bleu), 3 (vert), 4 (Rouge) et 5 (proche infra-rouge).

Deuxième chose à savoir: le principe utilisé est de calculer la diffusion pour la bande 4 (rouge) et d’utiliser une formule de correspondance pour déterminer la correction pour les autres trois bandes.

Troisième chose: nous allons travailler avec les bandes 2,3,4 et 5 en réflectance TOA mais aussi en radiance brute (DN) pour la bande 4 (rouge). Dans le mode d’emploi cité plus haut ceci n’apparaît pas du premier coup d’œil. Voyons la procédure pas à pas:

1- on calcule la réflectance TOA pour les quatre bandes avec la formule

TOAbande =((DN * 0.00002)-0.1)/sin(élévation solaire)

2- en analysant les valeurs de radiance brute (DN) de la bande rouge, nous allons déterminer une radiance due à la diffusion atmosphérique (DNred)

3- nous transformons cette radiance en réflectance TOA

CORred=((DNred * 0.00002)-0.1)/sin(élévation solaire)

4- avec un abaque on entre cette diffusion de la bande rouge et on obtient la diffusion pour les bandes 2,3 et 5 (CORbleu, CORvert, CORir)

5-on soustrait les valeurs obtenues des réflectances de TOA pour obtenir les valeurs de réflectance des 4 bandes corrigées de l’effet atmosphérique.

Rrouge= TOA(rouge)-CORred

Rbleue= TOA(bleu)-CORbleu

Rvert= TOA(vert)-CORvert

Rir= TOA(ir)-CORir

Voyons en détail les étapes 2 et 4.

Calcul de la diffusion de la bande rouge

Il y a plusieurs méthodes pour déterminer la correction à apporter à la bande rouge. La méthode conseillée est celle dite de la Bin 5 – 0.008. Vous trouverez ici une description détaillée de la méthode.

En pratique, avec QGis, noius allons utiliser l’histogramme de la bande pour déterminer la valeur de seuil de la diffusion atmosphérique.

1- Vous devez d’abord exécuter l’outil Statistiques de la couche raster sur les données de la bande 4 (rouge).

Une fois exécuté, vous aurez un ligne de résultat dans le panneau Visualiseur de résultats

Double-cliquez sur cette ligne pour ouvrir les résultats dans votre navigateur:

Notez la valeur minimale, ici 3916.

Nous allons rechercher la valeur à affecter à l correction de la diffusion atmosphérique dans l’histogramme de l’image.

Cliquez sur la couche de l’image -> Propriétés -> Histogramme

Entrez la valeur minimale trouvée dans les statistiques (ici 3916) dans la zone Min, puis cliquez sur la zone Max ci-dessous pour déplacer la valeur minimale de la ligne pointillée sur l’histogramme.

Vous pouvez utiliser le format de ligne pour trouver la base de l’histogramme, mais pour cet exemple, le format de barre sera utilisé et il est recommandé. Si l’histogramme est au format ligne, cliquez sur Préférences / Actions et décochez la case Tracer comme lignes.

Un histogramme au format barre apparaîtra. L’histogramme QGIS ne représente pas / ne montre pas toutes les valeurs comme l’histogramme ArcGIS, mais il affiche bien la forme qui peut être appropriée pour établir une valeur de dispersion qui représente la base de l’histogramme (en cas d’augmentation relativement abrupte des fréquences).  Les étapes suivantes montrent comment l’histogramme QGIS peut être utilisé pour établir la dispersion.

Passez le curseur de la souris sur l’histogramme pour afficher une loupe. Faites glisser la loupe sur l’extrémité inférieure de l’histogramme pour lancer le zoom sur la base de l’histogramme. Si vous effectuez un zoom trop important ou si vous souhaitez recommencer pour une raison quelconque, cliquez sur Préférences / Actions et Recalculer l’histogramme. L’histogramme est zoomé ci-dessous.

Vous verrez que la barre avec la valeur la plus basse est encore assez supérieure à la valeur d’image basse que vous avez précédemment entrée dans la zone Min (3619); l’histogramme QGIS ne montre pas beaucoup de valeurs dans une queue statistique (Landsat 8 a généralement de longues queues statistiques). La base de l’histogramme est la zone utilisée pour établir la dispersion.En zoomant sur les trois premières barres de l’histogramme:

Dans cet exemple, les valeurs comprises entre environ 5570 et 5595 peuvent se rapprocher de la base de l’histogramme – ce sont des valeurs très basses et ne sont pas trop sur la queue statistique (ce processus peut être quelque peu subjectif). Pour les besoins de ce document, le DN 5569 sera utilisé comme base de l’histogramme. L’histogramme sera zoomé une fois de plus pour cet exemple.

La règle à utiliser est de prendre comme valeur de dispersion la valeur minimum de la première barre dépassant une fréquence de 5 (d’où le nom Bin 5) et n’ayant pas de fréquences inférieures à 5 à sa droite. Dans l’exemple suivant:

Les deux premières barres dépassent la valeur 5, mais comme la troisième est inférieure, c’est la valeur min de la quatrième barre qui sera retenue.

Dans notre exemple, la valeur de radiance à retenir pour la correction atmosphérique de la bande rouge sera 5569.

On peut passer à la troisième étape,

3- nous transformons cette radiance en réflectance TOA

CORred=((DNred * 0.00002)-0.1)/sin(élévation solaire)

 CORred= (( 5569 * 0.00002)-0.1)/0.42631886  =0.026694

Cette valeur correspond à la Bin 5. Comme la méthode est Bin 5 – 0.008, nous enlevons 0.008 de cette valeur: 0.018694

C’est cette valeur qu’il faut soustraire aux valeur de TOA du rouge précédemment calculées pour avoir les valeurs corrigées de la diffusion atmosphérique.

Correction des autres bandes

A partir de la valeur de correction du rouge on obtient les valeurs à appliquer à chacune des autres bandes:

Allez sur la page http://www.gisagmaps.com/l8-s2-relative-scatter-calc/

Vous trouverez une fenêtre où vous pouvez rentrer la valeur de correction du rouge

Rentrez votre valeur (ici 0.018647), vous aurez automatiquement les valeurs des autres corrections qui vont s’afficher

Nous obtenons ainsi les corrections de

  • 0.06309 pour la bande Bleue (2)
  • 0.03442 pour la bande Verte (3)
  • 0.00626 pour la bande infrarouge proche (5)

Pour avoir les bandes corrigées on utilise la calculatrice raster en soustrayant à chaque bande corrigée TOA la constante de correction que nous avons trouvé.

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *