Towards a composite of the most recent clear-sky observations from the Theia Snow collection

We are thinking of a new « live » product which would allow anyone to display the most recent snow observations from Sentinel-2 in a web browser. This service would simply display a composite of the latest Theia Snow products over a given time window (e.g. the past 20 days).

This can be easily coded thanks to the Orfeo toolbox (OTB) band math application. My colleagues at CNES taught me that it is preferable to use OTB in operational context because it enables to specify the available memory. A typical Python code would load the full images into memory to process them as a stack of numpy arrays. Here we need to process an unknown number of snow products to compute the composite, because the number of acquisitions over a time window depends on the Sentinel-2 tile (revisit time is 5 day or less!). Therefore we cannot easily anticipate the required memory to do the processing with numpy arrays. Should the memory usage of our process exceed the allocated resources, it would be instantaneously killed by the merciless job scheduler of the CNES HPC…

The Theia snow products are coded as follows:

• 0: No-snow
• 100: Snow
• 205: Cloud including cloud shadow
• 255 (or 254 before LIS 1.6): No data

If we have a list of snow products in reversed chronological order (im1 is the most recent, then im2, etc.), then we can define the compositing expression using if-then-else operator: (condition ? value_true : value_false)

Two products:

im1b1 <= 100 ? im1b1 : im2b1

where im1b1 means « image 1 band 1 ».

In other words: if the value is snow or no-snow in the most recent product, then keep it, otherwise take the value of next product.

Three products:

im1b1 <= 100 ? im1b1 : (im2b1 <= 100 ? im2b1 : im3b1)

Four products:

im1b1 <= 100 ? im1b1 : (im2b1 <= 100 ? im2b1 : (im3b1 <= 100 ? im3b1 : im4b1))

etc..

We can even get rid of the parentheses, which greatly simplifies the code to build the expression for any number of products.

Here is an example using my favorite tile in the Pyrenees (snow: cyan, no-snow: grey, clouds: white, no-data: black)

The above images are all the available products from Theia since May 1st (as of yesterday) obtained with:

./theia_download.py -t 'T31TCH' -c Snow -a config_theia.cfg -d 2020-05-01 -f 2020-05-20

The composite of the most recent clear-sky observations is:

Because the above image was actually generated at 20 m resolution we can use it to check the snow conditions on hiking trails, or on the highways..

We are still in the development phase, let us know in the comments below if you think this product would be useful and if you have suggestions to improve it!