Algorithm to reconstruct a snow cover timeseries with high spatiotemporal resolution

Original reference: https://ieeexplore.ieee.org/document/9511108

Data Sources

High spatial, low temporal resolution: Landsat and Sentinel-2

• Sentinel-2 SCA can be obtained from SCL, but better is to use ‘SNOWFLAKE’ algorithm

Low spatial, high temporal resolution: MODIS snow cover fraction:

• https://stac.eurac.edu/collections/MOD10A1v61
• https://radiantearth.github.io/stac-browser/#/external/stac.openeo.vito.be/collections/modis-10A1-061

Long timeseries needed to have sufficient input for statistical methods.

Step 1: High resolution gap filling

For each partially clouded high-resolution image, compare it to every other HR observation in the timeseries and compute similarity. Fill in gaps if similarity is above a certain threshold.

Original code: https://github.com/vpremier/dailyHR_SCA/blob/07f0bc37eb55a1dd13606b325ce801790c2cc5f6/gapfilling.py#L40

Proposed openEO implementation

We need the full scene for spatial similarity metric and we need the full timeseries. So a callback using apply_neighborhood seems in order.

Using a UDF will allow us to start from the original code, maintaining similar performance. Biggest drawback is higher memory use, but other methods seem to imply a high degree of data movement, making it hard to achieve similar performance.

Step 2: Downscaling of low-resolution time series to high resolution

For each high-resolution pixel, the probability distribution of snow cover based on the low-resolution snow cover fraction is computed. So observed high-resolution SCF is first resampled to low-resolution, obtaining LR SCF. Then for all low resolution pixels with a similar SCF, the snow presence probablity is computed.

The distribution computation can be done once, and only changes when sufficient new data is available.

Original code for computing distribution:: https://github.com/vpremier/dailyHR_SCA/blob/07f0bc37eb55a1dd13606b325ce801790c2cc5f6/downscaling.py#L115

Original code to fill in high-resolution gaps based on computed distribution: https://github.com/vpremier/dailyHR_SCA/blob/07f0bc37eb55a1dd13606b325ce801790c2cc5f6/historical_block.py#L97

Proposed openEO implementation

Distribution computation already exists: https://github.com/SNOWCOP/openeo_mountains_snow/blob/9df6ed19dd5904cd82d764867daeb92c8249af3d/src/openeo_mountains_snow/snowcoverarea_reconstruction/downscaling_distribution.py#L174-L173

Actual downscaling is a matter of merging the distribution cube with MODIS time series upsampled to 20m, apply mask for gaps to be filled, and then filling the actual gaps.

Step 3: Improving LR time series and running again

Original paper runs the ‘historic block’ two times.

Step 4: run pyGAM to fill remaining gaps

This step fills remainging gaps. It can also be avoided by simply relaxing the threshold for step 2.

Run inference based on input features:

• NDVI
• NDSI
• CP
• SCA

https://github.com/vpremier/dailyHR_SCA/blob/c0728dfdb9d1a9b46b9ff33a710d178dab7a1baa/gam.py#L484