Les noms des niveaux de produits, comment ça marche ?

=>

Simulations des produits Sentinel-2 à partir d'images Formosat-2

 

Vous êtes nombreux à ne pas retenir les noms des produits d'observation de la terre, et si je n'étais pas tombé dedans quand j'étais petit, je ne ferais pas mieux.  Peut-être qu'en détaillant la logique qui a conduit à leur définition, il sera plus facile de les retenir. Voici donc comment les noms des produits du pôle THEIA ont été définis :

 

  • nous avons décidé d'utiliser la nomenclature de Sentinel-2 : 1C, 2A,  3A, car nous sommes convaincus que Sentinel-2 deviendra rapidement la référence en imagerie spatiale
    • Le produit de Niveau 1C est une image monodate ortho-rectifiée, exprimée en réflectance au sommet de l'atmosphère.
    • Le produit de Niveau 2A est une image monodate ortho-rectifiée, exprimée en réflectance de surface, accompagnée d'un masque de nuages/ombres de nuages/neige/eau
    • le produit de Niveau 3A est une synthèse mensuelle ortho-rectifiée des pixels non nuageux provenant des Niveau 2A.
  • cette nomenclature, définie par l'ESA et le CNES respecte les normes du Comittee on Earth Observation Satellites (CEOS)

 

Les normes du CEOS sont assez difficiles à trouver, je les avais cherchées à plusieurs reprises. Heureusement, j'ai fini par trouver quelques membres du comité travaillant sur l'interopérabilité des systèmes, et deux d'entre eux (Frédéric Baret (INRA), et Kenneth McDonald (NOAA)) ont répondu immédiatement à mes questions.

 

La norme du CEOS reprend en fait les termes d'une norme de la NASA, datant de 1986, qui elle est disponible sur Wikipedia :

 

Data Level NASA-EOSDIS Definition

http://science.nasa.gov/earth-science/earth-science-data/data-processing-levels-for-eosdis-data-products/
Level 0 Reconstructed, unprocessed instrument and payload data at full resolution, with any and all communications artifacts (e.g., synchronization frames, communications headers, duplicate data) removed. (In most cases, the EOS Data and Operations System (EDOS) provides these data to the data centers as production data sets for processing by the Science Data Processing Segment (SDPS) or by a SIPS to produce higher-level products.)
Level 1A Reconstructed, unprocessed instrument data at full resolution, time-referenced, and annotated with ancillary information, including radiometric and geometric calibration coefficients and georeferencing parameters (e.g., platform ephemeris) computed and appended but not applied to Level 0 data.
Level 1B Level 1A data that have been processed to sensor units (not all instruments have Level 1B source data).
Level 2 Derived geophysical variables at the same resolution and location as Level 1 source data.
Level 3 Variables mapped on uniform space-time grid scales, usually with some completeness and consistency.
Level 4 Model output or results from analyses of lower-level data (e.g., variables derived from multiple measurements).

 

Je ne sais pas ce que vous en pensez, mais je trouve cette description assez floue, la notion de "sensor units" n'est pas claire, et surtout, il faut attendre le produit de niveau 3 pour que les données soient rééchantillonnées sur une grille cartographique. Si l'utilisateur ne veut pas avoir à se soucier de rééchantillonnage géométrique, il doit utiliser le produit de Niveau 3.

 

La norme du CEOS (ci dessous), clairement inspirée de celle de la NASA, permet de rééchantillonner les données dès le Niveau 1 et de les exprimer en grandeur physiques et non pas pas en "sensor units". La norme du CEOS ne détaille pas les sous niveaux (1A, 1B). Il me semble que cette norme ne détaille cependant pas bien la différence entre les grandeurs physiques et variables géophysiques (par exemple, pourquoi une réflectance au sommet de l'atmosphère serait -elle une grandeur physique et une réflectance de surface une variable biophysique), et laisse beaucoup de marge d’interprétation.

 

Data Level CEOS Definition 

http://www.ceos.org/images/WGISS/Documents/Handbook.pdf

Level 0 Reconstructed unprocessed instrument data at full space time resolution with all available supplemental information to be used in subsequent processing (e.g., ephemeris, health and safety) appended.
Level 1 Unpacked, reformatted level 0 data, with all supplemental information to be used in subsequent processing appended. Optional radiometric and geometric correction applied to produce parameters in physical units. Data generally presented as full time/space resolution. A wide variety of sub level products are possible.
Level 2 Retrieved environmental variables (e.g., ocean wave height, soil moisture, ice concentration) at the same resolution and location as the level 1 source data.
Level 3 Data or retrieved environmental variables which have been spatially and/or temporally re-sampled (i.e., derived from level 1 or 2 products). Such re-sampling may include averaging and compositing.
Level 4 Model output or results from analyses of lower level data (i.e., variables that are not directly measured by the instruments, but are derived from these measurements).

Ayant eu bien du mal à trouver ces fameuses normes du CEOS, je m'en étais fait ma propre idée, probablement à partir de discussions avec les collègues qui travaillaient avec moi sur le projet POLDER. Voici ma version personnelle de la norme (sachant qu'une norme personnelle ne sert à rien bien sûr...)

 

Data Level La nomenclature des produits selon OH
Level 0 Données reconstruites mais non traitées, fournies à pleine résolution, accompagnées de toutes les informations nécessaires pour les traitements de niveaux supérieurs.
Level 1 Le produit de Niveau 1 est un produit dont tous les pixels ont été acquis à la même date (en quelques instants, au cours d'un passage du satellite), et dont le traitement ne fait pas d'hypothèses sur la nature du paysage observé. Chaque pixel de l'image est traité de la même manière, quelle que soit sa nature (nuage, forêt, mer...). Le produit exprime les données en grandeur physique ou fournit toutes les informations nécessaires pour les convertir. Il peut être rééchantillonné sur une grille cartographique ou doit fournir toutes les informations nécessaires à cette conversion.
Level 2 Le produit de Niveau 2 doit être lui aussi acquis en quelques instants, mais cette fois, des traitements différenciés par type de pixel sont autorisés, et l'on peut faire des hypothèses sur l'atmosphère, l'état de la végétation ou de la mer, et traiter les différentes classes de pixels d'une manière différente.
Level 3 Le produit de Niveau 3 est constitué de données acquises à des dates différentes, souvent avec des emprises différentes, et, comme sur le Niveau 2, les traitements peuvent faire des hypothèses sur ce qui est observé. Les produits de niveau 3 sont en général des synthèses périodiques (hebdomadaires, décadaires, mensuelles) des produits de Niveau 2.

 

 

La norme du CEOS ne définit les sous-niveaux de produits, mais dans ce domaine, la norme de la NASA a une certaine influence. Pour Sentinel-2, les produits de niveau 1A et 1B sont des produits internes assez proche de la norme de la NASA. mais le produit de niveau 1C a fournit des données ortho-rectifiées dans une projection cartographique. Seul le niveau 1C devrait être distribué aux utilisateurs normaux.

 

 

La spécification de mission de Sentinel-2 définit également un produit de Niveau 2A, exprimé en réflectance de surface après correction atmosphérique, et un produit de niveau 2B fournissant des variables biophysiques (LAI, FAPAR...). Aucun de ces produits ne sera cependant fourni par le segment sol de Sentinel-2. Ces produits existeront en revanche dans le centre de traitement de THEIA, et nous les avons déja produits pour SPOT4 (Take5).

 

Enfin, notons que de nombreux produits d'observation de la terre ne respectent pas ces normes, souvent parce que la nomenclature a été définie avant l'adoption de la norme : c'est le cas de SPOT et de ses successeurs, y compris les plus récents (Pleiades, Rapid Eye).

 

The terrain effect correction : how it works

=>


Caution, this post contains formulas.

 

The topographic (or terrain) effects on the observed reflectances are due to several phenomena, illustrated below :

  • the closer the surface is perpendicular to sun direction, the more energy it receives per surface unit (we talk about irradiance). If the surface is parallel to sun direction, it does not receive direct sunlight. We can model it this way :
    • For an horizontal surface :  E_h= E_0.T_{dir}^\downarrow.cos(\theta_s)

      Definition of sun zenith angle and sun incidence angle

    • For a sloped surface  E_i=E_0.T_{dir}^\downarrow.cos(\theta_i)
    •  E_0 is the Top of Atmosphere irradiance, and  T_{dir}^\downarrow is the downward direct transmission, i.e. the proportion of the light that reaches directly the surface without being scattered by the atmosphere.
    • Assuming that all the irradiance is direct, the measured reflectance if the surface was horizontal is calculated from the following formula:  \rho_h=\rho_i \frac{cos(\theta_s)}{cos(\theta_i)} . However, the above assumption is not true and this formula tends to over correct terrain effects
  • The surfaces also receive a diffuse sun irradiance scattered by the atmosphere. If the surface is not horizontal, a part of the sky is obscured by the slope reducing the diffuse irradiance. Moreover, the diffuse irradiance depends on the amount of aerosols (and clouds) in the atmosphere. In addition, the surrounding terrain can also hide a part of the sky, but we do not take this effect into account here in our modelling. We use the following approximation, which is equivalent to assuming that the slope is alone in a horizontal region.
    • If surface is horizontal, the visible sky fraction is 1, if it is vertical, this fraction is 1/2
    •  \displaystyle F_{sky}= \frac{1+cos(slope)}{2}
  • Finally, the slope can receive light from surrounding surfaces, which become directly visible. In our simplified model, we always assume the entire environment of our slope is flat and,for instance, we do not take the effect of the opposite side in a valley into account :
    • If surface is horizontal, the visible ground fraction is 0, if it is vertical, it is 1/2.
    •  \displaystyle F_{fround}= \frac{1-cos(slope)}{2}

 

Finally, we use the following formula  to compute the reflectance that would be observed if the surface was horizontal   \rho_{h} , as a function of the slope (inclined) reflectance   \rho_{i} :

  \rho_{h}=\displaystyle \rho_{i}.\frac{T^{\downarrow}}{T_{dir}^{\downarrow}.\frac{cos(\theta_i)}{cos(\theta_s)} + T_{dif}^{\downarrow} F_{sky} + T^{\downarrow} F_{ground} \rho_{env}}

 T^{\downarrow} is the downward transmission, sum of direct and diffuse irradiances :  T^{\downarrow}= T_{dir}^{\downarrow}+ T_{dif}^{\downarrow} , and \rho_{env} is the average reflectance of the neighbourhood.

 

Finally, we can also account for bidirectional reflectance effects, but this correction is tricky since directional effects depend on the surface cover type. See for instance : Dymond, J.R.; Shepherd, J.D. 1999: Correction of the topographic effect in remote sensing. IEEE Trans. Geosci. Remote Sens. 37(5): 2618-2620.

 

It is very difficult to validate the correction of directional effects : , we could compare the correction results for satellite overpasses at different times in the day. But all the satellite optical imagers have nearly the same overpass time. A qualitative way of estimating the accuracy is to check that similar land covers on opposite slopes in a valley ( a meadow, a forest) have a similar reflectance after correction. The most suitable points are North-South valleys. Here are some examples of terrain effects correction results.

Formosat-2 image in the Alps, before (left) and after (right) terrain effect correction. The image on the right seems to have been flattened, and the opposite slope reflectances seem much more alike after correction.

Finally, an essential part of the method's accuracy is the availability of a highly accurate digital  elevation model (DEM),  up to now, only the SRTM DEM is available globally, and it only has a 90 meter resolution. Its accuracy is somewhat inadequate and sometimes leaves artefacts if the slope changes are poorly located.

 

 

La correction des variations d'éclairement dues au relief

=>


Attention, cet article contient des formules.

 

Les effets du relief sur les réflectances observées proviennent de plusieurs phénomènes, schématisés dans l'illustration ci-dessous.

  • plus la surface est perpendiculaire à la direction solaire, plus elle va recevoir d'énergie par unité de surface (on parle d'éclairement). Si la surface est parallèle à la direction solaire, elle ne reçoit pas de lumière directe du soleil, l'éclairement direct reçu est donc nul.
    • Pour une surface horizontale :  E_h= E_0.T_{dir}^\downarrow.cos(\theta_s)

      Définition de l'angle zenithal solaire et de l'angle d'incidence

    • Pour une surface inclinée  E_i=E_0.T_{dir}^\downarrow.cos(\theta_i)
    •  E_0 est l'éclairement solaire en haut de l'atmosphère, et  T_{dir}^\downarrow est la transmission directe descendante, c'est à dire la proportion de cet éclairement qui rejoint directement la surface sans être diffusé par l'atmosphère.
    • Si on suppose que tout l'éclairement est direct, la réflectance mesurée si la surface avait été horizontale se calcule à partir de la formule suivante :  \rho_h=\rho_i \frac{cos(\theta_s)}{cos(\theta_i)} . Cette formulation a tendance à surcorriger l'effet des pentes.
  • Les surfaces reçoivent aussi un éclairement solaire diffusé par l'atmosphère. Si la surface n'est pas horizontale, une partie du ciel est masquée par la pente: l'éclairement diffusé par l'atmosphère et qui parvient à la surface est donc réduit. Cet éclairement va varier en fonction de la quantité d'aérosols (et de nuages) présents dans l'atmosphère. Par ailleurs, les reliefs environnants peuvent aussi masquer une partie du ciel, mais nous ne prenons pas cet effet là en compte dans notre modélisation. Nous utilisons l'approximation suivante, qui revient à supposer que la pente est le seul endroit incliné dans un environnement horizontal.
    • Si la surface est horizontale, la fraction de ciel Fciel vue est 1, si la surface est verticale, la fraction de ciel vue est 1/2
    •  \displaystyle Fciel= \frac{1+cos(pente)}{2}
  • enfin, la pente peut recevoir de la lumière des surfaces environnantes, qui deviennent directement visibles. Dans notre modélisation simplifiée, nous supposons toujours que tout l'environnement de notre pente est plat et nous ne prenons pas en compte l'effet du versant opposé dans une vallée.
    • Si la surface est horizontale, la fraction de sol vue est 0, si la surface est verticale, la fraction de sol vue est 1/2.
    •  \displaystyle Fsol= \frac{1-cos(pente)}{2}

 

Au final, pour calculer la réflectance mesurée si la surface avait été horizontale   \rho_{h} , en fonction de la réflectance mesurée sur la surface inclinée   \rho_{i} , on obtient donc la formule suivante :

  \rho_{h}=\displaystyle \rho_{i}.\frac{T^{\downarrow}}{T_{dir}^{\downarrow}.\frac{cos(\theta_i)}{cos(\theta_s)} + T_{dif}^{\downarrow} F_{ciel} + T^{\downarrow} F_{sol} \rho_{env}}

 T^{\downarrow} est la transmission descendante composée d'une partie directe et diffuse :  T^{\downarrow}= T_{dir}^{\downarrow}+ T_{dif}^{\downarrow} , et \rho_{env} est la réflectance moyenne autour du pixel à corriger.

 

Enfin, on peut aussi prendre en compte un facteur de correction lié aux effets directionnels, mais cet aspect est délicat, car les effets directionnels varient en fonction du type de surface observé. Voir par exemple  : Dymond, J.R.; Shepherd, J.D. 1999: Correction of the topographic effect in remote sensing. IEEE Trans. Geosci. Remote Sens. 37(5): 2618-2620.

 

Il est très difficile de valider cette correction : pour bien le faire, il faudrait comparer des images acquises à différentes heures de la journée. Malheureusement, tous les satellites à haute résolution passent à peu près à la même heure. On peut essayer cependant de vérifier que des surfaces ayant un couvert semblable (une prairie, une forêt) sur des versants opposés, présentent bien la même réflectance après correction. Les points les plus propices à ce genre de vérifications sont des vallées Nord Sud. Voici quelques exemples de résultats après correction.

Image Formosat-2 acquise dans les alpes, avant (à gauche) et après (à droite) correction des effets du relief. L'image de droite donne l'impression d'avoir été aplatie, les versants de part et d'autre des vallées présentent des réflectances bien plus semblables après correction.

Enfin, un élément essentiel de la précision de cette méthode est la disponibilité d'un modèle numérique de terrain très précis, pour le moment, le seul MNT disponible globalement est SRTM à 90 mètres de résolution. Sa précision est un peu insuffisante et laisse parfois des artefacts, si les ruptures de pentes sont mal localisées.

 

The adjacency effects, how they work.

As explained in the post about atmospheric effects, the scattering of light by molecules and aerosols in the atmosphere brings about several effects : scattering adds some haze on the images (the atmospheric reflectance), lessens the signal from the surface (the atmospheric transmission), and blurs the images (the adjacency effects). This post is about the adjacency effects, the other aspects have already been quickly explained in the above post.

 

The figure on the right shows the types of paths that light can follow before getting to the satellite. Path 1 corresponds to the atmospheric reflectance, path 2 is path that interacts with the target, it is the one which is useful to determine the surface reflectance, paths 3 and 4 contribute to the total reflectance but interact with the surface away from the target. These paths are thus the cause of adjacency effects and they blur the images.

 

 

If not corrected, adjacency effects may cause large errors. Let's take the case of a fully developed irrigated field surrounded by bare soil. For such a case, the second figure on the right shows the relative percentage of errors for reflectances and NDVI as a function of aerosol optical thickness, if adjacency effect is not corrected.

 

 

 

An approximate correction can be applied, but it thus requires to know the aerosol optical thickness. In our MACCS processor, here is how it works :

 

  1. We first correct the images under the assumption that the Landscape is uniform. We obtain a surface reflectance under uniform absorption which is noted  \rho_{s,unif} .
  2. We compute the neighbourhood reflectance (  \rho_{s,adj} ) using a convolution filter with a 2km radius, that computes the average neighborhood reflectance weighted by the distance to the target. To be fully rigorous, this filter should depend on the optical thickness and on the viewing and sun angle (The less aerosols, the larger radius), but as we did not work on an accurate model, we used a constant radius.
  3. We correct for the contribution of paths 3 and 4 using :

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

  • where  T^{\uparrow}=T_{dif}^{\uparrow}+T_{dir}^{\uparrow} is the total upward transmission, sum of diffuse and direct upward transmissions, and s is the atmosphere spheric albedo. These quantities depend on the wavelength, on the aerosol model and on the AOT. They are computed using Look up Tables based on radiative transfer calculations.

 

As this processing uses convolution with a large radius, it takes quite a large part of the atmospheric processing time.

 

Result Exemples

The images below show 3 stages of the atmospheric processing, for 2 Formosat-2 images obtained over Montreal (Canada) with a 2 days interval. The first image was acquired on a hazy day (aerosol optical thickness (AOT) of 0.47 according to MACCS estimate); and the second one on a clear day (AOT=0.1).

  • The first line corresponds to the Top Of Atmosphere images, without atmospheric correction. The left image is obviously blurred compared to the right image.
  • The second line corresponds to the atmospheric correction under uniform landscape assumption (as in step 1). The left image is still obviously blurred compared to the right image.
  • the third line show the same images after adjacency effect correction. In that case, the left image is not blurred any more, it is even maybe a little over corrected as it seems somewhat sharper that the right image.

TOA Images (On the left, the hazy image)


Surface reflectance under uniform landscape assumption (on the left, the hazy image)

 

Surface reflectance after adjacency effect correction (on the left, the hazy image)

 

The pixel wise comparison of reflectances is also a way to show the enhancement due to the adjacency effect correction. The plot below compares the images of both dates corrected under the uniform landscape assumption (on the left), and after adjacency effect correction (on the right). You may observe that the dots are closer the the black diagonal on the right. On the hazy image (May 27th), the high reflectances are a little too low, while the low reflectances are a little too high, which is the symptom of a loss of contrast.

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.

How to estimate Aerosol Optical Thickness

=>

Caution ! This post contains formulas !


Aerosols play a great role in the atmospheric effects. Aerosols are particles suspended in the atmosphere, which can be of several types: sand or dust, soot from combustion, sulfates or sea salt, surrounded by water... Their size ranges between 0.1 micron and a few microns, depending on the type of aerosol or on the air moisture. Their quantity is also extremely variable : rain can suddenly reduce their abundance (known as "aerosol optical thickness"). The abundance variations result in great variations of observable reflectances from one day to the next, and it is therefore necessary to know the quantity and type of aerosols, in order to correct their effects.

 

Unfortunately, to correct the effects of aerosols, there is no global aerosol observation network, and the only available data are local observations from the few hundred points of Aeronet network. Therefore, this network can not be used operationally to correct the satellite images over large areas.

Weather forecast models just start predicting the amounts of aerosols, based on satellite observations and modeling of sources and sinks and of the transport of aerosols by the winds, but these data do not seem to have sufficient accuracy yet to be used for the atmospheric correction of images.

 

Our atmospheric correction method, named MACCS, is therefore based on an estimate of aerosol optical depth from the images themselves. To understand how this method works, one must already understand the effects of aerosols on radiation. We have seen in this post, that the effects of diffusion can be modelled as follows (assuming the corrected gas absorption):

ρTOA = ρatm +Td ρsurf

The reflectance at the top of the atmosphere ρTOA (Top of Atmosphere) is the sum of the atmospheric reflectance  ρatm and of the surface reflectance ρsurf transmitted by the atmosphere. We seek to know the surface reflectance, but for each measurement made at the top of the atmosphere, there are three unknowns to be determined. To separate the effects of the atmosphere and surface effects, we must use other information.

 

Dark pixel method

When the image includes a surface whose surface reflectance is nearly zero, the reflectance observed at the top of the atmosphere becomes ρTOA = ρatm. We can therefore deduce the atmospheric reflectance and using a radiative transfer model, the aerosols optical thickness (AOT). Finally, knowing the AOT, we can compute the diffuse transmission, and finally calculate ρsurf. An even simpler and more approximate version of this method consists in subtracting directly the reflectance of the dark pixel (or ρatm) to the entire image (neglecting the transmission) [Chavez, 1988].

 

However, this method assumes that there is a very dark area in the image (which is not always the case), and that the reflectance of the dark surface is known. The method also assumes that the amount of aerosols is constant over the image and it neglects the effect of terrain. The results obtained by this method can be quite inaccurate. In our method (MACCS), however, we use the method of black pixel determine the maximum value of the optical thickness in the area.

 

Multi Spectral Method, called "DDV"

If you know the type of aerosols in the atmosphere, it is possible to deduce the properties of aerosols in a spectral band from the optical properties in another spectral band.

 

If there are two spectral bands, there are two measures ρsurf and three unknowns (both surface reflectance in these bands, and the amount of aerosols). An additional equation can be obtained if we know the relationship between the surface reflectance of the two bands.

 

The method named "Dark Dense Vegetation" (DDV) is based on assumptions about relationships between surface reflectances of the dense vegetation exploiting the fact that the spectrum of dense green vegetation is quite constant. The most famous version of this method is that used by NASA for MODIS project [Remer 2005]. It connects the surface reflectance in the blue and red with those in the SWIR. This provides two equations for estimating the type of aerosol optical thickness. This method works well in temperate and boreal zones, but not in arid areas where it is difficult to find the dense vegetation. Early versions used the following equations:

 

ρBlue = 0.5 ∗ ρSWIR

ρRed = 0.25 ∗ ρSWIR

 

The following versions of the MODIS DDV algorithm are a bit more complicated but follow the same principle. Our work has shown that using the equation below allows a more accurate determination of the optical thickness, for less dense vegetation cover (NDVI to a 0.2) because bare soil brown also respect this relationship.

 

ρBlue = 0.5 * ρRed

(the exact value of the coefficient is adjusted according to the spectral bands of the instrument)


This version of the  method, however, does not allow to determine the aerosol model. In the case of SPOT4 (Take5), the absence of a blue band does not allow us to use this equation, resulting in a slight loss in accuracy.

This diagram shows that the correlation between surface reflectance above vegetation is much better for the (blue, red) couple of spectral bands than for couples including using (SWIR).

 

 

 

Multi Temporal Method

In most cases, the reflectance of the land surface changes slowly over time, while the aerosol optical properties vary rapidly from one day to another. We can therefore consider what changes from one image to another (apart from special cases often linked to human intervention) is associated with aerosols, and deduce the properties of aerosols and then correct for atmospheric effects. This method is too complex to be explained in detail here, interested readers can refer to [Hagolle 2008].

 

So that surface reflectance be nearly constant from one image to another, however, it is required that images be acquired at a constant angle. Indeed, the reflectance depend on the viewing angles : this is what we call directional effects. This method therefore applies only to satellite observations obtained with constant angle. It does not apply to standard SPOT data, but this condition is true for SPOT4 (Take5) data. It will also apply to Landsat Venμs and Sentinel-2.

 

Finally :

 

Validation of aerosol optical thickness (AOT) from time-series of FORMOSAT-2 images, depending on the method (multi-spectral, multi-temporal, combined), compared with the measurements provided by the Aeronet network of in-situ measurements. The multi-spectral method works best on sites covered with vegetation and is much less accurate on arid sites, while the multi-temporal method performs a little worse on green sites, but much better on dry sites. The combination of the two methods retains the best of the two basic methods.

The MACCS method, used for SPOT4 (Take5) experiment, and also for LANDSAT, VENμS and Sentinel-2 data, combines the three methods described above to obtain robust estimates of aerosol optical thickness. These methods work in many cases, but sometimes fail when the assumptions on which they are based prove to be incorrect. They generally tend to work better on vegetated areas rather than in arid areas. for now, they assume the model known aerosol and in the coming years, we will look for reliable ways to identify the type of aerosols.

 

References :
Chavez Jr, P. S. (1988). An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment, 24(3), 459-479.

Remer, L. A., and Coauthors, 2005: The modis aerosol algorithm, products, and validation. J. Atmos. Sci., 62, 947–973.

Hagolle, O and co-authors, 2008. « Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images ». Remote sensing of environment 112 (4)

Hagolle, O.; Huc, M.; Villa Pascual, D.; Dedieu, G. A Multi-Temporal and Multi-Spectral Method to Estimate Aerosol Optical Thickness over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENμS and Sentinel-2 Images. Remote Sens. 2015, 7, 2668-2691.

Les aérosols jouent un rôle prépondérant dans les effets atmosphériques. Les aérosols sont des particules en suspension dans l'atmosphère, qui peuvent être de plusieurs types : grains de sable ou poussières, suies issues de combustion, sulfates ou sels marins entourés d'eau... Leur taille peut varier de 0.1 µm à quelques microns, en fonction du type d'aérosols ou de l'humidité de l'air. Quant à leur quantité, elle est extrêmement variable, une pluie pouvant réduire brutalement leur abondance (on parle d'"épaisseur optique d'aérosols"). Ils peuvent faire varier fortement d'un jour à l'autre les réflectances observables depuis le sommet de l'atmosphère et il est donc nécessaire de connaître leur quantité et leur type afin de pouvoir corriger leurs effets.

 

Malheureusement, pour corriger les effets des aérosols, on ne dispose pas de réseau global d'observation des aérosols, seulement d'observations locales, sur les quelques centaines de points du réseau Aeronet. Ce réseau ne peut donc pas être utilisé pour corriger opérationnellement les images de satellites sur de grandes étendues.
Des modèles météorologiques commencent à prédire les quantités d'aérosols, en se basant sur les observations de satellites et la modélisation des sources et du transport des aérosols par les vents, mais ces données ne semblent pas encore avoir une précision suffisante pour être utilisées pour la correction atmosphérique des images.

 

Notre méthode de correction atmosphérique (MACCS) repose donc sur une estimation de l'épaisseur optique des aérosols à partir des images elles-mêmes. Pour bien comprendre le fonctionnement de cette méthode, il faut déjà comprendre les effets des aérosols sur le rayonnement. On a vu, dans ce billet, que les effets de la diffusion peuvent être modélisés ainsi (on suppose l'absorption gazeuse corrigée) :

ρTOA = ρatm +Td ρsurf

La réflectance au sommet de l'atmosphère ρTOA (Top of Atmosphere) est la somme de la réflectance atmosphérique ρatm et de la réflectance de surface ρsurf transmise par l'atmosphère. On cherche à connaître la réflectance de surface, mais à chaque mesure réalisée au sommet de l'atmosphère, on a trois inconnues à déterminer. Pour séparer les effets de l'atmosphère et les effets de la surface, il faut donc utiliser d'autres informations.

 

Méthode du pixel noir

Lorsque l'image contient une surface dont la réflectance est quasi nulle, la réflectance observée au sommet de l'atmosphère devient ρTOA= ρatm. On peut donc en déduire la réflectance atmosphérique, et en utilisant un modèle de transfert radiatif, l'épaisseur optique des d'aérosols. On peut enfin en déduire la transmission diffuse, et finalement calculer ρsurf. Une version encore plus simple et plus approximative consiste à soustraire directement la réflectance du pixel sombre (soit ρatm) à toute l'image. [Chavez, 1988]

 

Cependant, cette méthode revient à supposer qu'il existe bien une surface très sombre dans l'image (ce qui n'est pas toujours le cas), et que la réflectance de cette surface sombre est connue. La méthode suppose aussi que la quantité d'aérosols est constante dans l'image et elle néglige les effets du relief. Les résultats obtenus par cette méthode peuvent donc être assez imprécis. Dans notre méthode (MACCS), nous utilisons cependant la méthode du pixel noir déterminer la valeur maximale de l'épaisseur optique dans la zone.

 

Méthode Multi Spectrale, dite "DDV"

Si on connaît le type d'aérosols présent dans l'atmosphère, il est possible de déduire les  propriétés des aérosols dans une bande spectrale, à partir des propriétés optiques dans une autre bande spectrale.

 

Si on dispose de deux bandes spectrales, on dispose de deux mesures ρsurf et de trois inconnues( les deux réflectances de surface dans ces bandes, et la quantité d'aérosols). Une équation supplémentaire peut être obtenue si on connaît la relation entre les réflectances de surface des deux bandes.

 

La méthode  méthode "Dark Dense Vegetation" (DDV ) est basée sur des hypothèses de relations entre réflectances de surface sur la végétation dense exploitant le fait que le spectre de la végétation dense et verte est un peu toujours le même. La version la plus connue de cette méthode est celle utilisée par la NASA pour le projet MODIS [Remer 2005]. Elle relie les réflectances de surface dans le bleu et dans le rouge avec celles dans le moyen infra-rouge. On dispose ainsi de deux équations qui permettent d’estimer le type d’aérosols et l’épaisseur optique. Cette méthode fonctionne bien en zones tempérées et boréales, mais pas en zones arides, où il est difficile de trouver de la végétation dense. Les premières versions utilisaient les équations suivante

ρBleu = 0.5 ∗ ρSWIR

ρRouge = 0.25 ∗ ρSWIR

Les versions suivantes ont un peu compliqué ces équations, sans en modifier le principe. Nos travaux ont montré que l’utilisation de l'équation ci dessous  (la valeur exacte du coefficient est à ajuster en fonction des bandes spectrales de l'instrument):

ρBleu = 0.5 ∗ ρRouge

permet une détermination plus précise de l’épaisseur optique, pour des couverts végétaux moins denses (jusqu’à un NDVI de 0.2), car les sols nus de couleur marron respectent aussi cette relation. La méthode ne permet pas, par contre, de déterminer le modèle d’aérosols. Dans le cas de SPOT4 (Take5) l'absence d'une bande bleue ne nous permet pas d'utiliser cette dernière équation, d’où une légère perte en précision.

Ce diagramme montre que la corrélation entre réflectances de surface au dessus de la végétation est bien meilleure pour le couple de bandes spectrales (bleu, rouge) que pour les couples incluant le moyen infra rouge. (SWIR)

 

Méthode Multi Temporelle

On observe dans la plupart des cas que les réflectances de la surface terrestre évoluent lentement avec le temps, alors que le propriétés optiques des aérosols varient très rapidement, d'un jour à l'autre. On peut donc considérer que ce qui change d'une image à l'autre (en dehors de cas particuliers souvent liées à des interventions humaines) est lié aux aérosols, et donc en déduire les propriétés des aérosols pour ensuite corriger les effets atmosphériques. Cette méthode est un peu trop complexe pour être expliquée en détails ici, les lecteurs intéressés pourront se reporter à [Hagolle 2008].

 

Pour que les réflectances de surface soient quasi constantes d'une image à l'autre, il faut cependant que les images soient acquises sous un angle de vue constant. Les changements d'angles d'observation font en effet varier les réflectances (ce phénomène sera prochainement expliqué dans un autre article). Cette méthode ne s'applique donc qu'aux seuls satellites permettant des observations à angle constant.  Elle ne s'applique donc pas aux données SPOT normales mais par contre convient parfaitement aux données SPOT4 (Take5). Elle s'appliquera aussi à Landsat, Venµs et Sentinel-2.

En résumé :

Performance de l'estimation de l'épaisseur optique des aérosols sur des séries temporelles d'images Formosat-2,, en fonction de la méthode (multi-spectrale, multi-temporelle, combinée), par comparaison avec les mesures fournies par le réseau de mesures in-situ Aeronet. La méthode multi spectrale fonctionne mieux sur des sites couverts de végétation et moins bien sur des sites arides, la méthode multi-temporelle marche un peu moins bien sur les sites verts, mais beaucoup mieux sur les sites arides. La combinaison des deux méthodes garde le meilleur des deux méthodes élémentaires.

 

Notre méthode MACCS, utilisée pour l'expérience SPOT4 (Take5), et pour les données LANDSAT, VENµS et Sentinel-2, combine les trois méthodes présentées ci-dessus pour obtenir des estimations robustes des épaisseurs optiques d'aérosols. Ces méthodes fonctionnent dans un grand nombre de cas, mais peuvent parfois échouer quand les hypothèses sur lesquelles elles reposent s'avèrent fausses. Elles ont en général tendance à mieux fonctionner sur des zones couvertes de végétation plutôt que dans des zones arides. pour le moment, elles supposent le modèle d'aérosol connu, et dans les prochaines années, nous chercherons des manières fiables d'identifier le type d'aérosols.

 

References :
Chavez Jr, P. S. (1988). An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment, 24(3), 459-479.

Remer, L. A., and Coauthors, 2005: The modis aerosol algorithm, products, and validation. J. Atmos. Sci., 62, 947–973.

Hagolle, O and co-authors, 2008. « Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images ». Remote sensing of environment 112 (4)

L'estimation du contenu atmosphérique en aérosols

=>

Attention, cet article contient des formules !

 

Les aérosols jouent un rôle prépondérant dans les effets atmosphériques. Les aérosols sont des particules en suspension dans l'atmosphère, qui peuvent être de plusieurs types : grains de sable ou poussières, suies issues de combustion, sulfates ou sels marins entourés d'eau... Leur taille peut varier de 0.1 µm à quelques microns, en fonction du type d'aérosols ou de l'humidité de l'air. Quant à leur quantité, elle est extrêmement variable, une pluie pouvant réduire brutalement leur abondance (on parle d'"épaisseur optique d'aérosols"). Ils peuvent faire varier fortement d'un jour à l'autre les réflectances observables depuis le sommet de l'atmosphère et il est donc nécessaire de connaître leur quantité et leur type afin de pouvoir corriger leurs effets.

 

Malheureusement, pour corriger les effets des aérosols, on ne dispose pas de réseau global d'observation des aérosols, seulement d'observations locales, sur les quelques centaines de points du réseau Aeronet. Ce réseau ne peut donc pas être utilisé pour corriger opérationnellement les images de satellites sur de grandes étendues.
Des modèles météorologiques commencent à prédire les quantités d'aérosols, en se basant sur les observations de satellites et la modélisation des sources et du transport des aérosols par les vents, mais ces données ne semblent pas encore avoir une précision suffisante pour être utilisées pour la correction atmosphérique des images.

 

Notre méthode de correction atmosphérique (MACCS) repose donc sur une estimation de l'épaisseur optique des aérosols à partir des images elles-mêmes. Pour bien comprendre le fonctionnement de cette méthode, il faut déjà comprendre les effets des aérosols sur le rayonnement. On a vu, dans ce billet, que les effets de la diffusion peuvent être modélisés ainsi (on suppose l'absorption gazeuse corrigée) :

ρTOA = ρatm +Td ρsurf

La réflectance au sommet de l'atmosphère ρTOA (Top of Atmosphere) est la somme de la réflectance atmosphérique ρatm et de la réflectance de surface ρsurf transmise par l'atmosphère. On cherche à connaître la réflectance de surface, mais à chaque mesure réalisée au sommet de l'atmosphère, on a trois inconnues à déterminer. Pour séparer les effets de l'atmosphère et les effets de la surface, il faut donc utiliser d'autres informations.

 

Méthode du pixel noir

Lorsque l'image contient une surface dont la réflectance est quasi nulle, la réflectance observée au sommet de l'atmosphère devient ρTOA= ρatm. On peut donc en déduire la réflectance atmosphérique, et en utilisant un modèle de transfert radiatif, l'épaisseur optique des d'aérosols. On peut enfin en déduire la transmission diffuse, et finalement calculer ρsurf. Une version encore plus simple et plus approximative consiste à soustraire directement la réflectance du pixel sombre (soit ρatm) à toute l'image. [Chavez, 1988]

 

Cependant, cette méthode revient à supposer qu'il existe bien une surface très sombre dans l'image (ce qui n'est pas toujours le cas), et que la réflectance de cette surface sombre est connue. La méthode suppose aussi que la quantité d'aérosols est constante dans l'image et elle néglige les effets du relief. Les résultats obtenus par cette méthode peuvent donc être assez imprécis. Dans notre méthode (MACCS), nous utilisons cependant la méthode du pixel noir déterminer la valeur maximale de l'épaisseur optique dans la zone.

 

Méthode Multi Spectrale, dite "DDV"

Si on connaît le type d'aérosols présent dans l'atmosphère, il est possible de déduire les  propriétés des aérosols dans une bande spectrale, à partir des propriétés optiques dans une autre bande spectrale.

 

Si on dispose de deux bandes spectrales, on dispose de deux mesures ρsurf et de trois inconnues( les deux réflectances de surface dans ces bandes, et la quantité d'aérosols). Une équation supplémentaire peut être obtenue si on connaît la relation entre les réflectances de surface des deux bandes.

 

La méthode  méthode "Dark Dense Vegetation" (DDV ) est basée sur des hypothèses de relations entre réflectances de surface sur la végétation dense exploitant le fait que le spectre de la végétation dense et verte est un peu toujours le même. La version la plus connue de cette méthode est celle utilisée par la NASA pour le projet MODIS [Remer 2005]. Elle relie les réflectances de surface dans le bleu et dans le rouge avec celles dans le moyen infra-rouge. On dispose ainsi de deux équations qui permettent d’estimer le type d’aérosols et l’épaisseur optique. Cette méthode fonctionne bien en zones tempérées et boréales, mais pas en zones arides, où il est difficile de trouver de la végétation dense. Les premières versions utilisaient les équations suivante :

 

ρBleu = 0.5 ∗ ρSWIR

ρRouge = 0.25 ∗ ρSWIR

 

Les versions suivantes ont un peu compliqué ces équations, sans en modifier le principe. Nos travaux ont montré que l’utilisation de l'équation ci dessous  (la valeur exacte du coefficient est à ajuster en fonction des bandes spectrales de l'instrument):

ρBleu = 0.5 ∗ ρRouge

 

permet une détermination plus précise de l’épaisseur optique, pour des couverts végétaux moins denses (jusqu’à un NDVI de 0.2), car les sols nus de couleur marron respectent aussi cette relation. La méthode ne permet pas, par contre, de déterminer le modèle d’aérosols. Dans le cas de SPOT4 (Take5) l'absence d'une bande bleue ne nous permet pas d'utiliser cette dernière équation, d’où une légère perte en précision.

Ce diagramme montre que la corrélation entre réflectances de surface au dessus de la végétation est bien meilleure pour le couple de bandes spectrales (bleu, rouge) que pour les couples incluant le moyen infra rouge. (SWIR)

 

Méthode Multi Temporelle

On observe dans la plupart des cas que les réflectances de la surface terrestre évoluent lentement avec le temps, alors que le propriétés optiques des aérosols varient très rapidement, d'un jour à l'autre. On peut donc considérer que ce qui change d'une image à l'autre (en dehors de cas particuliers souvent liées à des interventions humaines) est lié aux aérosols, et donc en déduire les propriétés des aérosols pour ensuite corriger les effets atmosphériques. Cette méthode est un peu trop complexe pour être expliquée en détails ici, les lecteurs intéressés pourront se reporter à [Hagolle 2008].

 

Pour que les réflectances de surface soient quasi constantes d'une image à l'autre, il faut cependant que les images soient acquises sous un angle de vue constant. Les changements d'angles d'observation font en effet varier les réflectances : c'est ce qu'on appelle les effets directionnels. Cette méthode ne s'applique donc qu'aux seuls satellites permettant des observations à angle constant.  Elle ne s'applique donc pas aux données SPOT normales mais par contre convient parfaitement aux données SPOT4 (Take5). Elle s'appliquera aussi à Landsat, Venµs et Sentinel-2.

 

En résumé :
Performance de l'estimation de l'épaisseur optique des aérosols sur des séries temporelles d'images Formosat-2,, en fonction de la méthode (multi-spectrale, multi-temporelle, combinée), par comparaison avec les mesures fournies par le réseau de mesures in-situ Aeronet. La méthode multi spectrale fonctionne mieux sur des sites couverts de végétation et moins bien sur des sites arides, la méthode multi-temporelle marche un peu moins bien sur les sites verts, mais beaucoup mieux sur les sites arides. La combinaison des deux méthodes garde le meilleur des deux méthodes élémentaires.

 

Notre méthode MACCS, utilisée pour l'expérience SPOT4 (Take5), et pour les données LANDSAT, VENµS et Sentinel-2, combine les trois méthodes présentées ci-dessus pour obtenir des estimations robustes des épaisseurs optiques d'aérosols. Ces méthodes fonctionnent dans un grand nombre de cas, mais peuvent parfois échouer quand les hypothèses sur lesquelles elles reposent s'avèrent fausses. Elles ont en général tendance à mieux fonctionner sur des zones couvertes de végétation plutôt que dans des zones arides. pour le moment, elles supposent le modèle d'aérosol connu, et dans les prochaines années, nous chercherons des manières fiables d'identifier le type d'aérosols.

 

References :
Chavez Jr, P. S. (1988). An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment, 24(3), 459-479.

Remer, L. A., and Coauthors, 2005: The modis aerosol algorithm, products, and validation. J. Atmos. Sci., 62, 947–973.

Hagolle, O and co-authors, 2008. « Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images ». Remote sensing of environment 112 (4)

Hagolle, O.; Huc, M.; Villa Pascual, D.; Dedieu, G. A Multi-Temporal and Multi-Spectral Method to Estimate Aerosol Optical Thickness over Land, for the Atmospheric Correction of FormoSat-2, LandSat, VENμS and Sentinel-2 Images. Remote Sens. 2015, 7, 2668-2691.

Land cover map production: how it works

=>

Land cover and land use maps

Although different, the terms land use and land cover are often used as synonymous. From Wikipedia Land cover is the physical material at the surface of the earth. Land covers include grass, asphalt, trees, bare ground, water, etc. There are two primary methods for capturing information on land cover: field survey and analysis of remotely sensed imagery. and Land use is the human use of land. Land use involves the management and modification of natural environment or wilderness into built environment such as fields, pastures, and settlements. It also has been defined as "the arrangements, activities and inputs people undertake in a certain land cover type to produce, change or maintain it" (FAO, 1997a; FAO/UNEP, 1999).


A precise knowledge of land use and land cover is crucial for many scientific studies and for many operational applications. This accurate knowledge needs frequent information updates, but may also need to be able to go back in time in order to perform trend analysis and to suggest evolution scenarios.

 

Satellite remote sensing offers the possibility to have a global point of view over large regions with frequent updates, and therefore it is a very valuable tool for land cover map production.

 

However, for those maps to be available in a timely manner and with a good quality, robust, reliable and automatic methods are needed for the exploitation of the available data.

 

 

 

Classical production approaches

The automatic approaches to land cover map production using remote sensing imagery are often based on image classification methods.

 

This classification can be:

  • supervised: areas for which the land cover is known are used as learning examples;
  • unsupervised: the image pixels are grouped by similarity and the classes are identified afterwards.

Supervised classification often yields better results, but it needs reference data which are difficult or costly to obtain (field campaigns, photo-interpretation, etc.).

 

 

 

What time series bring

Until recently, fine scale land cover maps have been nearly exclusively produced using a small number of acquisition dates due to the fact that dense image time series were not available.

 

The focus was therefore on the use of spectral richness in order to distinguish the different land cover classes. However, this approach is not able to differentiate classes which may have a similar spectral signature at the acquisition time, but that would have a different spectral behaviour at another point in time (bare soils which will become different crops, for instance). In order to overcome this problem, several acquisition dates can be used, but this needs a specific date selection depending on the map nomenclature.

 

For instance, in the left image, which is acquired in May, it is very difficult to tell where the rapeseed fields are since they are very similar to the wheat ones. On the right image, acquired in April, blooming rapeseed fields are very easy to spot.

 

May image. Light green fields are winter crops, mainly wheat and rapeseed. But which are the rapeseed ones?

April image. Blooming rapeseed fields are easily distinguished in yellow while wheat is in dark green.

 

If one wants to build generic (independent from the geographic sites and therefore also from the target nomenclatures) and operational systems, regular and frequent image acquisitions have to be ensured. This will soon be made possible by the Sentinel-2 mission, and it is right now already the case with demonstration data provided by Formosat-2 and SPOT4 (Take 5). Furthermore, it can be shown that having a high temporal resolution is more interesting than having a high spectral diversity. For instance, the following figure shows the classification performance results (in terms of  \kappa index, the higher the better) as a function of the number of images used. Formosat-2 images (4 spectral bands) and simulated Sentinel-2 (13 bands) and Venµs (12 bands) data have been used. It can be seen that, once enough acquisitions are available, the spectral richness is caught up by a fine description of the temporal evolution.

kappaVFS.png

 

 

What we can expect from Sentinel-2

Sentinel-2 has unique capabilities in the Earth observation systems landscape:

  • 290 km. swath;
  • 10 to 60 m. spatial resolution depending on the bands;
  • 5-day revisit cycle with 2 satellites;
  • 13 spectral bands.

Systems with similar spatial resolution (SPOT or Landsat) have longer revisit periods and fewer and larger spectral bands. Systems with similar temporal revisit have either a lower spatial resolution (MODIS) or narrower swaths (Formosat-2).

 

The kind of data provided by Sentinel-2 allows to foresee the development of land cover map production systems which should be able to update the information monthly at a global scale. The temporal dimension will allow to distinguish classes whose spectral signatures are very similar during long periods of the year. The increased spatial resolution will make possible to work with smaller minimum mapping units.

 

However, the operational implementation of such systems will require a particular attention to the validation procedures of the produced maps and also to the huge data volumes. Indeed, the land cover maps will have to be validated at the regional or even at the global scale. Also, since the reference data (i.e. ground truth) will be only available in limited amounts, supervised methods will have to be avoided as much as possible. One possibility consists of integrating prior knowledge (about the physics of the observed processes, or via expert rules) into the processing chains.

 

Last but not least, even if the acquisition capabilities of these new systems will be increased, there will always be temporal and spatial data holes (clouds, for instance). Processing chains will have to be robust to this kind of artefacts.

 

 

Ongoing work at CESBIO

 

Danielle Ducrot, Antoine Masse and a few CESBIO interns have recently produced a a large land cover map over the Pyrenees using 30 m. resolution multi-temporal Landsat images. This map, which is real craftsmanship, contains 70 different classes. It is made of 3 different parts using nearly cloud-free images acquired in 2010.

 

70-class land cover map obtained from multi-temporal Landsat data.

In his PhD work, Antoine works on methods allowing to select the best dates in order to perform a classification. At the same time, Isabel Rodes is looking into techniques enabling the use of all available acquisitions over very large areas by dealing with both missing data (clouds, shadows) and the fact that all pixels are not acquired at the same dates.

 

These 2 approaches are complementary: one allows to target very detailed nomenclatures, but needs some human intervention, and the other is fully automatic, but less ambitious in terms of nomenclature.

 

A third approach is being investigated at CESBIO in the PhD work of Julien Osman: the use of prior knowledge both quantitative (from historical records) and qualitative (expert knowledge) in order to guide the automatic classification systems.

 

We will give you more detailed information about all those approaches in coming posts on this blog.

La production de cartes d'occupation du sol, comment ça marche?

=>

Les cartes d'occupation du sol

D'après Wikipédia, l'occupation du sol désigne pour la FAO (1998) "la couverture (bio-)physique de la surface des terres émergées" et donc le type d'usage (ou de non-usage) fait des terres par l'Homme. La mosaïque paysagère est cartographiée en identifiant les types homogènes de milieux (ex : zones artificialisées, zones agricoles, forêts ou landes, zones humides, etc.).


La connaissance précise de cette occupation du sol est un enjeu crucial pour beaucoup de travaux de recherche et pour de nombreuses applications opérationnelles. Une connaissance précise demande une mise à jour fréquente de ces informations, mais peut aussi nécessiter de remonter dans le temps pour faire une analyse des tendances et proposer des scénarios d'évolution.

 

La possibilité offerte par la télédétection spatiale d'accéder à une vue d'ensemble de grandes régions de façon récurrente constitue donc un atout majeur pour la production de cartes d'occupation du sol.

 

Cependant, pour que ces cartes soient disponibles dans des délais raisonnables et avec une qualité suffisante, il est nécessaire de disposer de méthodes automatiques robustes et fiables, capables d'exploiter de façon efficace les données disponibles.

 

 

Les approches classiques de production

Les approches automatiques de production de cartes d'occupation du sol à partir d'images de télédétection sont souvent basées sur des méthodes de classification d'images.

 

Cette classification peut être :

  • supervisée : on utilise des zones pour lesquelles on connaît l'occupation du sol comme des exemples pour un apprentissage;
  • non supervisée : on regroupe les pixels de l'image par similarité et on reconnait les classes ensuite.

La classification supervisée fournit souvent de meilleurs résultats, mais elle nécessite des données de référence pour l'apprentissage qui sont coûteuses à obtenir (campagnes sur le terrain, photo-interprétation, etc.). C'est cependant cette approche qui est utilisée dans les travaux actuels du CESBIO, comme par exemple l'édition d'une carte d'occupation des sols annuelle sur la France (avec LANDSAT 8, en attendant Sentinel-2).

 

 

L'apport du multi-temporel

Jusqu'à récemment, les cartes d'occupation du sol à échelle cartographique fine ont été presque exclusivement produites à partir d'un petit nombre de dates et ceci principalement à cause du manque de séries multi-temporelles denses fournies par des capteurs à haute résolution spatiale. L'accent était donc mis sur la richesse spectrale des images pour distinguer les différentes classes d'occupation du sol.

 

Cependant, cette approche "monodate" ne permet pas de distinguer des classes qui auraient la même signature spectrale à une date d'acquisition donnée, mais une signature différente à une autre date (des sols nus qui deviendront des cultures différentes plus tard). Pour pallier à cette difficulté, plusieurs dates peuvent être utilisées, mais cela demande une sélection spécifique de dates en fonction de la nomenclature visée.

 

Par exemple, dans l'image de gauche, acquise au mois de mai, il est très difficile de dire où sont les parcelles de colza et quelles sont les parcelles de blé. Sur l'image de droite, acquise au mois d'avril, les parcelles de colza en fleur sont très faciles à distinguer des parcelles de blé bien vert.

 

]

Image du mois d'avril. Les parcelles de colza en pleine floraison sont parfaitement visibles, elles apparaissent en jaune Image du mois de mai. Les parcelles vert clair sont des cultures d'hiver, blé ou colza principalement. Où sont les champs de Colza ?

 

Si l'on souhaite mettre en place des systèmes opérationnels et génériques (indépendants des sites cartographiés et donc des nomenclatures visées), il faut assurer une acquisition d'images fréquente et régulière. Ceci sera rendu possible par la mission Sentinel-2, et déjà, sur les données de démonstration issues de Formosat-2 et SPOT4 (Take 5). En plus, on peut montrer que le fait de disposer d'une haute résolution temporelle peut être plus intéressant que de disposer d'une grande diversité spectrale. Par exemple, la figure suivante montre des résultats de performances de classification (indice  \kappa ; plus il est élevé, mieux c'est) en fonction du nombre de dates utilisées pour la classification. On a utilisé des images Formosat-2 (4 bandes spectrales) et des simulations Vénµs (12 bandes) et Sentinelle-2 (13 bandes). On constate qu'à partir d'un nombre suffisant de dates utilisées, la richesse spectrale de Vénµs et Sentinelle-2 est rattrapée par une description fine du comportement temporel obtenu avec le simple capteur Formosat-2.

kappaVFS.png

 

 

Ce qui peut être attendu de Sentinelle-2

Sentinelle-2 a des caractéristiques uniques dans le paysage des systèmes d'observation de la Terre :

  • fauchée de 290 km.;
  • résolution spatiale de 10 à 60 m. en fonction des bandes spectrales;
  • revisite de 5 jours (avec 2 satellites);
  • 13 bandes spectrales.

Les systèmes de résolution spatiale comparable (SPOT ou Landsat) ont des revisites plus faibles et moins de bandes spectrales. Les systèmes de revisite similaire, ont une résolution spatiale plus faible (MODIS) ou des fauchées réduites (Formosat-2).

 

Avec le type de données fournies par Sentinelle-2 il est possible d'envisager le développement de systèmes de production de cartes d'occupation du sol capables d'actualiser les informations une fois par mois à l'échelle globale. La dimension temporelle, permettra de distinguer des classes dont les signatures spectrales sont très proches pendant une grande partie de l'année. La résolution spatiale améliorée permettra de travailler avec des unités minimales de cartographie plus fines.

 

Cependant, la mise en oeuvre opérationnelle de tels systèmes nécessitera une attention particulière aux besoins de validation des produits générés et aux énormes volumes de données à traiter.

 

Les cartes d'occupation produites par un tel système devront suivre une validation à échelle régionale, voire globale. De plus, comme les données de référence seront limitées, il faudra se passer au maximum de techniques d'apprentissage et essayer d'intégrer des connaissances a priori (physiques ou expertes) dans les chaînes de traitement.

 

Enfin, même si la capacité d'acquisition des nouveaux systèmes spatiaux sera améliorée, il y aura toujours des trous dans les données (nuages, par exemple). Les chaînes de traitement devront donc savoir combler ces trous, ou en tout cas y être robustes.

 

 

Les travaux du CESBIO

Danielle Ducrot, Antoine Masse et de nombreux stagiaires du CESBIO ont fabriqué récemment une grande carte d'occupation des sols sur la chaîne des Pyrénées à partir de données multi-temporelles de LANDSAT à 30 mètres de résolution. Cette carte, qui représente un vrai travail d'orfèvre, contient 70 classes. Elle a été réalisée en trois parties à partir des images peu nuageuses collectées par les satellites Landsat au cours de l'année 2010.

 

 

Carte d'occupation des sols à 70 classes obtenue à partir de séries temporelles d'images LANDSAT.

Dans sa thèse, Antoine travaille sur les méthodes qui permettent de sélectionner les meilleures dates pour réaliser une classification. De son côté, Isabel Rodes s'intéresse aux méthodes qui permettent d'utiliser toutes les images disponibles sur des zones très étendues tout en gérant les données manquantes (nuages, ombres) et le fait que tous les pixels ne sont pas vus aux mêmes dates. Ces 2 approches sont complémentaires : l'une permet de travailler avec des nomenclatures très détaillées, mais demande l'intervention d'opérateurs humains, l'autre est complètement automatique, mais moins ambitieuse en termes de détails de la classification.

 

Une troisième approche est explorée au CESBIO dans le cadre de la thèse de Julien Osman : l'utilisation de connaissances a priori de type quantitatif (à partir de données historiques) et qualitatif (connaissances d'experts thématiques) pour guider les systèmes de classification automatique.

 

Nous vous décrirons plus en détails ces différentes approches dans des billets à venir.

 

 

The atmospheric effects : how they work.

=>

Earth surface observations by space-borne optical instruments are disrupted by the atmosphere. Two atmospheric effects combine to alter the images :

  • the light absorption by air molecules
  • the light scattering by molecules and aerosols

Here are two SPOT4 (Take5) images, acquired with a time gap of 5 days above Morocco. Because of atmospheric effects, the second image has less contrast and is"hazier" than the first one.

 

 

Light Absorption :
Atmospheric absorption : in blue, the surface reflectance of a vegetation pixel, as a function of wavelength. In red, the reflectance of the same pixel at the top of atmosphere.

The air molecules absorb the light within thin absorption bands. Within these absorption bands, the reflectance measured by the satellite is lessened, and in some cases, the light may be completely absorbed and the apparent reflectance at the top of atmosphere (TOA) is zero.  (for instance, at 1.4µm, in the figure on the right. Such a property is used to detect high clouds with Sentinel-2 or Landsat-8).

Thankfully, the satellite designers usually choose to locate the spectral bands away from strong absorption bands (but beware of satellite designers ;-) ). Within the satellite channels, the absorption is generally sufficiently low so that an approximate knowledge of the absorbent abundance is enough to obtain an accurate correction of absorption. Information on absorbing gases (ozone, water vapour) concentration may be found in weather analyses.

 

Light scattering

The air molecules scatter the light. A photon that passes close to a molecule will be deflected in another direction. As the air molecules are very small compared to visible light wavelengths, they will mainly scatter short wavelengths (in the blue range). The blue sky results from the scattering of sun light by air molecules, since the blue light in the sun spectrum is much scattered while the other wavelengths are mainly transmitted to the ground. A cloud also scatters the light, but its large particles (droplets, crystals) scatter all wavelengths, which explain its white colour.

 

Apart from clouds and air molecules, scattering may be due to aerosols. Aerosols are particles of diverse nature (sulphates, soot, dust...), suspended in the atmosphere. Their abundance, type and size are extremely variable, and their effect on light is also variable. Small aerosols will mostly scatter blue light, while larger aerosols will scatter all wavelengths. Some aerosols may also absorb light. All this variability makes the correction of their effect very tricky.

The above video, provided by NASA, gives an idea of the way aerosol properties may change from one day to the other, within a two years period. The colour gives an idea of aerosol types, while the colour intensity provides the aerosol optical thickness.

Simplified model :

In a very simplified way, atmospheric effects may be modelled as follows :

ρTOA= Tgatm +Td ρsurf)

where :

  • ρTOA is the Top of Atmosphere reflectance
  • ρsurf is the earth surface reflectance
  • ρatm is the atmospheric reflectance
  • Tg is the air molecules (gazeous) transmission (Tg<1)
  • Td is the transmission due to scattering (Td<1)

When aerosol quantity increases, ρatm increases while Td decreases. These two variables also depend on view and sun angles. The closer to vertical, the lower value of ρatm and the higher value of Td .

 

Adjacency effects :

The above model should only be applied to a uniform landscape. But above a standard landscape, a heavy loaded atmosphere will also blur the images. This is explained in another post.

Models, corrections.

Several models may be used to perform atmospheric corrections. For, approximate corrections, the SMAC model should be one of the simplest. SMAC be downloaded from the CESBIO site. The difficulty in using any atmospheric correction model lies in providing the necessary information on aerosol properties. We will talk about that in another post.

Other more accurate models may be used. In our case, in the MACCS processor, we pre-compute "Look-up Tables " using an accurate radiative transfer code (Successive Orders of Scattering), that simulates the light propagation through the atmosphere. But the use of a complex model is only justified if it is possible to obtain an accurate knowledge of the aerosol optical properties.