MUSCATE S2 product versions

V1_0 was a preliminary processing performed with early version of MACCS adapted to Sentinel-2. It was plagued with a bug in the cloud shadow mask : when more than 255 clouds were present in one image, shadow detection went completely wrong and shadows were detected anywhere. V1_0 was not produced in MUSCATE operationnal context, but in an earlier validation context. Production was available over France only.

 

V1_1 replaces V1_0 over France, Production started mid November, but was only released in February because of many difficulties encountered by MUSCATE. It corrects for the bug related to cloud shadows observed with V1_0, and  was fully processed by the operational MUSCATE center.

 

V1_2 was used in another context not related to MUSCATE, you will not find any product with this version number within MUSCATE server.

 

V1_3 was applied starting from February 2017, to data sets produced above Reunion, Burkina, Senegal, Tunisia, Morocco.The aerosol estimate is improved compared to V1_1. Together with Bastien Rouquié, from CESBIO, we worked on the tuning of the blue-red ratio which is used in the multi-spectral method to estimate aerosols (which is combined with a multi-temporal method). Initially, we used bands B2 (Blue) and B4 (Red), with a ratio of 0.45. We found out that better results were obtained with B1 (Blue) and B4 (Red), still with a ratio of 0.45. More accurate studies tend to recommend a higher value, closer to 0.5

 

As can be seen in the figure below, the estimates obtained with V1_3 are not biased anymore, and have a reduced standard deviation. As a consequence of getting a lower aerosol optical thickness, the surface reflectances of V1_3 products are 2% higher than those of the earlier versions (some user feedbacks from V1_0 said reflectances were sometimes too low).

 

Using B2/B4 ratio

Using B1/B4 ratio

 

V1_4 will be provided with a better shadows mask, as it has been found that the shadows masks are too severe and contain too many commission errors. It is in its final stages of validation.

 

 

Sentinel-2A unveils two new lakes in Borneo

The remote sensing company Planet has published images, which were captured by its fleet of nano-satellites, of Lake Bakun on the island of Borneo [1]. Planet tells us that this artificial lake is the largest lake in Southeast Asia since its flooding in 2011. Yet it is not visible on mapping web services like Bing Maps or Google Maps! It does show in OpenStreetMap at low resolution. This is due to the fact that the island of Borneo is one of the most cloudy places in the world so the usable optical remote sensing images are very rare.
Continue reading

Radiometric quantities : irradiance, radiance, reflectance

=>

Beware, this document contains equations

Radiometric quantities are numerous and may seem rather complex at first sight. Here is a little guide of the various quantities you might stumble upon when using optical remote sensing data. It might be a little boring, it is full-up with integral formulas, but well, to really understand what images mean, it might be useful. Come on, let's go, and cheers !

Continue reading

New shadows for SPOT5 (Take5)

=>

Don't worry, everything is fine for SPOT5 (Take5) ! While SPOT5 is peacefully getting some well deserved rest, THEIA just reprocessed the Take5 data to solve a little list of issues detected after the first processsing. The data should be released very soon. The most significant change involves the cloud shadows mask computed by MACCS processor . The masking method was completely changed, from floor to ceiling. It is the last work done by Mireille Huc on this processor before leaving for other projects.THEIA's users should be grateful to Mireille for all the excellent work she did on this method and software.

 

Shadow detection is even more difficult than cloud detection, because dark zones are quite frequent in most landscapes, and are more frequent than bright zones. As explained in this post, MACCS detects the "potential" shadows via a threshold on the darkening of pixels, then checks if the "potential shadows" may be traced back to a cloud.

 

Old method

The previous version was roughly based on the following steps :

1 - cloud detection

2 - compute darkening compared to reference image (a composite image made from the most recent cloud free observation for each pixel)

3 - for each altitude level between 500 and 10 000 m

  • geometric projection of the cloud shadow and computation of the mean darkening over the zone covered by the projection (excluding the zones covered by clouds).
  • select the H0 altitude with the maximal darkening

4. This altitude H0 is affected to all clouds and allows to compute the location of the shadows of all clouds

an additional step was added to account for cases when all the clouds within image are not at the same altitude

5. for all clouds :

  • check if the shadow zones correspond to a real darkening
  • for all couds for which darkening is too low
    • reiterate step 3 to find the cloud

 

This method worked quite well for small clouds, but could frequently fail for big clouds which masked a part of the shadow which lied beneath them. In this case, it was difficult to match the cloud and the shadow, as for instance in the image below :

On all the images showed here, the detected clouds are circled in green and the shadows are circled in black. On this image obtained with the previous method, most of the shadows are not detected.

New method

To avoid these issues, we have decided to stop trying to match directly clouds and shadows. The search of shadows in only restricted to zones where the detected clouds can cas shadows, and a final verification can solve the cases when the surface of the shadow is larger than the surface of the cloud.

 

 

The new version roughly works as follows :

1 - Clouds are detected

2 - for each cloud, search zone where clouds between 500 10000m altitude can cast shadows

  • search for darkened pixels within this zone
  • If total surface of darkened pixels is greater than cloud surface, adjust detection threshold to select darkest zones and reduce shadow surface to about the surface of the cloud

These images show the shadow detection from the previous version on the left, and from the new version on the right. Please click on the images to see them at a higher resolutions and observe the enhancement which is quite general. From top to bottom, the images come from Baotou, in China, Belgium and Central African Republic.

Image edges

On the image edges, shadows may have been projected by a sneaky cloud which is hidden outside the image. Up to now, within the zone where this can happen, we used to only do a simple threshold on darkening, with a conservative value to avoid false detections, which led to frequent omissions. Here, Mireille added a clever improvement which consists in guessing the darkening threshold from the darkening of the shadows obtained in the middle of the image, where we are more confident in the shadow detection quality.

In this case (Baotou, China), the shadow in the North West corner, which was undetected by the previous version has been detected now. But, because of a bug found too late, as the production was already running, sometimes the processor omits the last shadow in the list, after running through the images from left to right and from top to bottom. It is the case for the last clouds within this image.

Generalisation

This improvement, first implemented in MACCS prototype, which is used to produce the take5 images, is also being used to process the LANDSAT images produced to level 2A over France since the beginning of 2016, but we will not reprocess the older data. The improvement is also being implemented within MACCS operational version 5.0 and it will be applied to Sentinel-2 within MUSCATE.

Special case

We were unable to solve the issue of the image below, for which the apparent shadow is in fact a black cloud due to the explosion of a petrol storage facility !

 

De nouvelles ombres pour SPOT5 (Take5)

=>

Ne vous inquiétez pas, tout va bien pour SPOT5 (Take5) ! Alors que le satellite SPOT5 vit une retraite paisible, Theia vient de retraiter les données Take5 pour résoudre un petit nombre de problèmes détectés lors du premier traitement. Les données devraient être mises en ligne dans les jours qui viennent. La modification la plus importante concerne le masque des ombres de nuages, obtenu par la chaîne MACCS. Ce masque a été complètement revu, du sol au nuage. Il s'agit du dernier travail effectué par Mireille Huc sur cette chaîne avant son départ pour d'autres projets (les utilisateurs de MACCS pourront dire un grand merci à Mireille pour toute son oeuvre, dont cette dernière amélioration).

 

La détection des ombres est encore plus difficile que la détection des nuages, car les zones sombres sont beaucoup plus présentes que les images claires dans les images. Comme expliqué dans cet article, MACCS détecte les ombres "potentielles" par un seuillage sur l'assombrissement des images, puis vérifie qu'il arrive à faire la correspondance entre chaque ombre et le nuage dont elle provient.

 

Ancienne méthode

La version précédente fonctionnait (en gros) de la manière suivante :

1 - détection des nuages

2 - calcul de l'assombrissement de l'image par rapport à l'image de référence (image composite formée du dernier pixel clair disponible)

3 - pour chaque altitude entre 500 et 10 000 m

  • on calculait la projection de l'ombre des nuages et on calculait l'assombrissement moyen sur cette zone (en excluant les zones où l'ombre est masquée par le nuage)
  • on sélectionnait l'altitude H0 présentant le plus fort assombrissement moyen

4. on attribuait cette altitude H0 à tous les nuages et on projetait leur ombre

#étape supplémentaire pour tenir compte des cas où tous les nuages ne sont pas à la même altitude

5. pour tous les nuages :

  • on vérifiait si l'altitude correspondait bien à une ombre, sinon
  • pour tous les nuages dont l'ombre ne correspondait pas
    • On cherchait la bonne altitude par la méthode 3

 

Cette méthode marchait très bien pour de petits nuages, mais échouait fréquemment pour les grands nuages dont une partie de l'ombre est masquée par le nuage. Il était en effet difficile d'apparier les deux zones, comme par exemple dans l'image ci-dessous :

Sur toutes les images présentées ici, les nuages détectés sont entourés en vert et les ombres détectées en noir. Sur cette image obtenue avec l'ancienne méthode, la plupart des ombres du grand nuage ne sont pas détectées.

Nouvelle méthode

Pour éviter ces problèmes, nous avons décidé de ne plus essayer d'associer directement nuages et ombres. On réduit cependant la recherche des ombres aux zones où les nuages détectés peuvent produire des ombres, et une vérification finale permet de résoudre les cas où la surface des ombres détectées est plus grande que la surface du nuage.

 

 

La nouvelle version fonctionne (en gros toujours) de la manière suivante :

1 - on détecte les nuages

2 - pour chaque nuage, on recherche la zone où peuvent se cacher des ombres (pour une altitude allant de 250 à 10 000 m)

  • On cherche les zones qui se sont assombries à l'intérieur de ces zones
  • Si la surface de ces zones est supérieure à la surface des nuages, on ajuste le seuil de détection pour sélectionner les zones les plus sombres et donc réduire la surface

Ces images présentent, à gauche, l'ancienne version, et à droite, la nouvelle version.. N'hésitez pas à cliquer sur les images pour les voir à pleine résolution et pour observer l'amélioration, qui est assez générale. De haut en bas, les images choisies viennent des sites de Baotou (Chine), Belgique et République Centre Africaine.

Bords d'image

Sur les bords est de l'image, les ombres ont pu être projetées par un nuage sournois qui se cache en dehors de l'image. Jusqu'ici, à l'intérieur de cette zone, nous nous contentions donc de la détection de l'assombrissement, et pour éviter de nombreuses fausses détections, il fallait un fort assombrissement pour détecter une ombre. Ce qui provoquait donc de fréquentes omissions. Ici aussi, Mireille a apporté une élégante amélioration. Puisqu'il est probable que l'on ait détecté des ombres dans le reste de l'image par la méthode exposée ci-dessus, on peut mesurer leur assombrissement et s'inspirer de cette valeur pour définir le seuil de détection.

L'ombre de nuages dans le coin Nord ouest, qui n'était pas détectée dans la version précédente, l'est maintenant. Par contre en raison d'un bug détecté trop tard (la production était déjà lancée), il arrive que les ombres du dernier nuage, en parcourant l'image de gauche à droite puis de haut en bas, soient oublié.  C'est le cas sur cette image.

Généralisation

Cette amélioration appliquée dans notre prototype, qui sert à générer les images Take5, a également été mise en service pour LANDSAT 8 pour les images produites sur la France depuis le début de l'année, mais nour n'allons pas retraiter les anciennes données dans l'immédiat. L'amélioration est en train d'être reportée dans la version 5.0 de MACCS opérationnel et sera donc disponible pour traiter Sentinel-2.

 

Cas particulier

Nous n'avons rien pu faire pour le cas particulier ci dessous. L'ombre qui semble présente sur cette image est en fait un nuage de fumée noire, suite à l'explosion d'un dépôt de carburant !

Les grandeurs radiométriques : éclairement, luminance, réflectance

=>

Attention, cet article contient beaucoup de formules.

Les grandeurs radiométriques sont nombreuses et au premier abord, il n'est pas facile de s'y retrouver. Voici un petit guide des différentes grandeurs que rencontrera tout utilisateur de données de télédétection. C'est un peu rébarbatif, et c'est bourré de formules compliquées, avec des intégrales et tout ;-) , mais bon, il faut en passer par là pour bien comprendre nos images. Allez courage !

Radiométrie

Les détecteurs des instruments optiques sont sensibles à l’énergie lumineuse qu’ils reçoivent pendant un temps d’observation. Sur un appareil photo, ce temps de pose correspond à l’intervalle entre l’ouverture et la fermeture de l’obturateur. Sur un instrument d’observation numérique, on utilise davantage le terme “temps d’intégration” : celui-ci est souvent fixé par l’électronique de l’instrument.

 

L’énergie reçue est reliée aux propriétés de la scène observée, mais dépend aussi de nombreux autres paramètres liés à l’instrument lui même, ce qui n’est pas pratique lorsqu’on cherche, par exemple, à comparer des mesures prises par des instruments différents. Nous allons définir successivement, dans les paragraphes qui viennent, des grandeurs physiques qui permettent peu à peu de s’affranchir de ces paramètres dépendant de l’instrument pour obtenir une caractéristique de la surface observée seulement. Ceci va nécessiter de définir toute une série de grandeurs physiques : la description de cette série de grandeurs peut être ressentie comme rébarbative par les étudiants (et même par leur professeurs...), mais s’avère finalement bien utile un jour ou l’autre pour comprendre ce qui est observé.

Energie, Energie spectrale, Sensibilité spectrale, Energie Equivalente

Un détecteur est sensible à une énergie ε, exprimée en joules (j), reçue au cours du temps d’intégration t i. Cette énergie est apportée par des photons de différentes longueurs d’onde. En fonction de la scène observée, l’énergie spectrale ε(λ) reçue dans chaque longueur d’onde varie. L’énergie spectrale est exprimée en joules par unité de longueur d’onde (souvent en j/µm). Il faut, de plus, tenir compte du fait que l’instrument présente une sensibilité spectrale différente S(λ) en fonction des longueurs d’ondes. S(λ) est une grandeur sans unité.

Sensibilités spectrales des 5 bandes du satellite Formosat-2

 

L’énergie reçue pendant le temps d’intégration par un instrument dont la sensibilité est non nulle dans l’intervalle de longueurs d’onde [λ1, λ2] s’exprime donc de la manière suivante :

\epsilon=\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\epsilon(\lambda)\, d\lambda

Pour un instrument donné, plus la bande spectrale [λ1, λ2] est large, plus l’énergie collectée va être grande. Pour éviter de devoir manipuler des intégrales dans les formules, il est utile de définir l’énergie spectrale équivalente comme la moyenne de l’énergie spectrale dans la bande [λ1, λ2], pondérée par la sensibilité spectrale.

\epsilon_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\epsilon(\lambda)\, d\lambda}{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda}

L’avantage de cette unité est que l’énergie spectrale équivalente ne dépend plus de la largeur de bande si l’énergie spectrale ne varie pas dans la bande.

Le flux

Pour une longueur d’onde donnée, l’énergie spectrale est l’intégrale de la puissance lumineuse qui parvient sur le détecteur pendant le temps d’intégration, cette puissance lumineuse est appelée le Flux Spectral et s’exprime en Watts par unité de longueur d’onde (souvent en W/µm).

\phi(\lambda)=\frac{d\epsilon(\lambda)}{dt}

 

On peut, comme pour l’énergie, définir le flux spectral équivalent, de la manière suivante :

\phi_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\phi(\lambda)\, d\lambda}{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda}

et il suffit de développer cette formule pour vérifier que les grandeurs équivalentes sont reliées de la même manière que les grandeurs spectrales :

\phi_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{d\epsilon_{eq\,[\lambda_{1},\lambda_{2}]}}{dt}\label{eq:flux-energie-equivalente}

 

Le flux qui parvient à l’instrument reste cependant une grandeur bien éloignée des propriétés optiques de la scène que nous cherchons à mesurer. Il dépend par exemple de la surface du détecteur. Pour s’affranchir de ce paramètre, il est nécessaire d’introduire une nouvelle grandeur : l’éclairement.

L’éclairement

L’éclairement (spectral) est le flux (spectral) qui parvient au détecteur par unité de surface ds du détecteur, il se définit donc par :

E(\lambda)=\frac{d\phi(\lambda)}{ds}

 

De même, l’éclairement solaire équivalent se définit par :

E_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)E(\lambda)\, d\lambda}{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda}=\frac{d\phi_{eq\,[\lambda_{1},\lambda_{2}]}}{ds}

 

L’éclairement permet par exemple de caractériser la puissance lumineuse qui parvient sur une surface perpendiculaire à la source lumineuse, par unité de surface. L’éclairement solaire reçu sur terre varie avec la distance terre soleil, qui varie au cours de l’année. L’éclairement solaire intégré au jour j pour une surface perpendiculaire au soleil vaut

E(j)=E_{0}\frac{d^{2}(j)}{d_{0}^{^{2}}} , avec E0 = 1367 w / m2 et \frac{d^{2}(j)}{d_{0}^{2}}=(1-0.01673*cos(2\pi/365.3*(j-j0-2))) , et j0 est le 01/01/1950. L’éclairement solaire spectral présente de fortes variations avec la longueur d’onde.

Eclairement solaire en haut de l'atmosphère en fonction de la longueur d'onde. A ma connaissance, la mesure la plus précise disponible de nos jours est celle de Thuillier et al.

 

L’éclairement est une grandeur bien utile, mais qui dépend quand même du champ de vue du détecteur. Plus celui-ci sera réduit, plus l’éclairement sera réduit. Il dépend également de l’orientation du détecteur par rapport à la direction d’où provient la lumière. Par exemple, si le détecteur est orienté parallèlement à la direction des photons, il ne recevra pas de lumière. Inversement, lorsque le détecteur est perpendiculaire à la direction lumineuse, il reçoit un éclairement maximal. Pour obtenir une unité indépendante des caractéristiques de l’instrument, nous allons donc devoir définir une nouvelle grandeur (oui, encore une...).

La Luminance

La Luminance spectrale est le flux spectral qui parvient à l’instrument par unité de surface ds et par unité d’angle solide dΩ, et ce perpendiculairement à la surface du détecteur.

Définition de l’angle solide d\Omega_{s}=\frac{dS_{s}cos\theta_{s}}{r^{2}} . L’angle solide s’exprime en stéradians (sr)

La luminance s’exprime en w / m2 / sr

L(\lambda)=\frac{d^{2}\phi(\lambda)}{ds\, d\Omega\, cos\,\Theta_{s}}

 

La grandeur d2G = ds. dΩ . cos Θ s est très utilisée par les opticiens et s’appelle l’étendue géométrique.

 

Le flux spectral obtenu en entrée d’un instrument de surface S et dont le champ couvre l’angle solide Ω s’écrit donc :

\phi(\lambda)=\intop_{S}\intop_{\Omega}L(\lambda)cos\,\Theta_{s}ds\, d\Omega\,=\iint_{G}L(\lambda)d^{2}G

 

Comme pour les autres grandeurs, on définit la luminance spectrale équivalente de la manière suivante :

L_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)E(\lambda)\, d\lambda}{\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda}=\frac{d^{2}\phi_{eq\,[\lambda_{1},\lambda_{2}]}}{ds.d\Omega.cos\,\Theta_{s}}=\frac{d^{2}\phi_{eq\,[\lambda_{1},\lambda_{2}]}}{d^{2}G}

 

Lorsque l’étendue géométrique est petite, ce qui est en général le cas en télédétection, on peut écrire :

L_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{\phi_{eq\,[\lambda_{1},\lambda_{2}]}}{G}\label{eq:flux-luminance-=0000E9tendue}

 

La luminance est une unité très utilisée en télédétection, puisqu’elle combine plusieurs avantages :

  • l’énergie mesurée par un détecteur est proportionnelle à la luminance équivalente. En effet, pour un instrument d’étendue géométrique G, de sensibilité spectrale S(λ) dans la bande [λ1, λ2], qui effectue des mesures pendant le temps d’intégration ti, la mesure X obtenue par le détecteur de l’instrument s’exprime de la manière suivante :

X=\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\epsilon(\lambda)\, d\lambda=\left(\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda\right).\epsilon_{eq\,[\lambda_{1},\lambda_{2}]}.

 

Et en utilisant les formules [eq:flux-energie-equivalente] et [eq:flux-luminance-=0000E9tendue], on obtient donc

X=\left(\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda\right).\phi_{eq\,[\lambda_{1},\lambda_{2}]}.t_{i}=\left(\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda\right).L_{eq\,[\lambda_{1},\lambda_{2}]}.G.t_{i}=S{}_{[\lambda_{1},\lambda_{2}]}.L_{eq\,[\lambda_{1},\lambda_{2}]}.G.t_{i}

avec S_{_{[\lambda_{1},\lambda_{2}]}}=\intop_{\lambda_{1}}^{\lambda_{2}}S(\lambda)\, d\lambda , sensibilité du détecteur dans la bande [λ1, λ2].

 

Ce qui revient donc à écrire que :

X=A.L_{eq\,[\lambda_{1},\lambda_{2}]} .

où A est un coefficient constant pour un instrument donné, qu'on appelle coefficient d'étalonnage absolu.

  • la luminance ne dépend plus des caractéristiques instrumentales.

 

Dans cette équation, tous les éléments sont des caractéristiques de l’instrument, exceptée la luminance spectrale équivalente, qui elle est indépendante des caractéristiques de l’instrument. Les mesures obtenues sont donc proportionnelles à la luminance du pixel observé, à un facteur(S.G.t) près qui ne dépend que de l’instrument. Les concepteurs de l’instrument vont donc chercher à faire en sorte que S.G.t soit :

  • le plus grand possible, pour augmenter la sensibilité et le rapport signal sur bruit
  • le plus uniforme possible d’un pixel à l’autre, pour éviter de faire apparaître dans l’image des variations qui ne proviendraient pas de la scène observée.

La réflectance

Pour les longueurs d’ondes inférieures à 3 µm, la luminance des surfaces terrestres provient essentiellement de la réflexion par la terre de l’éclairement solaire (à l’exception de quelques sources lumineuses très intenses : projecteurs orientés vers le satellite, coulées de lave (dans le moyen infra rouge)...). Afin de s’abstraire de l’éclairement solaire, qui dépend de la distance entre la terre et le soleil (variable avec la saison), l’introduction d’une nouvelle unité est nécessaire : il s’agit de la réflectance (c’est la dernière grandeur à définir, nous sommes au bout du calvaire).

 

La réflectance spectrale se définit par :

\rho(\lambda)=\frac{\pi.L(\lambda)}{Es(\lambda,t).cos(\theta_{s})}\label{eq:def reflectance spectrale}

 

On définit par ailleurs de la même manière la réflectance équivalente :

\rho{}_{eq\,[\lambda_{1},\lambda_{2}]}=\frac{\pi.L{}_{eq\,[\lambda_{1},\lambda_{2}]}}{Es{}_{eq\,[\lambda_{1},\lambda_{2}]}(t).cos(\theta_{s})}

 

La formulation ci-dessus peut surprendre au premier abord, on aurait pu s’attendre à une formule beaucoup plus simple, la réflectance étant égale à la luminance divisée par l’éclairement. La formule adoptée pour définir la réflectance provient de considérations géométriques. Le flux reçu sur une surface dépend de son orientation par rapport à la source lumineuse. Si la surface est perpendiculaire aux rayons lumineux, le flux est maximum, alors que si elle est parallèle, il devient nul. Ceci explique la présence du cos(θs). Le facteur de normalisation π fait en sorte que la réflectance soit égale à 1 si la surface réfléchit vers le ciel la totalité de l’éclairement incident. Pour nous en convaincre, il suffit de calculer l’intégrale de la luminance réfléchie dans toutes les directions par une surface dont la réflectance est égale à 1 ( sans oublier que pour intégrer sur la sphère, il faut utiliser le facteur sin \theta_{s} d\theta_{s} d\phi et non pas seulement  d\theta_{s} d\phi .

\intop_{\theta=0}^{\frac{\pi}{2}}\intop_{\phi=0}^{2\pi}Lsin\theta_{s}d\theta_{s}d\varphi=\intop_{\theta=0}^{\frac{\pi}{2}}\intop_{\phi=0}^{2\pi}\frac{Es.cos\theta_{s}}{\pi}sin\theta_{s}d\theta_{s}d\varphi=\frac{E_{s}}{\pi}\intop_{\theta=0}^{\frac{\pi}{2}}\intop_{\phi=0}^{2\pi}\cdot cos\theta_{s}.sin\theta_{s}.d\theta_{s}d\varphi=E_{s}

puisque :
\intop_{\theta=0}^{\frac{\pi}{2}}\intop_{\phi=0}^{2\pi} cos\theta_{s}sin\theta_{s}d\theta_{s}d\varphi

 

=2\pi.\intop_{\theta=0}^{\frac{\pi}{2}}sin\theta_{s}d(sin\theta_{s})=2\pi\frac{sin^{2}\frac{\pi}{2}}{2}=\pi

 

La réflectance est une caractéristique de la surface et ne dépend ni de l’instrument qui l’observe ni de l’éclairement reçu. Cette réflectance dépend cependant de l’angle sous lequel la surface est éclairée et observée et aussi, bien sûr, de la longueur d’onde. On entend parfois parler de l'albédo spectral qui est la moyenne de la réflectance dans toutes les directions, ou de l'albédo, qui est la moyenne de l'albédo spectral sur tout le spectre.

a(\lambda)=\frac{1}{\pi}\intop_{\theta=0}^{\pi}\intop_{\phi=0}^{2\pi}\rho(\lambda).d\theta_{s}d\varphi

 

Cette fois, nous sommes au bout de nos peines, et nous avons démontré que les mesures de nos satellites sont proportionnelles à une grandeur caractéristique de la surface, la réflectance. C'est en réflectance que les données Sentinel-2, ou LANDSAT sont exprimées, même s'il nous a fallu discuter avec les tenants de la luminance. La réflectance présente en effet de nombreux avantages :

  • La réflectance est une caractéristique de la surface, elle ne dépend que de la surface et des angles d'éclairement et d'observation (cf les effets directionnels)
  • Ses valeurs sont en général comprises entre zéro et un (même si des valeurs supérieures à  1 peuvent être observées, c'est assez rare). Ces valeurs sont faciles à mémoriser, contrairement aux valeurs de luminances, qui varient d'une bande à l'autre en fonction de l'éclairement solaire
  • Il est possible de comparer les réflectances d'une bande spectrale à l'autre.
  • Il est possible de comparer directement les réflectances mesurées en hiver et en été, sans avoir à prendre en compte la variation de l'élévation solaire.

Réflectances d'une parcelle de blé, dans les 4 bandes (bleu, vert rouge, et proche infra-rouge), en fonction du temps. La même courbe exprimée en luminance serait beaucoup moins intéressante.

 

Let it snow! Development of an operational snow cover product from Sentinel-2 and Landsat-8 data

=>

 

"Winter is coming" ― George R.R. Martin, A Game of Thrones

 

As Christmas holidays are approaching you might want to know if there is snow in your favorite spot of ski touring? A good knowledge of the snow cover variability is important - not only to plan your next week-end, but also because the snow is a key water resource in many regions, including here in south west France. Continue reading

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

 

The Sentinel-2 tiles, how they work ?

=> 

The Sentinel-2 Level 1C products will be split into tiles. and I have searched a little how the naming convention of Sentinel-2 tiles works. I have found some interesting resources in Sentinel-2 website :

  • A User Handbook for Sentinel-2 products (which does not explain the naming convention of the tiles)
  • A kml file that provides the footprint of all the tiles, and their name.

 

 

 

Here is  what I have guessed :

  • The first 2 numbers of a tile name (such as 31TCJ) correspond to the UTM zone. The world is divided in 60 UTM zones of 8 degrees width in longitude, with numbers increasing towards the East. Zone 1 is over the Pacific Ocean.
  • Each zone is divided in latitude, by chunks of 6 degrees. This is represented by a letter, which increases from South to North.
  • And finally, each chunk is divided in 110 km tiles, with a 10 km overlap, from West to East, 4th Letter, and South to North, 5th Letter. What is not easy to guess is the fact that the forth and fifth letters are not reset for each chunk, but continue to increase when changing zone and chunk, but are reset to A when the counter exceeds Z.

 

For instance :

  • the tile in Toulouse, id 31TCJ. 31 is the UTM zone, T is the latitudinal chunk. C denotes the West-East tile position within the chunk and J the South-North position.
  • If you want to see the tile in the East of Toulouse, which includes Castres, it will be 31TDJ.
  • If you want to see the tile in the North East of Toulouse, which includes Aurillac, it will be 31TDK.
  • If you want to see the tile at the West of Toulouse, which includes the Armagnac region, it will be... 30TYP. 30T because it is a different UTM zone, and YP ... because this numbering is not practical  ?
  • The tile at the South of that one, for which we processed our first Level 2A products is  30TYN (because the O letter is excluded, to avoid confusion with 0)

 

Sorry for that, I would like to have a formula which guesses the tile name as a function of the longitude and latitude, but I am not able to figure out one. Has any of you heard of one ?

 

What I am happy with is the fact that tiles overlap by 10km. This ovelap means an increase of 21 % of the data volume, which is already large, but enables to process Level 2A products independently. It is very difficult to detect a cloud shadow when you do not see the corresponding cloud. And on the edge of the tiles, you will have shadows on clouds you do not see. But thanks to the overlap, the shadows will be well detected on the other overlapping tile.

 

A colleague of mine told me this grid was born in the brain of an American solldier (yes, soldiers have brains ). Here is a description, which sadly does not tell how the two last letters were defined.