Sentinel-2 MACCS L2A reflectance stability

As MACCS was used to produce a large amount of images (all the images acquired over France), we were able to make some more validation tests, such as checking the surface reflectance stability. To do that, the region of Arles in Provence (France), tile 31TFJ, is a good site, with a lot of clear skies, and the chance to be observed from two different orbits by Sentinel-2.

 

Here is an example for a coniferous forest (Pines probably): The reflectances of the 4 Sentinel-2 10m spectral bands are stable indeed, as it may be seen on the top plot, while the bottom plot shows the variation of the aerosol optical thickness measured by MACCS. You may have noticed that towards the summer, we can see a few couples of points, acquired 3 days apart, with small differences in surface reflectance, noticeable in the NIR infrared (black dots). These differences are probably due to directional effects as the observation from two different orbits are made with different viewing angles. These effects may be corrected. We do that for instance in our level 3 products.

Top : surface reflectances corrected with MACCS as a function of time (blue : B2 (blue), green B3 (green), red B4 (red), black B8 (Near infrared)). Bottom, Aerosol Optical Thickness measured by MACCS.

The overall stability of surface reflectances is clearly a good sign of the quality of MACCS processing.

 

MACCS/MAJA, how it works

=> 


MACCS (Multi-sensor Atmospheric Correction and Cloud Screening) is a level 2A processor, which detects the clouds and their shadows, and estimates aerosol optical thickness (AOT), water vapour and corrects for the atmospheric effects. The processor was jointly developed by CESBIO and CNES. CESBIO developed the methods and a prototype, while CNES funded the operational version of the processor, with a strong support from CESBIO for the validation.

 

More recently CNES+CESBIO and DLR joined their efforts to develop a joint processor named MAJA, for MACCS-ATCOR Joint Algorithm. MAJA is an evolution  of MACCS, in which a couple of methods inspired by ATCOR software have been added. MAJA V1_0 could have been called MACCS V6.0, but we wanted to celebrate the association of both entities with a new name.

MAJA runs at CNES within the Muscate ground segment of Theia (data are available here) and within the Venµs ground segment. And finally, MAJA is freely available in binary version for non commercial use.

 

MAJA's distinctive feature is its dedication to high resolution time series and its wide use of multi-temporal methods. For this reason, MAJA can only be applied to the optical missions which observe the earth under constant viewing angles. However, MAJA was already applied to several satellites :

 

MAJA is briefly described in the joined figure.

Atmajaspheric absorption

In the case of Sentinel-2 and Venµs, which include a water vapour channel at 940 nm (resp. 910nm) in a strong water vapour absorption band, a first step consists in estimating the atmospheric water vapour content. For the other satellites, weather analysis data can be used. After that, the processor can correct for the gaseous absorption using the SMAC model.

 

Composite imajage

The next steps deeply involve multi-temporal methods. Of course, to do that, a time series must be processed in chronological order. After each processing, a composite image is updated with with the unclouded pixels from the processed date. This composite image is used as a reference for the cloud detection and the AOT estimate.

 

The cloud majasks

Our cloud detection method is based on a large number of tests, the most efficient of which are :

  • a test based on the  cirrus band (at 1380 nm), available on Landsat 8 and Sentinel-2, which detects very well the high clouds (above 2000m)
  • a multi-temporal test, which detects a steep increase of the blue surface reflectance, which is the sign of presence of a cloud.
  • and finally, to avoid over detections of clouds, for each potential detected by one of the previous tests, a last test measures the correlation of the pixel neighbourhood with the previous images. As it is unlikely that two different clouds at the same location on successive dates have the same shape, if a large correlation is observed, the pixel is finally not declared as a cloud.

Having detected the clouds, we can follow with the detection of  cloud shadows, water, and snow.

Aerosol optical thickness estimajate

The aerosol optical thickness (AOT) estimate combines several criteria in the computation of a global cost function, which is then inverted using non linear least mean squares inversion.

  • A multi-temporal criterion : after atmospheric correction, two successive observations of the same neighbourhood should provide nearly the same surface reflectances. The squared residuals after atmospheric correction are inserted in the cost function.
  • A multi-spectral criterion : above vegetation, and also above many bare soils, the surface reflectance in the blue is close to half the reflectance in the red. The squared residuals to this relation after atmospheric correction are also added to the cost function.
  • Optical Thickness minimum and maximum : AOT cannot be negative, and should not get higher that the one measured using the dark pixel method. When the AOT values are above maximum or under minimum, a high cost is added to the cost function.

 

The cost function evaluation is evaluated using neighbourhoods of coarse resolution pixels (240m), spreading over 2 kilometres. The obtained AOT images are then smoothed, the gaps are filled to obtain finally an AOT map with a 5 km resolution. The aerosol type is not estimated, it is a processing parameter which can be fixed per geographic region.

Atmajaspheric correction

One of he quicklooks we produce with each image for visual verification, here for Chiapas site in Mexico, with the TOA reflectance, top left, the AOT and cloud mask, bottom left, surface reflectance after adjacency effect correction, top right, and the same with slope correction, bottom right.

 

Once the AOT is known, we can retrieve the surface reflectances. To do that, we are using look-up tables (LUT) which are computed using the SOS radiative transfer code (Successive Orders of Scattering, Lenoble, 2007). These LUT are also used in the AOT estimation.

 

The surface reflectance of the cloud free pixels obtained there are used to update the composite image, which will be used for the processing of the next image in the time series.

 

Before editing the output product, we still need to correct for two other points, already described in this blog : the adjacency effects and the effects of terrain slopes on the illumination.

 

MAJA development started in 2005, and the contributor list is starting to be quite long :

  • at CESBIO : H.Tromp, V. Debaecker, M. Huc, P.Gely, Bastien Rouquié and O.Hagolle,
  • at CNES : B. Petrucci, D.Villa-Pascual, Camille Desjardins, Pierre Lassalle
  • at DLR : A. Makarau and R.Richter
  • at CS-SI : T.Feuvrier, C.Ruffel, A.Bricier and many others
  • at CAP GEMINI : M.Farges, G. Rochais, E.Durand
  • at Magellium : E. Hillairet

For more details, we have published 4 papers about MACCS methods and validation :

  • A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENµS, LANDSAT and SENTINEL-2 images, O Hagolle, M Huc, D. Villa Pascual, G Dedieu, Remote Sensing of Environment 114 (8), 1747-1755
  • Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images, O Hagolle, G Dedieu, B Mougenot, V Debaecker, B Duchemin, A Meygret, Remote Sensing of Environment 112 (4), 1689-1701

 

SAT-IRR: Satellite for Irrigation Scheduling

=>

Suite à l'expérience de pilotage d'irrigation menée au Maroc lors de l'expérience SPOT4-Take5 (Le Page et al, 2014), un outil Web d'aide à la prise de décision d'irrigation est en cours de développement (http://osr-cesbio.ups-tlse.fr/Satirr). L’outil est fonctionnel sur trois tuiles Landsat8 situées à Marrakech au Maroc, Kairouan en Tunisie, Toulouse en France.

 

L'outil s'adresse à des irrigants :  après avoir dessiné sa parcelle sur un fond cartographique, l’utilisateur répond à 4 questions. Il choisit sa culture parmi 7 options actuellement renseignées (maïs, blé, olivier…), son sol parmi les 12 sols type de l’USDA, sa date de semis et son mode d’irrigation (gravitaire, aspersion ou goutte à goutte). Cette initialisation sommaire est suffisante pour lancer le service, mais l’utilisateur pourra modifier à tout moment les contours de sa parcelle ou affiner la paramétrisation s’il connaît bien le sol, les particularités de sa culture, etc.

 

Dans un premier temps le serveur se charge de faire une approximation d’un comportement

satirr_screenshot

Screenshot from SAT-IRR web interface. The four icons allow modifying the plot parameters and contours, input irrigation, and consulting the results as graphs or tables. The graph results show a small Openlayers window with the last NDVI image, the sequence of NDVI thumbnail images, and 4 graphics: The “atmospheric part (rainfall, Reference Evapotranspiration and actual evapotranspiration), the second graph shows the status of the soil water content separated in three layers, the third graph shows the evolution of Basal Crop Coefficient and Fraction Cover, and the last graph is NDVI. The blue square at the right of the graphs are projections for the next month, including the green bars which are irrigation recommendations

moyen de la plante. Pour cela, une climatologie mensuelle est compilée (moyenne multi-annuelle de paramètres météo) puis interpolée au pas de temps journalier, alors que le comportement moyen de la plante est tiré des tables du document FAO-56 « FAO Irrigation and Drainage n°56: Guidelines for Computing Crop Water Requirements » (Allen et al, 1998). Dans un second temps, les images satellites déjà présentes sur le serveur sont examinées puis les relations entre NDVI et Coefficient Cultural de Base (Basal Crop Coefficient, Kcb) et le pourcentage de la couverture du sol par la végétation (Fraction cover, Fc) sont déterminées à chaque date disponible en faisant une moyenne sur la parcelle.

 

La météo passée est renseignée par les mesures effectuées sur la station synoptique de l’Organisation Mondiale Météorologique la plus proche, et synthétisée quotidiennement sous la forme de l’évapotranspiration de référence (ET0) et de la pluie. Enfin, des prévisions météo sont obtenues grâce à l’API de l’Institut Météorologique Norvégien.

 

Finalement, un bilan hydrique très proche de celui décrit dans la méthode FAO-56 est calculé en combinant ainsi comportement cultural et climatologie type, imagerie satellitaire, mesures et prévisions météo ainsi que projection dans le futur du développement de la culture. Le but étant bien évidemment de proposer une date et dose d’irrigation.

 

En plus de mettre à jour la météo (mesures et prévisions), le serveur vérifiera chaque jour la disponibilité de nouvelles images (uniquement Landsat8 pour le moment). Si une nouvelle image est disponible, elle est téléchargée, corrigée des effets atmosphériques en utilisant les informations fournies par le photomètre du réseau Aeronet le plus proche en utilisant le code SMAC (Rahman & Dedieu, 1994), puis un masque de nuage est créé et le NDVI est calculé. Cette image est stockée alors que le fichier original est jeté pour ne pas encombrer le serveur.

 

L’ensemble paramétrisation/mesures/prévision est stocké sur une base postgres/postgis qui fait le lien avec une interface web. L’utilisateur peut consulter les résultats sous forme de tableaux ou de graphes, et rajouter ses propres irrigations dans une autre interface dédiée.

 

Bien que l’interface soit encore un peu fruste, nous envisageons surtout des développements du côté serveur:

  • Adaptation à Sentinel-2 : à priori le passage à S2 ne devrait pas poser de soucis. Il faudra cependant adapter le calcul des tuiles à télécharger, le code de téléchargement, ainsi que la lecture du format.
  • Utilisation de Sentinel-1: Dans l’état actuel, le bon fonctionnement du bilan hydrique repose sur l’information réelle de l’irrigation que doit fournir l’utilisateur. Nous prévoyons de tester l'utilisation d’images S1 pour déterminer les dates d'irrigation.
  • Accès à des stations agro-météo locales : Dans le cadre du développement du Système d’Information Environnemental au Cesbio, la télémétrie de plusieurs stations météo se met petit à petit en place (par exemple voir http://trema.ucam.ac.ma (Jarlan et al, 2015)), nous comptons rendre ces stations accessibles à travers un service web normalisé du type Sensor Web.
  • Introduction de réseau d’irrigation collective. Les travaux de thèse de Kharrou (2013) et Belaqziz (2013, 2014) ont montré que la télédétection spatiale peut servir à optimiser les tours d’eau sur un secteur irrigué. Nous comptons donc offrir la possibilité d’introduire un ensemble de parcelle pour l’associer à un réseau de distribution et proposer in fine un arrangement optimisé du tour d’eau. Cependant, à l’heure actuelle, cet objectif est plutôt de l’ordre du défi !
  • Nous travaillons actuellement sur une procédure d'estimation du rendement du blé avec la télédétection spatiale (Thèse J. Toumi) et espérons ainsi introduire une estimation précoce du rendement dans cet outil.

Si vous souhaitez essayer l'outil, inscrivez-vous, c'est gratuit. Si les régions de test ne vous conviennent pas, contactez-moi!

Références:

  1. Le Page M., J. Toumi, S. Khabba, O. Hagolle, A. Tavernier, M. Kharrou, S. Er-Raki, M. Huc, M. Kasbani, A. Moutamanni, M. Yousfi, and L. Jarlan, “A Life-Size and Near Real-Time Test of Irrigation Scheduling with a Sentinel-2 Like Time Series (SPOT4-Take5) in Morocco,” Remote Sens., vol. 6, no. 11, pp. 11182–11203, Nov. 2014.
  2. Allen R., L. Pereira, D. Raes, and M. Smith, FAO Irrigation and Drainage n°56: Guidelines for Computing Crop Water Requirements, no. 56. FAO, 1998, pp. 273–282.
  3. Rahman H. and G. Dedieu, “SMAC: a simplified method for the atmospheric correction of satellite measurements in the solar spectrum,” Int. J. Remote Sens., vol. 15, no. 1, pp. 123–143, 1994.
  4. Kharrou M.H., M. Le Page, A. Chehbouni, V. Simonneaux, S. Er-Raki, L. Jarlan, L. Ouzine, S. Khabba, and A. Chehbouni, “Assessment of Equity and Adequacy of Water Delivery in Irrigation Systems Using Remote Sensing-Based Indicators in Semi-Arid Region, Morocco,” Water Resour. Manag., vol. 27, no. 13, pp. 4697–4714, Sep. 2013.
  5. Belaqziz S., S. Mangiarotti, M. Le Page, S. Khabba, S. Er-Raki, T. Agouti, L. Drapeau, M. H. Kharrou, M. El Adnani, and L. Jarlan, “Irrigation scheduling of a classical gravity network based on the Covariance Matrix Adaptation – Evolutionary Strategy algorithm,” Comput. Electron. Agric., vol. 102, pp. 64–72, Mar. 2014.
  6. Belaqziz S., S. Khabba, S. Er-Raki, L. Jarlan, M. Le Page, M. H. Kharrou, M. El Adnani, and A. Chehbouni, “A new irrigation priority index based on remote sensing data for assessing the networks irrigation scheduling,” Agric. Water Manag., vol. 119, pp. 1–9, Mar. 2013.
  7. Jarlan L., S. Khabba, S. Er-raki, M. Le Page et al, “Remote sensing of water resources in semi-arid Mediterranean basins: The Joint International Laboratory TREMA,” Int. J. Remote Sens., vol. (under review), 2015.

How MACCS estimates Aerosol Optical Depth.

=>

I already explained in this blog  the principles of estimation of the aerosol optical thickness that we use to process the LANDSAT or SPOT5 (Take5) data, and soon Venµs or Sentinel-2 data, within the MACCS method developed at CESBIO, and used at CNES by THEIA. We started writing an article in 2010 to explain the method details and show the validation results, but we only found a sufficiently quiet period this autumn to finish it. The paper has just been published in remote sensing (MDPI), with open access. Enjoy your reading !

 

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.

From left to right, validation results for Aerosol Optical Thickness (AOT) measures by the multi-temporal method, the multi-spectral method and the combination of both. The combination of both methods allows to measure AOT in a much larger range of cases without degrading accuracy.

USGS releases LANDSAT 8 data corrected from atmospheric effects

=>

Since Christmas (almost, December 23rd), USGS distributes LANDSAT 8 data corrected from atmospheric effects, as they already did for LANDSAT 5 and 7. These data have been long awaited and will be very useful.

 

Of course, USGS processing addresses the whole world. To access these products, you should use the earthexplorer.server. You will find the data by clicking on LANDSAT_CDR, as shown in the image on the right. You will have to order for their processing before being able to download them. Usually this processing is very quick.

 

We did not yet compare the USGS products developed by NASA (E.Vermote) with our Level 2A LANDAT 8 data acquired over France and produced and distributed by THEIA with CESBIO processor, but we will work into it shortly.

 

LANDSAT 5 & 7 acquired above France since 2009 soon released by THEIA

=>

As for LANDSAT 8 a few weeks ago, we just produced the level 2A products for the LANDSAT 5 et LANDSAT 7 data acquired above France from 2009 to 2011. This data set will be released in a few days, when its transfer to the distribution server has ended. The MUSCATE team took charge of the processing for the THEIA Land data center, using CNES computing center. The data will be available on the following site :

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

Example of a Level 2A product from LANDSAT 5 over the Atlantic coast of France. The clouds are circled in green, the water mask in blue, and the snow in pink. Sometimes, the water turbidity provides a signal similar to snow in the infra-red, which turns the now flag on...

The processing methods and the data format are similar to the LANDSAT 8 data set described here. However there are also a few differences which are detailed below :

 

Starting point.

The starting point is not the same for LANDSAT 5 and LANDSAT 7 :

  • For LANDSAT 7, as for LANDSAT 8, we start from Level 1T products provided by USGS. These products have a huge defect, with black stripes appearing away from the center axis of the image. These stripes are due to the breakdown of a mirror in 2003. The origin of this defect is described here. In our case, we decided to use only the central part of the images, doing a slight interpolation when the stripes are thin, and removing the data when the lack stripes are too large.
  • LANDSAT 5 data acquired above France are not yet available at USGS. ESA owns these data and agreed to provide them  (Thanks to Bianca Hoersh and Riccardo Biasutti from ESA, and to the SERCO company who processed the data). As a result, this data set is a unique data set, only available online here and nowhere else ! These data are provided at level 1G, for which the data have not yet been ortho-rectified.  We had to ortho-rectify them at Theia, using CNES's SIGMA tool, as for SPOT4 (Take5).
  • Having a different approach for both sensors has a drawback. The grounc control point data base used at USGS seems to have some errors in France, and for instance, the location errors near Toulouse have a bias of about one pixel southward. It is not the case for the LANDSAT5 data ortho-rectified by THEIA, and therefore, one may observe registration errors in a time series involving LANDSAT 5 and LANDSAT 7 images. ESA's data also have some defects, which are presented at the end of this post.
Resampling to Lambert'93 projection

Level 1T data are provided with the UTM projection. This projections uses three different zones over France, for which the registration of data is not direct. We decided to resample the data on a Lambert'93 projection, which is the official French projection. Of course, the LANDSAT5 have been directly projected in Lambert 93.

Tiling of products

We chose to tile the data in 110*110 km tile s spaced with a 100 km interval, as it will be done for Sentinel-2. The (1,1) tile is in the SouthWest corner of France. The tile of Toulouse is the 5th to the West, and the 2nd to the North. It is named D0005H0002 (D for "droite", H for "Haut"). For Corsica, a different tiling made of 2 tiles was defined.

 

For each tile, we provide the whole set of dates for which a LANDSAT 5 or 7 image intersects the tile. A few date may be missing, for several reasons, in general related to the cloud cover :

  • The image was not acquired by LANDSAT 5 or 7 (when a 100% cloud cover is forecast, the image is not acquired).
  • The image was acquired but not processed to L1T by LANDSAT7 at USGS, or to L1C at CNES, because the cloud cover prevented from using a sufficient number of ground control points
  • The Level 2A processing rejects images with more than 90% of cloud cover.

 

Level 2A processing (atmospheric correction and cloud screening)

First of all, we would like to outline that our processor does not process the themal bands of LANDSAT

 

For the visible, near and short wave infrared bands, we use the same method as for SPOT4(Take5). It involves also the MACCS processor, developed and maintained by Mireille Huc at CESBIO. It is based on multi-temporal methods for cloud screening, cloud shadow detection, water detection as well as for the estimation of the aerosol optical thickness.

 

However, thanks to LANDSAT spectral bands, our processing was enriched compared to SPOT4 (Take5). Thanks to the blue band, we have an additional criterion to detect the aerosols, thanks to the quasi constant relationship between the surface reflectances in the blue and in the red above vegetation. The precision gain due to this criterion compensates for the precision loss due the lower repetitivity of  LANDSAT images. Finally, as there is no 1.38 µm band on LANDSAT 5 and 7, the detection of high clouds is much less easy than for LANDSAT 8.

 

Images of one of the Atlantic Coast tiles, coming from different LANDSAT Paths (left and middle, tracks 201 and 200). The viewing angles are slightly different as the left image was observed from the West and the rmiddle image from the East. On the right, a Landsat7 image from track 201 reduced to its central part.

To increase repetitivity of observation which is essential in our multi-temporal method, we decided to use time series that merge LANDSAT 5 and LANDSAT 7 data as well as LANDSAT 5 data coming from adjacent tracks. As these data are not observed under the same viewing angle (+/- 7 degrees), but the angle difference is small enough to increase precision on the overlap zones, even if it may cause the appearance of artefacts in the AOT images.

 

Data Format

For LANDSAT 5 and 7, we used the same format as for SPOT4 (Take5).

 

Known issues :

Here is a little list of known defects for THEIA's LANDSAT 5 and 7 L2A products :

Example of LANDSATV5 "afterglow" issue near a large cloud. This electronic issue takes the appearance of whiter stripes above vegetation.

- reference data for ortho-rectification at USGS may be biased by more than 30 m (38 m in Toulouse). The Landsat 8 data could be misregistered with the LANDSAT 5 data ortho-rectified at CNES using a national geographical reference.

- LANDSAT 5 TM instrument electronics have an "afterglow" issue, that causes the appearance of whiter stripes perpendicularly to the satellite track near very bright zones such as a large cloud.

- ESA's LANDSAT 5 products have some random bright spots that appear as colored spots in color composites.

Bright color spots observed on some ESA LANDSAT 5 images.

- in LANDSAT products, the "nodata" value that tells if a pixel is outside the image is 0, which is also a value observed within the image. Sometimes pixels may be identified as nodata when the are in fact within the image. It happens mainly over sea, where the medium infrared reflectance is often equal to zero. In this case, all the bands have the nodata value which, in our products is -10 000, to avoid the same difficulties for subsequent users.

Completion of the processing of LANDSAT 8 Level 2A products taken above France in 2013.

=>

That was fast ! The processing of all the LANDSAT 8 images taken above France in 2013 took less than 15 days. The first LANDSAT 8 images were taken in April 2013. The MUSCATE team processed the data for the THEIA land data center, using CNES computing center.

 

A few more days will be necessary to upload the data on the THEIA website and to check that the data are correct. Finally, the longest part in the processing is the downloading of the input Level 1T products from USGS earthexplorer website (equivalent to the Level 1C in THEIA's nomenclature).

 

The Level 2A data quality is quite good, as may be seen on the browse products on the right, as shown by images on the right, which come from the times series obtained on the tile of Paris. As usual on this blog, the clouds are circled in green, the shadows in black, the snow in pink and the water in blue. A few clouds are sometimes missed by our multi-temporal method, when the repetitivity of cloud free acquisition is too low, as in the image on the right which was acquired during a cloudy spring. The following images in the time series are not affected by this kind of defect.

 

This paper aims at describing the main steps of the processing.

 

Starting point : Level 1T

We download the input data from the earthexplorer website, using an enhanced version of the script described here. These products are ortho-rectified by USGS, using a global data base of ground control points.

 

The location requirement for LANDSAT 8 is 50 m, which seems to be met by the L1T products. We found location errors around 1.5 pixels near Toulouse, but most regions seem to have better performances. USGS confirmed a 38m bias Southward near Toulouse and will try to correct them.  Our processing does not correct for these small errors, and the next version of the USGS LANDSAT 8 processing only wil lcorrect for this bias.

 

Regarding LANDSAT 8 absolute radiometric calibration, we use the coefficient values recommended by LANDSAT 8 and provided with the L1T products.

 

Resampling to Lambert'93 projection


Level 1T data are provided with the UTM projection. This projections uses three different zones over France, for which the registration of data is not direct. We decided to resample the data on a Lambert'93 projection, which is the official French projection.

Tiling of products

We chose to tile the data in 110*110 km tiles spaced with a 100 km interval, as it will be done for Sentinel-2. The (1,1) tile is in the SouthWest corner of France. The tile of Toulouse is the 5th to the West, and the 2nd to the North. It is named D0005H0002 (D for "droite", H for "Haut")

 

For Corsica, a different tiling made of 2 tiles was defined.

 

For each tile, we provide the whole set of dates for which a LANDSAT 8 image intersects the tile. A few date may be missing, for several reasons, in general related to the cloud cover :

  • The image was not acquired by LANDSAT 8 (when a 100% cloud cover is forecast, the image is not acquired).
  • The image was acquired but not processed to level 1T by LANDSAT8, because the cloud cover prevented from using a sufficient number of ground control points
  • The Level 2A processing rejects images with more than 90% of cloud cover.

 

Level 2A processing (atmospheric correction and cloud screening)

First of all, we would like to outline that our processor does not process the themal bands of LANDSAT 8.

For the visible, near and short wave infrared bands, we use the same method as for SPOT4(Take5). It involves also the MACCS processor, developed and maintained by Mireille Huc at CESBIO. It is based on multi-temporal methods for cloud screening, cloud shadow detection, water detection as well as for the estimation of the aerosol optical thickness.

 

However, thanks to LANDSAT 8 spectral bands, our processing was enriched compared to SPOT4 (Take5) : LANDSAT8's 1.38µm band enables an enhanced detection of high and thin clouds. And thanks to the blue band, we have an additional criterion to detect the aerosols, thanks to the quasi constant relationship between the surface reflectances in the blue and in the red above vegetation. The precision gain due to this criterion compensates for the precision loss due the lower repetitivity of  LANDSAT8 images.

Level 2A images from Paris's tile, from 3 different LANDSAT 8 tracks (From left to right, tracks 200, 199, 198). The viewing angle differs as the image is from the west on the left image, at nadir in the center and from the east for the right image.

 

To enhance the cloud screening accuracy, we decided to use the data from adjacent satellite tracks within the same time series. These data are not acquired under exactly the same angle (+/- 7 degrees), which is the assumed by the multi-temporal method, but the difference is small enough to allow a large accuracy gain due to the enhanced repetitivity. However, because of this approximation, a few artefacts may be observed.

 

For a greater enhancement, we might also use LANDSAT 7 and LANDSAT 8 data in the same time series, but we will implement that later on...

 

Data Format

For LANDSAT 8, we used the same format as for SPOT4 (Take5), excepted a few details, that I will describe soon...

 

 

 

New satellites added to SMAC atmospheric correction

=>

New coefficients have been added to the CESBIO repository for SMAC coefficients. The new sensors taken into account are :

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

 

The Simplified Model for Atmospheric Correction (SMAC) is the perfect model to perform easy, quick and not too dirty atmospheric corrections. It is based on very simple analytic formulas, based on the 5S model. The 49 coefficients of this model are fitted using a large number of radiative transfer simulations with the 6S model (the old historic version, not the recent vector version). This software is not very accurate (much less than MACCS), and it requires in-situ measurements for the aerosol optical thickness, and weather analyses for ozone and water vapour. If these data are available,  in most cases, its accuracy is within 2 and 3 percent, if we do not account for adjacency effects and slope effects, and it may be worse for large viewing and solar angles (above 70°) or within strong absorption bands.

 

SMAC is very easy to use:

#read the 49 coefficients in smac_soefs table
nom_smac ='COEFS/coef_FORMOSAT2_B1_CONT.dat'
coefs=coeff(nom_smac)
 
#read the TOA reflectance image in r_toa variable
#depends on the file format
 
#read the angle values in the image metadata
theta_s=30
phi_s=180
theta_v=0
phi_v=0
# compute pressure at pixel altitude
pressure=PdeZ(1300)
 
#find the values of AOT, UO3, UH2O
AOT550=0.1
UO3=0.3
UH2O=3
 
#compute the atmospheric correction
r_surf=smac_inv(r_toa,theta_s,phi_s,theta_v,phi_v,pressure,AOT,UO3,UH2O,coefs)

where :

  • theta_s, phi_s are the solar zenith and azimuth angles
  • theta_v, phi_v are the viewing zenith and azimuth angles
  • AOT is the aerosol optical thickness at 550 nm which may be obtained from an Aeronet stations, or guessed, or equal to 0.1 for a really dirty atmospheric correction.
  • UO3 is the ozone content in cm.atm (0.3 is OK)
  • UH2O is the water vapour integrated content in kg/m². When I do quick and dirty atmospheric correction, I often use a value equal to 3, but I do not process spectral bands with strong water vapour absorption bands.
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.

 

SPOT4(Take5) aerosol optical thickness validation results

We are currently preparing a data reprocessing of all SPOT4 (Take5) data, to be released before the end of 2013. For this, I tested several aerosol models and compiled all the validation results for our multi-temporal Aerosol Optical Thickness (AOT) estimation method named MACCS. Our estimates are compared to AERONET in-situ AOT measurements.

The MACCS method applied to SPOT4(Take5) data, which lacks a blue band, uses two procedures to estimate AOT :

  • either the AOT is estimated by a multi-temporal method
  • or it is gap-filled. The presence of gaps may be due to clouds, water or snow, or because the pixel reflectance is too-high for an accurate estmate, or because of a too large variation of reflectance with time is detected.

 

Comparison of MACCS AOT estimates with the in-situ measurements from AERONET. The blue dots correspond to cases for which the atmosphere is stable and for which there are no clouds in the neighborhood of the AERONET site. The red dots correspond to situations when the AERONET optical thickness varies around the satellite overpass time, or when clouds are detected in the image neighbourhood (20*20 km).
On the left plot, only the dates and sites for which less than 60% of the pixels were gap-filled; wheras the right plot only tolerates 20% of gap-filled pixels. The gap-filling method does not seem to introduce large amount of errors in wases when the atmosphere is stable, but it is less accurate in unstable cases..

 

The aerosol estimates have been obtained with MACCS prototype which is developed and maintained by Mireille Huc at CESBIO. The aerosol model is not the same as the one used for SPOT4 (Take5) first processing. This model is based on greater particles (with a modal radius of 0.2µm, compared to 0.1µm in the initial processing), as it provides a better overall agreement with AERONET measurements. We will use this model for most sites for SPOT4(Take5) reprocessing.

 

The RMS error of AOT estimates is 0.06, which is a state of the art performance, obtained in a very difficult condition with no blue band available. Moreover, in order to show more validation points, a few validation sites (Bruxelles, Gwangju, Ouarzazate, Wallops, NASA_LaRC) are in fact distant by more that 60 kilometers from the image footprint, which tends to degrade the performances.

 

The AERONET sites used in this study are :

 

SPOT4 Take5 Site
Aeronet Site
Belgium Brussels
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

 

The worst results are obtained for the following sites :
  • Gwangju (Korea): The SPOT footprint in on the coast, while the AERONET site is 70 km inland, near a large town.
  • Ben Salem (Tunisia): this site was very cloudy in Spring, and large reflectance variations are observed between the remaining clear dates.
  • Palaiseau and Paris : In that case, the aerosol model seems to be inappropriate, and absorbing pollution aerosol should be introduced.

On the contrary, several sites provide very accurate results, for instance in Morocco (even the desertic Ouarzazate), Provence (including the Frioul Island where the AOT is extrapolated from the coast), and also Sudmipy, Wallops et Cart_site. Some SPOT4 (Take5) users reported inaccuracies on some tropical sites but we do not have an AERONET validation site near these SPOT4(Take5) sites.

 

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.