De l'importance d'un bon masque de nuage pour le traitement automatisé de séries temporelles

=>

Série temporelle de réflectance au sommet de l'atmosphère pour les 4 bandes à 10m de résolution de Sentinel-2, issues des produits L1C

 

Le graphique ci-dessus montre une série temporelle de réflectance TOA rassemblée par Sentinel-2 sur un pixel choisi au hasard dans une tuile au centre de la France (tuile 31TDK, pixel 3000-7000), à partir de produits L1C. En regardant la série chronologique, il est assez difficile de dire quel type de surface a été observé, même si un cycle végétatif semble être présent. Comme nous le verrons ci-dessous, la plupart du bruit observé est dû à la présence de nuages ​​ou d'ombres de nuages.

 

La courbe ci-dessous montre qu'après avoir retiré tous les nuages ​​et leurs ombres, la réflectance au sommet de l'atmosphère est déjà plus lisse, et il est ainsi beaucoup plus facile de comprendre le type de surface observée. Continue reading

MACCS/MAJA, comment ça marche ?

=>

MACCS (Multi-sensor Atmospheric Correction and Cloud Screening) est une chaîne de traitement de niveau 2A, qui, comme son nom l'indique, détecte les nuages et leurs ombres, estime l'épaisseur optique des aérosols, la quantité de vapeur d'eau, et corrige les effets atmosphériques. La chaîne a été développée conjointement par le CESBIO et le CNES. Le CESBIO a mis au point les méthodes et développé un prototype, tandis que le CNES a pris en charge la version opérationnelle de la chaîne, que le CESBIO a largement contribué à valider. Le développement de la chaîne opérationnelle a été confié par le CNES à la compagnie CS-SI, et depuis peu, sa validation est gérée par CAP Gemini et le CNES.

 

Plus récemment, le CNES+CESBIO et le DLR ont décidé de mettre en commun leurs efforts pour développer la chaîne MAJA (MACCS-ATCOR Joint Algorithm). Cette chaîne est une évolution de la chaîne MACCS dans laquelle des méthodes issues de la chaîne ATCOR du DLR seront ajoutées progressivement. MAJA V1.0 a MACCS aurait dû s'appeler MACCS V6.0, mais le CNES et le DLR ont préférer changer de nom pour célébrer leur entente dans ce domaine.

 

MAJA est exploité au CNES dans le segment sol MUSCATE de Theia (données disponibles ici) et dans le segment sol du satellite Venµs. Enfin MAJA est disponible gratuitement en version binaire pour utilisation non commerciale.

 

MAJA a la particularité de s'appliquer à des séries temporelles d'images de résolution décamétrique, et de faire une large utilisation de méthodes multi-temporelles. Pour cette raison, elle ne s'applique qu'aux mission dont les observations se font sous des angles quasi constants. Ceci dit, MAJA a déjà été appliqué à de nombreux satellites :

 

Macroscopiquement, le fonctionnement de MAJA est expliqué dans le schéma ci-contre.

Correction de l'absorption atmajasphérique

Dans le cas de Sentinel-2 et de Venµs, qui disposent d'un canal 940 nm (resp 910 nm), situé dans une forte bande d'absorption de la vapeur d'eau, une première étape estime la quantité de vapeur d'eau dans  la colonne d'air traversée par l'observation. Pour les autres satellites, ce sont des données météo qui sont utilisées. Puis on passe à la correction de l'absorption gazeuse. Cette phase est réalisée avec le code SMAC.

Imajage composite

Les étapes suivantes font un grand usage de méthodes multi-temporelles. Et pour cela, une série temporelle doit être traitée dans l'ordre chronologique. En sortie de chaque traitement, une image composite est mise à jour avec les derniers pixels non nuageux acquis. Cette image sert de référence au traitement de détection des nuages et d'estimation de l'épaisseur optique des aérosols.

Les majasques de nuages

La détection des nuage repose sur une batterie de tests, dont les plus importants et plus efficaces, sont :

  • celui qui utilise la bande Cirrus (à 1380 nm) présente sur Landsat 8 et Sentinel-2 qui détecte les nuages hauts,
  • le test multi-temporel, qui détecte une soudaine augmentation de la réflectance dans le bleu, signe de la présence d'un nuage.
  • enfin, pour éviter des sur-détections, pour chaque nuage potentiellement détecté, un dernier test mesure le taux de corrélation d'un voisinage du pixel avec les images précédentes. Comme il est peu probable que plusieurs nuages successifs se ressemblent, au même endroit, si une bonne corrélation est obtenue avec l'une des images, le pixel n'est finalement pas classé nuageux

Après la détection des nuages, il faut aussi détecter leurs ombres, l'eau, et la neige.

Estimajation de l'épaisseur optique des aérosols

L'estimation de l'épaisseur optique combine plusieurs critères, décrits dans cette page, dans le calcul d'une fonction coût, qui est ensuite inversée par la méthode des moindres carrés non linéaires.

  • Critère multi-temporel ; après correction atmosphériques, deux voisinages successifs non nuageux (celui de l'image en cours, et celui du satellite) doivent avoir quasiment les mêmes réflectances. Les écarts résiduels après correction atmosphérique sont inclus dans la fonction coût.
  • Critère multi-spectral ; au dessus de la végétation, et aussi au dessus de nombreux sols nus, la réflectance de surface dans le bleu est proche de la moitié de la réflectance dans le rouge. Les écarts à cette relation après correction atmosphérique sont donc inclus aussi dans la fonction coût
  • Minimum ou maximum de l'épaisseur optique : l'épaisseur optique ne peut pas être négative, et elle ne peut normalement pas dépasser la valeur estimée par la méthode du pixel sombre. Lorsque les erreurs dépassent ces seuils, une forte erreur est ajoutée dans la fonction coût.

 

La minimisation de cette fonction coût est réalisée sur un voisinage de pixels à 240 m,  s'étendant sur environ deux kilomètres. Les résultats obtenus sont ensuite lissés, et les trous dûs au nuages sont bouchés, pour obtenir finalement une carte d'aérosols à une résolution de 5 km. Le type d'aérosols n'est pas estimé, c'est un paramètre de la méthode qui peut être fixé par zone géographique.

Correction atmajasphérique

 

L'un des quicklooks que nous produisons systématiquement à chaque traitement de MACCS, pour vérifier le bon fonctionnement d'un coup d'oeil, ici pour le site SPOT5 (Take5) du Chiapas, au Mexique. En haut à gauche, l'image au sommet de l'atmosphère, en bas à gauche, l'AOT et le masque de nuages, en haut à droite, les réflectances de surface corrigées des effets d'environnement,,en bas à droite, les mêmes, corrigées aussi des effets de pente.

Une fois connue l'épaisseur optique des aérosols, nous pouvons calculer les réflectances de surface.  Pour celà, nous utilisons des tableaux précalculés, à partir du code de transfert radiatif SOS (Successive Orders of Scattering, Lenoble, 2007). Ces mêmes tableaux sont aussi utilisés pour l'estimation des aérosols.

Les réflectances ainsi obtenues peuvent être utilisées pour mettre à jour l'image composite.

Avant d'éditer le produit de sortie, il reste cependant à corriger deux autres points, déjà détaillés dans ce blog,  les effets d'environnement et les effets de pentes en présence de relief.

 

Le développement de MACCS a démarré depuis plus de 10 ans, et de très nombreuses personnes y ont contribué au cours de ces années :

  • au CESBIO (H.Tromp, V. Debaecker, M. Huc, P.Gely, B.Rouquié et O.Hagolle),
  • au CNES (B. Petrucci, D.Villa-Pascual, Camille Desjardins),
  • au DLR (A.Makarau, R.Richterà
  • Chez CS-SI (T.Feuvrier, C.Ruffel, A.Bricier et plusieurs autres personnes)
  • Chez Cap Gemini (M.Farges, G.Rochais, E.Durand)
  • Chez Magellium (E. Hillairet)

 

Nous avons publié quatre articles dans des revues à comité de lecture, sur les méthodes et la validation de MACCS :

L'USGS distribue (enfin) les données LANDSAT 8 corrigées des effets atmosphériques

=>

Depuis Noêl ou presque (le 23/12),  l'USGS distribue aussi des données LANDSAT 8 corrigées des effets atmosphériques, comme c'était déjà le cas pour LANDSAT 5 et 7.Ces données étaient très attendues et devraient connaître un grand succès.

 

Bien évidemment, les traitements de l'USGS portent sur le monde entier. Pour y accéder, il faut utiliser le serveur earthexplorer. Les données peuvent être trouvées en cliquant sur LANDSAT CDR (Climate Data Record), comme montré sur l'image ci-jointe. Les données ne sont pas en ligne et doivent être commandées.

 

Nous n'avons pas encore fait de comparaison entre les produits de Niveau 2A de LANDAT 8 mis à disposition par THEIA sur la France avec la chaîne du CESBIO, et les traitements de l'USGS développés par la NASA (E.Vermote), mais nous y travaillerons prochainement.

 

Les données LANDSAT 5 et 7 acquises sur la France de 2009 à 2011 bientôt disponibles au niveau 2A.

=>

Comme pour LANDSAT 8, il y a quelques semaines, nous venons de produire les données de LANDSAT 5 et LANDSAT 7 acquises au dessus de la France métropolitaine (sans la Corse) de 2009 à 2011. Les données seront mises en ligne d'ici quelques jours, le temps de les transférer sur le serveur.  Les traitements ont été effectués par l'équipe MUSCATE, pour le compte du pôle THEIA, en utilisant les moyens du centre informatique du CNES. Les données seront disponibles sur le site :

http://spirit.cnes.fr/resto/Landsat/

Exemple de produit de Niveau 2A obtenu avec LANDSAT5 sur la côte atlantique. Les nuages sont entourés en vert, le masque d'eau est entouré en bleu, la neige en rose. De temps en temps, la forte turbidité des eaux charentaises fait basculer les seuils de neige, et certains zones particulièrement turbides peuvent être identifiées dans le masque de neige...

Le traitement, le format et la présentation de ces données ont beaucoup de points communs avec ceux de Landsat 8, décrits ici. Cependant, ils s'en écartent par plusieurs aspects, qui font l'objet de cet article.

Point de départ.

Le point de départ est différent pour les satellites LANDSAT 7 et LANDSAT 5.

  • Pour LANDSAT 7, comme pour LANDSAT 8, nous démarrons des produits de niveau 1T fournis par l'USGS. Ces produits présentent cependant un gros défaut, avec la présence de stries dès que l'on s'écarte du centre de l'image. Ces stries sont dues à la panne d'un miroir sur ce satellite depuis 2013. Se reporter ici pour une description du problème. Dans notre cas, nous avons décidé de n'utiliser que la partie centrale de l'image, en interpolant un peu les données manquantes, tant que les trous sont inférieurs à 4 pixels, et en éliminant directement la partie de l'image où les trous ont une largeur supérieure à 4 pixels.
  • Les données acquises par LANDSAT 5 en France ne sont malheureusement pas disponibles sur les serveurs de l'USGS. C'est l'ESA qui dispose de ces données, et qui a bien voulu nous les fournir (Merci encore à Bianca Hoersh et  Riccardo Biasutti de l'ESA, et à la société SERCO qui nous a fait parvenir les données). Le jeu de données que nous vous fournissons ici est en fait un jeu unique, qui actuellement n'est disponible sur aucun autre serveur, même si l'ESA compte en faire une production prochaînement. Ces données nous parviennent au niveau 1G, un niveau de traitement intermédiaire, pour lequel les données n'ont pas été ortho-rectifiées. Nous avons donc dû les ortho-rectifier à Theia, en utilisant l'outil SIGMA du CNES, comme pour SPOT4 (Take5).

Cette  double approche différente pour les deux capteurs a des inconvénients. La base de données de référence géométrique de LANDSAT semble avoir quelques erreurs en France, et les données obtenues par l'USGS sur la région Toulousaine par exemple sont souvent décalées d'un pixel. Ce n'est pas le cas pour les données LANDSAT 5 ortho-rectifiées par THEIA, il peut donc parfois y avoir un pixel d'écart d'une date à l'autre selon qu'elle vient de LANDSAT 5 ou LANDSAT 7. Par ailleurs, les données de l'ESA présentent des défauts, comme par exemple la présence de points brillants colorés par-ci par là, de manière aléatoire.

 

 

Reprojection en Lambert 93

Les données de l'USGS sont fournies en projection UTM. Cette projection utilise trois fuseaux différents au dessus de la France, qui se divise donc en 3 zones différentes, l'Ouest de la France, le Centre et l'Est. Comme les données de deux fuseaux différents ne se superposent pas directement, nous avons donc décidé de reprojeter les données LANDSAT 7 en Lambert 93, qui est la projection officielle pour la France.Les données LANDSAT5 ont elles aussi, bien sûr, été directement projetées en Lambert 93.

Découpage des produits en tuiles

Nous avons pris le parti de suivre la même logique que celle utilisée par Sentinel-2, et de découper les données en tuiles de 110*110 km décalées de 100 km les unes par rapport aux autres. La tuile 1x1 se trouve au sud-ouest de la France, lorsqu'on va vers l'est (vers la droite), on incrémente la première coordonnée de D0001 à D0010 (D pour Droite), lorsqu'on va vers le Nord (vers le haut), on incrémente la seconde coordonnée, de H0001 à H0010 (H pour Haut). La tuile de Toulouse s'appelle donc D0005H0002 Le découpage en tuiles est visible sur l'image ci-jointe.

 

Pour chacune des tuiles, nous fournissons l'ensemble des dates pour lesquelles une image LANDSAT (5 ou 7) a une intersection non nulle avec la tuile. Quelques dates peuvent manquer, pour plusieurs raisons, liées en général à la couverture nuageuse :

  • l'image n'a pas été acquise par LANDSAT (quand les prévisions météo indiquent un temps très couvert, les images ne sont pas acquises).
  • l'image a été acquise mais s'est avérée trop nuageuse pour être traitée au niveau 1T par l'USGS (L7) ou par le CNES (L5)
  • l'image est trop nuageuse pour être traitée par la chaîne de Niveau 2A

 

Traitement de Niveau 2A (correction atmosphériques et détection des nuages)

Il est important de noter que notre chaîne ne traite pas les bandes thermiques pour l'instant. Une correction est à l'étude, mais celle-ci ne sera pas opérationnelle avant un ou deux ans.

La méthode utilisée pour les bandes visible, proche et moyen infra rouge est quasiment la même que pour SPOT4(Take5). Le traitement a été effectué avec la même chaîne, le prototype de MACCS, développé et maintenu au CESBIO par Mireille Huc. Notre méthode de base est une méthode multi-temporelle à la fois pour la détection des nuages, des ombres de nuages, de l'eau et pour l'estimation de l'épaisseur optique des aérosols.

 

Cependant, grâce à la richesse spectrale de LANDSAT, nous avons pu enrichir nos méthodes par rapport à la version utilisée pour SPOT4 (Take5) : grâce à la bande bleue, nous pouvons utiliser un critère complémentaire pour détecter les aérosols, grâce à la relation quasi constante observée entre les réflectances des bandes bleues et rouges au dessus de la végétation. Le gain de précision dû à la présence de cette bande permet de compenser la perte de précision de la méthode multi-temporelle due à la faible répétitivité des LANDSAT. Enfin, contrairement à LANDSAT 8, la bande 1.38 n'existe pas sur LANDSAT 5 et 7, la détection des nuages hauts n'est donc pas évidente.

Images de la tuile obtenue sur la côte atlantique, en provenance de deux traces différentes de LANDSAT (à gauche au milieu, les traces 201 et 200 ). Les angles de visée sont légèrement différents sur chacune des traces (visée depuis l'ouest sur l'image de gauche visée depuis l'est sur l'image de droite). A droite une image LANDSAT 7 de la trace 201, réduite à la portion centrale.

 

Pour augmenter la précision de la détection des nuages, nous avons décidé d'utiliser les données issues de traces adjacentes de LANDSAT (5 ou 7) dans les séries temporelles de niveau 2A. Ces données ne sont pas acquises exactement sous le même angle (+/- 7 degrés), mais la différence d'angle est suffisamment petite pour qu'il y ait un vrai gain de précision sur les zones d'intersection entre traces. En raison de cette approximation, quelques artefacts peuvent être observés.

 

Format des données

Le format des produits de Niveau 2A de LANDSAT 5 et 7 est le même que celui des données SPOT4 (Take5).

 

Défauts connus :

Exemple du phénomène de rémanence observé sur les données de LANDSAT 5 à proximité d’un gros nuage. Ce défaut électronique se traduit ici par une alternance de bandes sombres et claires au dessus de la végétation.

Voici une petite liste des défauts connus des données LANDSAT 5 ou 7 :

 

  • la donnée de référence pour l'ortho-rectification à l'USGS peut présenter des biais supérieurs à 30 mètres (38 mètres à Toulouse). Les données LANDSAT 7 peuvent donc parfois être légèrement décalées par rapport aux données LANDSAT 5 ortho-rectifiées au CNES.
  • L'instrument TM sur LANDSAT 5 présente un phénomène de rémanence qui se traduit par des bandes plus ou moins sombres, perpendiculaires à la trace du satellite, a proximité de zones contrastées, comme par exemple près d'une importante masse nuageuse.
  • Exemple de parasites pouvant apparaître sur les images LANDSAT 5, sous la forme de points colorés.

     

  • Les produits LANDSAT 5 que l'ESA nous a fournis voilà deux ans présentent parfois des "parasites" sous la forme de points colorés apparaissant de manière aléatoire sur l'une ou l'autre bande.
  • Dans les produits LANDSAT, la valeur qui indique si un pixel est en dehors de l'image est égale à 0. Or cette valeur peut aussi se retrouver à l'intérieur des données. Nous avons essayé de séparer les vraies valeurs hors image des valeurs normales dans l'image mais nous n'y arrivons pas toujours, et dans ce cas, toutes les valeurs de toutes les bandes sont mise à la valeur nodata, qui chez nous vaut -10 000 justement pour éviter ces problèmes.

 

Le traitement des données acquises par LANDSAT8 sur la France en 2013 est terminé.

=>

Ça torche* ! J'annonçais le démarrage du traitement dans un précédent article: il aura fallu moins de 15 jours pour traiter l'ensemble des données acquises par LANDSAT 8 en 2013 sur la France métropolitaine (y compris la corse). La période traitée s'étend d'avril 2013 à fin 2013. Les traitements ont été effectués par l'équipe MUSCATE, pour le compte du pôle THEIA, en utilisant les moyens du centre informatique du CNES.

 

Il suffira de quelques jours de plus pour mettre les données en ligne sur le site du pôle et vérifier que les données peuvent être distribuées. Finalement, le plus long dans ce traitement est le téléchargement des produits de Niveau 1T (l'équivalent du Niveau 1C dans la nomenclature de THEIA) sur le site earthexplorer de l'USGS.

 

Les données obtenues sont de bonne qualité, comme le montrent les exemples de quicklooks de niveau 2A, à droite, tirés de la série temporelle obtenue sur la tuile de la région parisienne. Comme d'habitude, les nuages détectés sont entourés en vert, les ombres en noir, la neige en rose et l'eau en bleu. Quelques nuages sont parfois oubliés par notre méthode multi-temporelle de détection nuageuse, lorsque la répétitivité des données est insuffisante, comme, sur l'image ci contre, lors de ce printemps 2013 extrêmement nuageux. Les images suivantes de la série ne présentent pas ce genre d'erreurs.

 

Cet article décrit brièvement les principales étapes du traitement qui ont permis d'obtenir ce résultat.

 

Point de départ : Niveau 1T

Nous téléchargeons les données depuis le site earthexplorer, en utilisant une version améliorée du script de téléchargement décrit ici. Ces données sont orthorectifiées par l'USGS, à partir d'une base de données mondiale de points d'appuis.

 

Les spécifications de localisation pour les données LANDSAT de l'USGS sont de 50 mètres. Nous avons constaté, sur la région de Toulouse, des erreurs de localisation égales à 1,5 pixels, mais d'autres régions sont mieux loties. L'USGS a confirmé les erreurs que nous observons du côté de Toulouse (38m vers le sud) et promet de les corriger. Notre traitement ne corrige pas ces erreurs qui restent modérées, et pour obtenir une performance améliorée, il faudra attendre le retraitement de ces données par l'USGS, ou rattraper le décalage éventuel par vous mêmes.

 

En ce qui concerne l'étalonnage absolu du satellite LANDSAT 8, nous utilisons celui recommandé par la NASA (et fourni avec les produits LANDSAT 8).

 

Reprojection en Lambert 93


Les données de l'USGS sont fournies en projection UTM. Cette projection utilise trois fuseaux différents au dessus de la France, qui se divise donc en 3 zones différentes, l'Ouest de la France, le Centre et l'Est. Comme les données de deux fuseaux différents ne se superposent pas directement, nous avons donc décidé de reprojeter les données en Lambert 93, qui est la projection officielle pour la France.

Découpage des produits en tuiles

Nous avons pris le parti de suivre la même logique que celle utilisée par Sentinel-2, et de découper les données en tuiles de 110*110 km décalées de 100 km les unes par rapport aux autres. La tuile 1x1 se trouve au sud-ouest de la France, lorsqu'on va vers l'est (vers la droite), on incrémente la première coordonnée de D0001 à D0010 (D pour Droite), lorsqu'on va vers le Nord (vers le haut), on incrémente la seconde coordonnée, de H0001 à H0010 (H pour Haut). La tuile de Toulouse s'appelle donc D0005H0002  Le découpage en tuiles est visible sur l'image ci-jointe.

 

Pour la Corse, nous avons défini une autre grille de tuiles, composée de deux tuiles.

 

Pour chacune des tuiles, nous fournissons l'ensemble des dates pour lesquelles une image LANDSAT8 a une intersection non nulle avec la tuile. Quelques dates peuvent manquer, pour plusieurs raisons, liées en général à la couverture nuageuse :

  • l'image n'a pas été acquise par LANDSAT 8 (quand les prévisions météo indiquent un temps très couvert, les images ne sont pas acquises).
  • l'image a été acquise mais s'est avérée trop nuageuse pour être traitée au niveau 1T par l'USGS
  • l'image est trop nuageuse pour être traitée par la chaîne de Niveau 2A

 

Traitement de Niveau 2A (correction atmosphériques et détection des nuages)

Il est important de noter que notre chaîne ne traite pas les bandes thermiques pour l'instant. Une correction est à l'étude, mais celle-ci ne sera pas opérationnelle avant un ou deux ans.

 

La méthode utilisée pour les bandes visible, proche et moyen infra rouge est quasiment la même que pour SPOT4(Take5). Le traitement a été effectué avec la même chaîne, le prototype de MACCS, développé et maintenu au CESBIO par Mireille Huc. Notre méthode de base est une méthode multi-temporelle à la fois pour la détection des nuages, des ombres de nuages, de l'eau et pour l'estimation de l'épaisseur optique des aérosols.

 

Cependant, grâce à la richesse spectrale de LANDSAT, nous avons pu enrichir nos méthodes par rapport à la version utilisée pour SPOT4 (Take5) :  la bande 1.38µm de LANDSAT8 permet de détecter les nuages hauts. Et grâce à la bande bleue, nous pouvons utiliser un critère complémentaire pour détecter les aérosols, grâce à la relation quasi constante observée entre les réflectances des bandes bleues et rouges au dessus de la végétation. Le gain de précision dû à la présence de cette bande permet de compenser la perte de précision de la méthode multi-temporelle due à la faible répétitivité de LANDSAT8.

Images de la tuile obtenue sur la région parisienne, en provenance de trois traces différentes de LANDSAT (de gauche à droite, les traces 200, 199 et 198). Les angles de visée sont légèrement différents sur chacune des traces (visée depuis l'ouest sur l'image de gauche, visée au nadir au centre, visée depuis l'est sur l'image de droite).

 

Pour augmenter la précision de la détection des nuages, nous avons décidé d'utiliser les données issues de traces adjacentes de LANDSAT8 dans les séries temporelles de niveau 2A. Ces données ne sont pas acquises exactement sous le même angle (+/- 7 degrés), mais la différence d'angle est suffisamment petite pour qu'il y ait un vrai gain de précision sur les zones d'intersection entre traces. En raison de cette approximation, quelques artefacts peuvent être observés.

 

Pour gagner davantage de répétitivité, nous pourrions aussi entrelacer des données LANDSAT 7 et LANDSAT 8 dans des séries temporelles communes, mais ceci ne sera implémenté que dans une prochaine version.

 

Format des données

Le format des produits de Niveau 2A de LANDSAT 8 est le même que celui des données SPOT4 (Take5), à quelques détails près que je fournirai dans une nouvelle page, dès que j'ai un moment...

 

* Ca va vite, oui, c'est comme ça qu'on parle chez les anciens jeunes du CNES...

 

 

Production of LANDSAT L2A data at THEIA to begin shortly.

At CNES, the prototype MUSCATE production facility of THEIA land data centre will soon start the production and distribution of Level 2A Landsat 5 and 7 data, and shortly after of LANDSAT 8 data covering the entire surface of France.

 

Mosaic of LANDSAT 5 & 7 data produced at CESBIO, from both ESA and USGS data. These data are cut in 110 x 110 km² tiles. For each tile, each LANDSAT acquisition with at least a little clear sky corner is provided.

 

For Landsat 5 and 7, we use data from both USGS and ESA : indeed, up to now only ESA has the LANDSAT 5 data that were acquired over the receiving stations of Mas Palomas (Canary Islands), Matera (Italy), and Svalbard (Norway). A transfer to USGS of ESA's data is expected, it may have started in Svalbard, but it has not yet begun for the Matera station, which covers France.

Level 1C

The USGS data are orthorectified, but those from ESA are not, so, as for SPOT4 (Take5), we set up an ortho-rectification processing using the SIGMA tool provided by CNES.  The ESA's products we received 2 years ago also have some flaws (which may have been corrected by now, but given it took months to obtain the data we did not ask for a reprocessing): the thermal band is unusable and you will find here and there colourful bright spots, such as those produced by your neighbour moped on your TV when you were a child. Nevertheless, we can produce correct Level 1 products, although we look forward to the reprocessing of these products by USGS. ESA has now it own processing of LANDSAT data, but it stops at level 1C.

For Landsat 7, this processing is not necessary because the data are already ortho-rectified. We interpolate only a small portion of the missing data due to LANDSAT 7 SLC breakdown, and then we discard the parts of the image where the gaps are too large. For LANDSAT 8, none of these processings are needed.

 

Level 2A

The Level 2A products (Cloud Mask, Atmospheric correction) will be produced by the prototype of MACCS software developed and maintained by Mireille Huc (CESBIO, CNRS). Two years ago, I had already produced such a data set on the most Southern part of France, at CESBIO. These products are already distributed on THEIA web site and are also used to illustrate this post.

 

LANDSAT 5 and 7 :

Starting in April, we will process the LANDSAT 5 & 7 data acquired above France from 2009 to 2011.

 

LANDSAT 8 :

From April or May, we will process the LANDSAT 8 data acquired since april 2013, and we will try to keep the pace so as to produce the incoming new LANDSAT 8 acquisitions with a short delay.

 

Data Format :

We will reuse the data format of SPOT4 (Take5). France will be split into 110*110 km² tiles with a 10 km overlap with their neighbours. (See the image mosaic obtained for the South of France).

Depending on the success of the distribution of these data, we will decide if it is worth producing older time periods or other regions. Please tell us if you need such data.

Exemple : available LANDSAT images from July to October 2009 for the tile centered on Toulouse. For each date, we provide the level 1C image tTOA reflectance), on the left, and the level 2A image on the right (surface reflectance). The detected clouds are circled in red.

SMAC prend en charge de nouveaux satellites

=>

De nouveaux coefficients ont été ajoutés au site du CESBIO . Les nouveaux satellites pris en compte sont :

  • Landsat 8, RapidEye
  • Quickbird, Worldview2, Ikonos
  • Pléiades1A (PHR1A)
  • DMC-DEIMOS1

 

Le Simplifié Modèle d'Atmosphérique Correction (SMAC) est parfaitement adapté à l'implémentation rapide et approchée de corrections atmosphériques. Il s'agit de fonctions analytiques dérivées du modèle 5S. Les 49 coefficients de ce modèle sont ajustés à partir de simulations de transfert radiatif obtenues avec le modèle 6S (l'ancienne version, pas la récente version vectorielle). SMAC n'est pas un modèle très précis (beaucoup moins que MACCS), et il faut lui fournir des données auxiliaires pour l'épaisseur optique des aérosols ou pour les contenus atmosphériques en ozone et vapeur d'eau. Quand ces données sont précisément connues, la précision des simulations est en général meilleure que deux à trois pour cent, sauf parfois pour les grands angles (au dessus de 70°) ou dans de fortes bandes d'absorption et si on ne prend pas en compte les effets d'environnement et les effets de pente.

 

SMAC a été conçu pour être facile à utiliser :

#lecture des 49 cofficients
nom_smac ='COEFS/coef_FORMOSAT2_B1_CONT.dat'
coefs=coeff(nom_smac)
 
#Lire la réflectance TOA de la bande à traiter
#(celà va dépendre du format de l'image)
 
#Lire les valeurs des angles dans les métadonnées (ici, on les fixe)
theta_s=30
phi_s=180
theta_v=0
phi_v=0
# calculer la pression atmosphérique à l'altitude du pixel
pressure=PdeZ(1300)
 
#trouver les valeurs des variables atmsophériques (ici, on les fixe arbitrairement)
AOT550=0.1
UO3=0.3
UH2O=3
 
#calculer la correction atmosphérique
r_surf=smac_inv(r_toa,theta_s,phi_s,theta_v,phi_v,pressure,AOT,UO3,UH2O,coefs)

Dans la dernière ligne ci-dessus :

  • theta_s, phi_s sont resp. l'angle solaire zenithal et azimuthal
  • theta_v, phi_v sont resp. l'angle de visée zenithal et azimutha
  • AOT est l'épaisseur optique à 550 nm qui peut provenir d'une station Aeronet, ou fixée au jugé, ou égale à 0.1 (si on se contente d'une correction très approchée).
  • UO3 est le contenu en ozone, en cm.atm (0.3 est souvent suffisamment précis)
  • UH2O est le contenu en vapeur d'eau, en kg/m². J'utilise souvent une valeur égale à 3, quand je ne cherche pas une grande précision, dans des bandes ou l'absorption est faible
References

[1] Rahman, H., & Dedieu, G. (1994). SMAC: a simplified method for the atmospheric correction of satellite measurements in the solar spectrum. REMOTE SENSING, 15(1), 123-143.
"[2]"Tanré, D., Deroo, C., Duhaut, P., Herman, M., Morcrette, J. J., Perbos, J., & Deschamps, P. Y. (1990). Technical note Description of a computer code to simulate the satellite signal in the solar spectrum: the 5S code. International Journal of Remote Sensing, 11(4), 659-668.

"[3]"Vermote, E. F., Tanré, D., Deuze, J. L., Herman, M., & Morcette, J. J. (1997). Second simulation of the satellite signal in the solar spectrum, 6S: An overview. Geoscience and Remote Sensing, IEEE Transactions on, 35(3), 675-686.>
"[4]"Kotchenova, S. Y., Vermote, E. F., Matarrese, R., & Klemm Jr, F. J. (2006). Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance. Applied Optics, 45(26), 6762-6774.
"[5]"Kotchenova, S. Y., & Vermote, E. F. (2007). Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part II. Homogeneous Lambertian and anisotropic surfaces. Applied Optics, 46(20), 4455-4464.

 

Validation des épaiseurs optiques d'aérosols obtenues pour SPOT4(Take5)

=>

Nous sommes en train de préparer un retraitement des données SPOT4 (Take5) pour la fin de l'année. Pour en choisir les paramètres, j'ai compilé tous les résultats de validation des épaisseurs optiques mesurées par notre méthode MACCS, en les comparant avec celles mesurées par le réseau Aeronet. Dans la méthode MACCS appliquée à SPOT4(Take5), en l'absence d'une bande spectrale dans le bleu, l'estimation de l'épaisseur optique peut être faite de deux manières :

  • soit par inversion de l'épaisseur optique par la méthode multi-temporelle,
  • soit par bouchage de trous, les trous pouvant être dus à la présence de nuages, d'eau ou de neige, ou à la présence de pixels de réflectance trop forte ou ayant varié trop fortement entre deux dates.

 

Comparaison des épaisseur optiques mesurées avec MACCS avec celles mesurées in-situ par des photomètres du réseau Aéronet. Les points bleus correspondent aux cas où l'atmosphère est stable et où il n'y a pas de nuages dans le voisinage. Les points rouges correspondent aux cas où l'atmosphère est moins stable (variation de l'épaisseur optique en une heure) ou aux cas où des nuages sont détectés dans le voisinage de la mesure (20*20 km).
La courbe de gauche ne retient que les dates et sites ayant moins de 60% de pixels dont l'épaisseur optique est obtenue par bouchage de trous et pour celle de droite, seules les dates et sites ayant moins de 20% de pixels bouchés sont retenus. On constate que le bouchage de trous n'introduit pas trop d'erreurs dans les cas stables, mais en introduit davantage dans les cas instables.

 

Les estimations d'aérosols ont été obtenues avec le prototype de MACCS développé et maintenu par Mireille Huc au CESBIO. Le modèle d'aérosols n'est pas celui qui a été utilisé pour la première production de SPOT4 (Take5). Il s'agit d'aérosols un peu plus gros (rayon modal 0.2 µm, contre 0.1µm dans le traitement initial), qui fournissent un meilleur accord global avec les mesures d'Aeronet. C'est ce modèle qui sera utilisé dans le futur retraitement des données SPOT4 (Take5), pour la plupart des sites.

 

L'écart-type des mesures d'épaisseur optique pour les cas stables est de 0.06, ce qui est du niveau de l'état de l'art, et constitue une performance remarquable, compte tenu de la difficulté à estimer l'épaisseur optique des aérosols sans bande bleue. De plus, certains des sites Aéronet (Bruxelles, Gwangju, Ouarzazate, Wallops, NASA_LaRC) utilisés sont assez éloignés de l'image SPOT, parfois de plus de 60 kilomètres.

Les sites aéronet utilisés ici sont :

Site SPOT4 Take5 Site Aeronet
Belgique Bruxelles
South Great Plains Cart_Site
Korea Gwangju_GIST
Chesapeake NASA_LaRC
Chesapeake Wallops
Versailles Paris
Versailles Palaiseau
Tunisia Ben Salem
Maroc Saada
Maroc Ouarzazate
Sudmipy-Est Seysses + Le Fauga
Sudmipy-Ouest Seysses
Provence Carpentras
Provence Frioul

 
Les plus mauvais résultats sont obtenus pour les sites :

  • Gwangju (Korée): le site SPOT est en bord de mer, alors que le site Aeronet est à l'intérieur des terres, dans une grande ville (notamment les points rouges avec une forte épaisseur optique sur la courbe "_60".)
  • Ben Salem (Tunisie): sur ce site très nuageux au début de la période, de fortes variations des réflectances de surface alors que les images à peu près claires sont très espacées dans le temps.
  • Palaiseau et Paris : dans ce cas, le modèle d'aérosols utilisé pour tous les sites n'a pas l'air d'être le bon modèle, il faudrait rajouter des aérosols carbonés et absorbants.

En revanche, les sites Maroc, Provence (même l'extrapolation à l'ile du Frioul), Sudmipy, Wallops et Cart_site donnent d'excellents résultats.

 

Enfin, certains utilisateurs nous on rapporté des problèmes de correction atmosphérique pour les les sites forestiers équatoriaux, mais malheureusement (ou heureusement pour nos statistiques), aucun d'entre eux ne dispose d'un instrument du réseau AERONET.

 

Les effets d'environnement, comment ça marche ?

=>

Comme expliqué dans l'article sur les effets atmosphériques, la diffusion de la lumière par les molécules et les aérosols présents dans l'atmosphère provoque plusieurs effets. La diffusion ajoute un voile aux données (la réflectance atmosphérique), atténue le signal en provenance de la surface (la transmission atmosphérique), et rend les images floues (les effets d'environnement). Cet article s'intéresse aux effets d'environnement, les autres aspects ont été abordés dans le lien fourni ci-dessus.

 

Le schéma ci-joint montre les différents types de trajets que peut suivre la lumière avant d'arriver au capteur. Le Trajet 1 correspond à la réflectance atmosphérique, le trajet 2 est proportionnel à la réflectance de la cible observée atténué par sa traversée de l'atmosphère, c'est celui qui nous intéresse et nous permet de retrouver la réflectance de surface. Les trajets 3 et 4 apportent au capteur une part de signal qui ne provient pas directement de la surface que le satellite observe mais de son voisinage (d'où le nom d'"effets d'environnement"). Ce sont ces trajets qui apportent du flou sur l'image.

 

Les effets d'environnement peuvent engendrer de fortes erreurs lorsqu'on observe une parcelle de végétation entourée de sols nus ou de végétation senescente. Pour un tel cas, la figure ci-dessous présente les erreurs en pourcentage de réflectance et de NDVI, si on ne prend pas en compte les effets d'environnement, en fonction de l'épaisseur optique des aérosols.

 

Il est possible de corriger ces effets de manière approchée, à condition de connaître la quantité d'aérosols. Dans les traitements de la chaîne MACCS, nous procédons de la manière suivante :

 

  1. Nous procédons à la correction atmosphérique en supposant que le paysage est uniforme. Nous obtenons une réflectance de surface sous hypothèse uniforme que nous notons  \rho_{s,unif} .
  2. Nous calculons la réflectance de l'environnement du pixel (  \rho_{s,env} ) en utilisant un filtre de convolution gaussien de 2 km de diamètre, qui calcule une moyenne pondérée de la réflectance environnante. En toute rigueur, ce filtre devrait dépendre de la quantité d'aérosols présents dans l'atmosphère (moins il y a d'aérosols, plus le rayon devrait être grand), et des angles de prise de vue, mais nous n'avons pas encore travaillé sur cet aspect, nous avons donc utilisé un filtre constant.
  3. Nous corrigeons finalement la réflectance des trajets 3 et 4 par la formule suivante :

 \rho_{s}=\frac{\rho_{s,unif}.T^{\uparrow}.\frac{1-\rho_{s,unif}.s}{1-\rho_{s,env}.s}-\rho_{s,env}. T_{dif}^{\uparrow}}{T_{dir}^{\uparrow}}

  •  T^{\uparrow}=T_{dif}^{\uparrow}+T_{dir}^{\uparrow} est la transmission atmosphérique montante totale, somme de la transmission atmosphérique diffuse et directe. s est l'albedo atmosphérique. Toutes ces grandeurs sont déduites de calculs de transfert radiatif et dépendent de la quantité et du type d'aérosols.

 

Cette correction qui implique l'utilisation de convolutions est assez lourde et prend près d'un quart du temps de correction atmosphérique.

 

Exemples de résultats

Les images ci dessous présentent 3 stades de la correction atmosphérique pour deux images Formosat-2 acquises au dessus du Canada, à deux jours d'intervalle, la première image est acquise un jour il y a beaucoup d'aérosols (épaisseur optique de 0.47 d'après nos calculs), alors que la seconde est acquise un jour très clair (épaisseur optique de 0.1 selon nos calculs).

 

  • La première ligne correspond aux images au sommet de l'atmosphère, sans correction atmosphérique. On voit bien que l'image de gaucheest plus floue.
  • La deuxième ligne correspond aux images corrigées en supposant le paysage uniforme. Il s'agit de l'image obtenue à l'issue de l'étape 1 dans la méthode décrite ci-dessus. L'image de gauche est toujours plus floue.
  • La troisième ligne présente ces mêmes images après la correction d'environnement. Dans ce cas, l'image de gauche n'est plus floue, elles est même légèrement trop nette (un peu de sur correction).

Images TOA (à gauche, l'image avec fort contenu en aérosols)

 

Images en réflectance de surface, en supposant le paysage uniforme (à gauche, l'image avec fort contenu en aérosols)

 

Images en réflectance de surface après correction des effets d'environnement (à gauche, l'image avec fort contenu en aérosols)

 

On peut aussi comparer point à point les réflectances pour juger de l'amélioration après correction des effets d'environnement. La courbe ci-dessous compare les images corrigées en supposant le paysage uniforme, et les images corrigées en tenant compte des effets d'environnement. On constate que les points se rapprochent de la diagonale après correction des effets d'environnement. Sur l'image du 27 mai, pour laquelle l'épaisseur optique est la plus forte, on note que les fortes réflectances sont un peu trop faibles, alors que les faibles réflectances sont un peu trop fortes, ce qui correspond bien à une perte de contraste.

L'USGS distribue des produits LANDSAT corrigés des effets atmosphériques

=>

L'USGS et de la NASA viennent de mettre en place la production et la distribution de produits LANDSAT de Niveau 2A (exprimés en réflectance de surface, corrigés des effets atmosphériques, et accompagnés d'un masque des nuages et de leurs ombres, et de masques de neige et d'eau). Pour le moment, les données thermiques sont exprimées en température de brillance au sommet de l'atmosphère. Ce produit est déjà disponible pour LANDSAT 5 et LANDSAT 7 et la NASA travaille actuellement à l'adaptation de ses chaînes à LANDSAT -8.

 

Depuis la semaine dernière, ces données sont disponibles sur le site earthexplorer, il faut chercher un peu pour les trouver: cliquer sur l'onglet "data sets", et choisir "Landsat CDR". Vous pouvez donc commander vos données, qui sont prêtes dès le lendemain (sauf peut-être si vous en commandez des centaines). Il est seulement dommage que l'USGS ne dispose pas d'une grande partie des données LANDSAT historiques (Landsat 4 et 5) acquises dans les stations de l'ESA. Un transfert de ces données est prévu, mais il n'a pas encore eu lieu (nous l'attendons avec impatience).

 

Comparaison des pourcentages de nuages détectés par MACCS (=MTCD) et LEDAPS sur des scènes LANDSAT aux Etats-Unis. L'accord est très bon dans la plupart des cas, même si on note dans quelques cas rares une sous détection de nuages par LEDAPS. L'un des points les plus éloignés de l'axe est illustré ci-dessous.

Les produits de réflectance de surface LANDSAT ont été obtenus avec la chaîne LEDAPS de [Vermote et Masek 2006]. La chaîne LEDAPS utilise la méthode DDV pour estimer l'épaisseur optique des aérosols, et une batterie de tests portant sur une seule date pour détecter les nuages. Nous avions fait l'an dernier une comparaison avec nos méthodes (MACCS), et les résultats étaient assez proches, même si la méthode mono date ne peut pas tout à fait rivaliser avec la précision des méthodes multi-temporelles (cf figures ci-dessous). LEDAPS est également un peu moins bon pour la détection des ombres, encore plus difficiles à détecter que les nuages. Mais au total, l'existence de ce produit va certainement simplifier le travail de beaucoup d'utilisateurs.

 

Le produit LANDSAT Surface Reflectance est pour le moment le seul produit de niveau 2A à haute résolution existant ou prévu, avec une couverture mondiale.  Cependant, le pôle THEIA a produit des données de niveau 2A à partir des données de l'expérience SPOT4 (Take5) (disponibles ici) et à partir  des données acquises par LANDSAT 5, 7 ou 8 sur la France depuis 2009 (disponibles là).

Selon les plans actuels, l'ESA ne fournira pas de produit 2A pour Sentinel-2. Il est dores et déjà prévu que THEIA fasse tourner la chaîne MACCS pour fournir des données de niveau 2A pour Sentinel-2 sur toute la France, et très probablement sur quelques autres pays ou régions à définir, couvrant au total une surface équivalente à 10 fois la France. L'application de MACCS au monde entier n'est pas encore décidée, même si THEIA se positionne pour le faire.

 

Comparaison des masques de nuages MACCS (gauche) et LEDAPS (Droite), dans un cas difficile avec des nuages très fins, mieux détectés (entourés en rouge) par MACCS. Dans les cas standards, les masques des deux méthodes sont plus proches.

 

Reference : Masek, J. G., Vermote, E. F., Saleous, N. E., Wolfe, R., Hall, F. G., Huemmrich, K. F., ... & Lim, T. K. (2006). A Landsat surface reflectance dataset for North America, 1990-2000. Geoscience and Remote Sensing Letters, IEEE, 3(1), 68-72.