# Theia Snow collection: High-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data

Simon Gascoin, M. Grizonnet, M. Bouchet, G. Salgues, Olivier Hagolle

► To cite this version:

Simon Gascoin, M. Grizonnet, M. Bouchet, G. Salgues, Olivier Hagolle. Theia Snow collection: High-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data. Earth System Science Data, 2019, 11 (2), pp.493–514. 10.5194/essd-11-493-2019. hal-02351620

HAL Id: hal-02351620

https://hal.science/hal-02351620

Submitted on 25 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Theia Snow collection: high-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data

Simon Gascoin¹, Manuel Grizonnet², Marine Bouchet¹, Germain Salgues²,³, and Olivier Hagolle¹,²

¹CESBIO, Université de Toulouse, CNRS/CNES/IRD/INRA/UPS, Toulouse, France

²CNES, Toulouse, France

³Magellium, Toulouse, France

**Correspondence:** Simon Gascoin (simon.gascoin@cesbio.cnes.fr)

Received: 22 November 2018 – Discussion started: 29 November 2018

Revised: 16 March 2019 – Accepted: 20 March 2019 – Published: 16 April 2019

**Abstract.** The Theia Snow collection routinely provides high-resolution maps of the snow-covered area from Sentinel-2 and Landsat-8 observations. The collection covers selected areas worldwide, including the main mountain regions in western Europe (e.g. Alps, Pyrenees) and the High Atlas in Morocco. Each product of the Theia Snow collection contains four classes: snow, no snow, cloud and no data. We present the algorithm to generate the snow products and provide an evaluation of the accuracy of Sentinel-2 snow products using in situ snow depth measurements, higher-resolution snow maps and visual control. The results suggest that the snow is accurately detected in the Theia snow collection and that the snow detection is more accurate than the Sen2Cor outputs (ESA level 2 product). An issue that should be addressed in a future release is the occurrence of false snow detection in some large clouds. The snow maps are currently produced and freely distributed on average 5 d after the image acquisition as raster and vector files via the Theia portal (https://doi.org/10.24400/329360/F7Q52MNK).

## 1 Introduction

The snow cover is an important driver of many ecological, climatic and hydrological processes in mountain regions and in high latitude areas. In these regions, in situ observations are generally insufficient to characterize the high spatial variability of the snowpack properties. The snow cover was included as one of the 50 essential climate variables (ECVs) to be monitored by satellite remote sensing by the Global Observing System for Climate (GCOS) in accordance with the Committee on Earth Observation Satellites (CEOS) agencies. ECVs are intended to support the work of the United Nations Framework Convention on Climate Change and the Intergovernmental Panel on Climate Change.

The snow cover is not a variable sensu stricto, but an object which can be characterized through many variables, including snow-covered area (SCA), fractional area (fSCA), albedo, liquid water content, snow depth and snow water equivalent (Frei et al., 2012). An international survey by the Cryoland consortium (http://cryoland.enveo.at/consortium,

last access: 12 April 2019) shed light on the user requirements for snow services based on satellite remote sensing. SCA and fSCA products were ranked as the most important by the respondents (Malnes et al., 2015). The survey also indicated the need for operational products provided on a daily basis, with latency times shorter than 12 h. High-resolution data down to 50 m resolution were sought by road and avalanches authorities (Malnes et al., 2015). The respondents requested regional products, e.g. on the scale of entire mountain ranges like the Alps or even the whole of Europe, and preferred the Universal Transverse Mercator (UTM) projection (Malnes et al., 2015). At the national level in France there is also a need for operational high-resolution snow-covered area maps as revealed by the recent roadmap for satellite applications issued by the French Government (*Plan d'applications satellitaires 2018*, Commissariat général au développement durable, 2018). Based on a wide panel of end users, this plan selected the monitoring of the snow-covered area in French national parks as a one of the key actions
which should be achieved using Earth observation satellites
in the near future (2018–2022).

Operational SCA maps have been generated from satellite
observations since the 1960s (Ramsay, 1998). Current SCA
and fSCA products are mostly derived from MODIS data
(Hall et al., 2002; Sirguey et al., 2008; Painter et al., 2009;
Metsämäki et al., 2015), but their spatial resolution (1 km to
250 m) can be too coarse for various applications, especially
in mountain regions where the snow cover properties vary on
scales of 10 to 100 m (Blöschl, 1999). High-resolution (30 m)
snow cover maps can be generated from Landsat images but
the low temporal revisit of the Landsat mission (16 d) is an
important limitation for snow cover monitoring, especially
considering that the cloud probability can exceed 50 % in
temperate mountains (Parajka and Blöschl, 2008; Gascoin
et al., 2015). Since 2017, with the launch of Sentinel-2B, the
Copernicus Sentinel-2 mission offers the unique opportunity
to map the snow cover extent at 20 m resolution with a re-
visit time of 5 d (cloud permitting) (Drusch et al., 2012). The
combination of Sentinel-2 and Landsat-8 data provides the
opportunity for even more frequent observations of the snow
cover with a global median average revisit interval of 2.9 d
(Li and Roy, 2017).

The principles of snow detection from multispectral op-
tical imagery is well established since the pioneering work
of Dozier (1989) with the Landsat Thematic Mapper. Today,
the challenge for scientists and end users is rather to cope
with the large amount of data that is generated by a mission
like Sentinel-2. The generation of a single snow map from
a Sentinel-2 level-1C tiled product (monodate orthorectified
image expressed in top-of-atmosphere reflectance) involves
downloading a 700+ Mb zip file (once uncompressed, the
product contains 12 folders and 108 files, including 15 raster
files and 13 XML metadata files). Since March 2018, ESA
began the operational processing of level 2A (L2A) prod-
ucts (monodate orthorectified images expressed in surface
reflectance, including a cloud mask). Each ESA L2A prod-
uct also includes a snow cover mask. However, the size of
a single L2A product before unzipping can exceed 1 Gb (15
folders, 137 files). In addition, the quality of the ESA L2A
snow mask can be improved for two reasons: (i) ESA's L2A
algorithm is a general-purpose algorithm, which was not op-
timized for snow detection, and (ii) because ESA's L2A pro-
cessor treats each image independently (i.e. mono-date ap-
proach, Louis et al., 2016), the output cloud mask has a lower
accuracy than a cloud mask generated by a multi-temporal
algorithm (Hagolle et al., 2010).

In this article we introduce the Theia Snow collection, a
new collection of snow cover maps, which are derived from
Sentinel-2 at 20 m resolution in an operational context. The
Theia Snow collection was recently upgraded to also provide
snow cover maps at 30 m resolution from Landsat-8 using
the same algorithm. The data are currently being produced

and freely distributed via the Theia land data centre.<fnref n="1" /> The
snow retrieval algorithm is based on the Normalized Differ-
ence Snow Index (Dozier, 1989) but also uses a digital ele-
vation model to better constrain the snow detection. We first
describe the algorithm (Sect. 2) and its implementation in the
let-it-snow processor (LIS, Sect. 2.4.4). Then, we provide a
detailed description of the product characteristics (coverage,
period, format, etc.; Sect. 3). In Sect. 4, we present an evalu-
ation of Theia snow products using in situ snow depth mea-
surements, very high-resolution clear-sky snow maps and vi-
sual control. We finally discuss the main limitations and po-
tential applications of the Theia Snow collection.

2 Algorithm

2.1 Scope

An algorithm was designed to determine the snow presence
or absence from Sentinel-2 and Landsat-8 observations out-
side areas of dense forest with the following requirements:

- It should be scalable, i.e. allow the processing of large areas (10⁴ km²) with a reasonable computation cost (typically less than 1 h for a single product).

- It should be robust to seasonal and spatial variability of the snow cover and land surface properties.

- It should maximize the number of pixels that are classified as snow or no snow.

- It is preferable to falsely classify a pixel as cloud than falsely classify a pixel as snow or no snow.

2.2 Input

The algorithm works with multispectral remotely sensed im-
ages, which include at least a channel in the visible part
of the spectrum and a channel near 1.5 µm (referred to as
shortwave-infrared or SWIR). It takes the following as input:

- a L2A product, including

- the cloud and cloud shadow mask (referred to as an
  “L2A cloud mask” in the following),

- the green, red and SWIR bands from the flat-surface
reflectance product. These images are corrected
for atmospheric and terrain slope effects (Hagolle
et al., 2017). The slope correction is important in
mountain regions since it enables us to use the same
detection thresholds whatever the sun-slope geom-
etry.

- A digital elevation model (DEM).

<footnote n="1">¹Theia is a French inter-agency organization designed to foster the use of Earth observation for land surface monitoring for academics and public policy actors.</footnote>
## 2.3 Pre-processing

The red and green bands are resampled with the cubic method to a pixel size of 20 m by 20 m to match the resolution of the SWIR band. The DEM is generated from the Shuttle Radar Topography Mission seamless DEM (Jarvis et al., 2008) by cubic spline resampling to the same 20 m resolution grid.

## 2.4 Algorithm description

### 2.4.1 Snow detection

The snow detection is based on the Normalized Difference Snow Index (NDSI, Dozier, 1989) and the reflectance in the red band. The NDSI is defined as

$$NDSI = \frac{\rho_{\text{green}} - \rho_{\text{SWIR}}}{\rho_{\text{green}} + \rho_{\text{SWIR}}}, \qquad (1)$$

where $\rho_{\text{green}}$ (or $\rho_{\text{SWIR}}$) is the slope-corrected surface reflectance in the green band (or SWIR at 1.6 µm). The NDSI expresses the fact that only snow surfaces are very bright in the visible and very dark in the shortwave infrared. Turbid water surfaces like some lakes or rivers may also have a high NDSI value; hence a additional criterion on the red reflectance is used to avoid false snow detection in these areas. A cloud-free pixel is classified as snow if the following condition is true:

$$(\text{NDSI} > n_i) \text{ and } (\rho_{\text{red}} > r_i), \qquad (2)$$

where $n_i$ and $r_i$ are two parameter pairs (i.e. $i \in \{1, 2\}$) as explained below. Otherwise the pixel is marked as “no snow”. The values of the parameters are provided in Table 1.

### 2.4.2 Snow line elevation

The algorithm works in two passes. The snow detection (Sect. 2.4.1) is performed a first time using parameters $n_1$ and $r_1$, which are set in such way as to minimize the number of false snow detections (Table 1). This first pass enables us to estimate the minimum elevation of the snow cover $z_s$, above which a second pass will be performed with less conservative parameters values $n_2$ and $r_2$ (Table 1). An example of the evolution of the snow line elevation over the tile 31TCH in the Pyrenees is shown in Fig. 1 and the corresponding snow maps are shown in Fig. 1. Based on regional analyses of the snow line elevation variability in mountain ranges, (e.g. Gascoin et al., 2015; Krajčí et al., 2016), we assume that large altitudinal variations of the snow line elevation are not likely on the scale of a tile of 110 km by 110 km. Therefore, the minimum snow elevation $z_s$ is considered uniform on the scale of the processed image.

Before proceeding to pass 2, the total snow fraction in the image after pass 1 is computed. If this snow fraction is below $f_t$ (Table 1), then pass 2 is skipped, as the sample of snow

pixels is not considered statistically significant for determining the snow line elevation. Otherwise, the pass 2 is activated. For that purpose, the DEM is used to segment the image in elevation bands with a fixed height of $d_z$ (Table 1). Then, the fraction of the cloud-free area that is covered by snow in each band (after pass 1) is computed. The algorithm finds the lowest elevation band $b$ at which the snow cover fraction is greater than a given fraction $f_s$ (Table 1). The value of $z_s$ is the lower edge of the elevation band that is two elevation bands below band $b$.

### 2.4.3 Cloud mask processing

The cloud mask in the input L2A product is conservative because (i) it is computed at a coarser resolution and (ii) it was developed to remove surface reflectance variations due to cloud contamination. However, the scattering of some thin clouds is low in the SWIR, green and red bands, which are used for snow detection (Sect. 2.4.1). Hence, the human eye can see the snow (or the absence thereof) through these semi-transparent clouds in a colour composite. In addition, the L2A cloud mask tends to falsely classify the edges of the snow cover as cloud. Therefore, the algorithm includes some additional steps to recover those pixels from the L2A cloud mask and reclassify them as snow or no snow. This step is important because it substantially increases the number of observations as specified in Sect. 2.1.

A pixel from the L2A cloud mask cannot be reclassified as snow or no snow if any of these conditions are satisfied:

- it is coded as “cloud shadow” in the L2A cloud mask;

- it is coded as “high-altitude cloud” (or “cirrus”) in the L2A cloud mask;

- it is not a “dark cloud” (see below).

The cloud shadows are excluded because the signal-to-noise ratio can be very low in these areas. The high clouds are excluded because they can have a similar spectral signature as the snow cover, i.e. a high reflectance in the visible and a low reflectance in the SWIR. This type of cloud are detected in Sentinel-2 and Landsat-8 L2A products based on the spectral band centred on the 1.38 µm wavelength (Gao et al., 1993).

We select only the dark clouds as possible reclassification areas, because the NDSI test is robust to the snow-cloud confusion in this case. The dark clouds are defined using a threshold in the red band after downsampling the red band by a factor $r_f$ using the bilinear method. This resampling is applied to smooth locally anomalous pixels, following the MAJA algorithm, which performs the cloud detection at 240 m for L2A products (Hagolle et al., 2017). Therefore, if a (non-shadow, non-high-cloud) cloud pixel has a red reflectance at this coarser resolution that is lower than $r_D$ (Table 1), then it is temporarily removed from the cloud mask
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Time series of snow and cloud cover area distribution by elevation band over tile 31TCH (tile location in Fig. 5). The dashed red line indicates the position of zs as determined by the LIS algorithm. The corresponding snow maps are shown in Fig. 2.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time series of snow maps over tile 31TCH (tile location in Fig. 5). The corresponding charts of snow and cloud cover area distribution by elevation band are shown in Fig. 1.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Flow chart of the snow detection algorithm. Table 1 gives the description and value of the algorithm parameters (written in red in this chart). MUSCATE is the scheduler which manages the L2A production for Theia.</figcaption></figure>

For example, the eastern snow product in Fig. 5 was extracted
from a file named

SENTINEL2A_20151130-105641-486_L2B-
-SNOW_T31TDH_D_V1-0.zip,

which indicates that this product was generated from a
Sentinel-2A image acquired on 30 November 2015 at

10:56:41.486 UTC and covers tile 31TDH. The version of
the product is 1.0 (see above).

The size of each zipped product varies depending on the
complexity of the snow and cloud mask but generally ranges
between 10 and 100 Mb. After extracting the product, the fol-
lowing files are created:
<table><caption>**Table 1.** LIS algorithm parameters description and default values for Theia's L2A Sentinel-2 (and Landsat-8) products.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Value</td></tr></thead><tbody><tr><td>r<sub>f</sub></td><td>Resize factor to produce the downsampled red band</td><td>12 (8)</td></tr><tr><td>r<sub>D</sub></td><td>Maximum value of the downsampled red-band reflectance used to define a dark cloud pixel</td><td>0.300</td></tr><tr><td>n<sub>1</sub></td><td>Minimum value of the NDSI for the pass 1 snow test</td><td>0.400</td></tr><tr><td>n<sub>2</sub></td><td>Minimum value of the NDSI for the pass 2 snow test</td><td>0.150</td></tr><tr><td>r<sub>1</sub></td><td>Minimum value of the red-band reflectance for the pass 1 snow test</td><td>0.200</td></tr><tr><td>r<sub>2</sub></td><td>Minimum value of the red-band reflectance for the pass 2 snow test</td><td>0.040</td></tr><tr><td>d<sub>z</sub></td><td>Size of elevation band in the DEM used to define z<sub>s</sub></td><td>100</td></tr><tr><td>f<sub>s</sub></td><td>Minimum snow fraction in an elevation band to define z<sub>s</sub></td><td>0.100</td></tr><tr><td>fc<sub>t</sub></td><td>Minimum clear pixels fraction (snow and no snow) in an elevation band used to define z<sub>s</sub></td><td>0.100</td></tr><tr><td>f<sub>t</sub></td><td>Minimum snow fraction in the image to activate the pass 2 snow test</td><td>0.001</td></tr><tr><td>r<sub>B</sub></td><td>Minimum value of the red band reflectance to return a non-snow pixel to the cloud mask</td><td>0.100</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Coverage of the Theia Snow collection tiles.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Map of the first two products in the Theia Snow collection and comparison with MOD10A1.005 snow product of the same day. The map background is the colourized version of the EU DEM.</figcaption></figure>

<table><caption>**Table 2.** Spatial resolution and central wavelength of the spectral bands used by LIS for each compatible sensor.</caption><thead><tr><td>Sensor</td><td>Green band</td><td>Red band</td><td>SWIR band</td></tr></thead><tbody><tr><td>SPOT-4 HRV</td><td>20 m, 0.55 µm</td><td>20 m, 0.65 µm</td><td>20 m, 1.6 µm</td></tr><tr><td>SPOT-5 HRG</td><td>10 m, 0.55 µm</td><td>10 m, 0.65 µm</td><td>10 m, 1.6 µm</td></tr><tr><td>Sentinel-2 MSI</td><td>10 m, 0.56 µm</td><td>10 m, 0.66 µm</td><td>20 m, 1.6 µm</td></tr><tr><td>Landsat-8 OLI</td><td>30 m, 0.56 µm</td><td>30 m, 0.65 µm</td><td>30 m, 1.6 µm</td></tr></tbody></table>

- `productID_SNW_R2.tif` is an 8 bit single-band GeoTiff image which provides the following classification for each pixel:

- 0: no snow,

- 100: snow,

- 205: cloud including cloud shadow,

- 254: no data.

- `productID_SNW_R2.{shp,shx,dbf,prj}` is a vector file in the ESRI shapefile format,

which contains a vectorized version of the product ID `SNW_R2.tif` using the polygon geometry.

- `productID_CMP_R2.tif` is an 8 bit RGB GeoTiff image which shows the outlines of the snow mask (green lines) and cloud mask (magenta lines) on a false-colour composite of the input L2A image (RGB image made with SWIR, red and green bands). This composition was chosen because RGB composites using the SWIR-band image are useful for discriminating the snow cover from the snow-free areas and from the clouds (Vidot et al., 2017). Hence, this image mainly allows the user to visually check the consistency of the snow and cloud mask.

- `productID_MTD_ALL.xml` is a metadata file.

- `productID_QKL_ALL.jpg` is a quick picture of the product ID `SNW_R2.tif` image made using this colour map (8 bit RGB code in parentheses): snow in cyan (0,255,255), cloud in white (255,255,255), no
<table><caption>**Table 3.** Non-exhaustive list of mountain ranges covered by the Theia Snow collection and corresponding tiles. The tiles marked with an asterisk were made available before the start of the operational production (see Sect. 3).</caption><thead><tr><th>Mountain range</th><th>Tiles</th></tr></thead><tbody><tr><td>Alps</td><td>32TMS, 32TMR, 32TMT, 32TNS, 32TNR, 32TLQ*, 32TLP, 32TLS*, 32TLR*, 32TQT, 32TQS, 33TUM, 31TGK*, 32TNT, 31TGJ, 32TPR, 31TGL*, 32TPT, 32TPS</td></tr><tr><td>Appenines</td><td>33TUH, 32TNQ, 33TUG, 33TVG, 33TWF, 32TPP, 33TWE, 32TPQ, 32TQP, 32TQN</td></tr><tr><td>Cantabrian</td><td>30TUN, 29TPH, 29TPG, 29TQH</td></tr><tr><td>Corsica</td><td>32TMM, 32TMN, 32TNM, 32TNM</td></tr><tr><td>Grant Range</td><td>11SPC</td></tr><tr><td>High Atlas</td><td>29RPQ*, 29RNQ**, 29SQR*, 29SPR*, 30STA*</td></tr><tr><td>Sistema Iberico</td><td>30TVM, 30TWK, 30TWM, 30TWL, 30TXK, 30TXL</td></tr><tr><td>Jura</td><td>31TGM*, 32TLT</td></tr><tr><td>Lebanon</td><td>T36SYD, T37SBU, T36SYC, T37SBT, T36SYB</td></tr><tr><td>Massif Central</td><td>31TDL, 31TDK, 31TEK, 31TEL</td></tr><tr><td>Pyrenees</td><td>30TXN*, 30TYM, 30TYN*, 31TCG, 31TDG, 31TCH*, 31TDH*</td></tr><tr><td>Sierra de Cazorla</td><td>30SWH</td></tr><tr><td>Sierra Nevada</td><td>30SVG, 30SWG</td></tr><tr><td>Sistema Central</td><td>30TUL, 30TUK, 30TVL, 30TTK</td></tr><tr><td>Vosges</td><td>32TLT, 32ULU</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Cumulative probability of the time lags between the Sentinel-2 acquisition and the distribution date of the corresponding snow product, computed for all snow products with acquisition date in April 2018.</figcaption></figure>

snow in grey (119,119,119) and no data in black (0,0,0) (see an example in Fig. 7).

- DATA/productID_HIS_R2.txt (since product version 1.4) is a text file indicating the cloud, snow and no-snow fraction area by elevation bins (Sect. 2).

- MASKS/productID_EXS_R2.tif (in product version 1.0, it is located in the root of the zipped file) is an 8 bit GeoTiff image for expert evaluation, where each

bit indicates the value of an intermediate computation mask (Sect. 2):

- bit 1: snow (pass 1),

- bit 2: snow (pass 2),

- bit 3: clouds (pass 1),

- bit 4: clouds (pass 2),

- bit 5: clouds (initial all cloud).

For most applications, only files with the SNW suffix should be useful.

The shapefile and raster images are referenced in WGS-84 UTM system with the zone number given by the first two digits of the tile name. All raster files have a 20 m resolution. Note that no-snow pixels can be any surface, including water surface. No data pixels are the pixels outside of the acquisition segment.

## 4 Evaluation

In this section we present an evaluation of the Theia Snow collection. We first present the methods developed to do this evaluation and then the results. The evaluation was based only on Sentinel-2 products, because the production for Landsat-8 had just started at the time of writing and thus the Landsat-8 products represented a very minor fraction of the Theia Snow collection.

<footnote>* These tiles can be viewed on an interactive map on this page: https://theia.cnes.fr/atdistrib/rocket/AT1	extbackslash slash#/sites (last access: 12 April 2019).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Example of a snow map (SNW file) and its corresponding colour composite (CMP file). In the snow map the snow is represented in cyan, no snow in grey and cloud in white. In the RGB composite, the snow-covered areas are delineated in magenta, while the clouds are delineated in green (Theia snow product of tile 32TNS on 27 May 2017).</figcaption></figure>

## 4.1 Method

### 4.1.1 Comparison with in situ snow depth measurements

We collected all available snow depth measurements within tiles 31TGM, 31TGL, 31TGK, 32TLS, 32TLR, 32TLQ, 31TDH, 31TCH, 30TYN and 30TNX from the Météo-France database between 1 September 2017 and 31 August 2018. These 10 tiles cover most of the French Alps and Pyrenees (Fig. 8). We obtained 120 snow depth time series with at least one measurement. We gathered all available Sentinel-2 snow products for these tiles over the same period, i.e. a total of 1134 products. Then, we extracted the pixel values at the location of each snow measurement station for all dates. When a station is located in two tiles we only kept the data from the first tile. We selected the snow depth measurements which were collected on the same day of a Theia snow product and discarded the measurements corresponding to a cloud detection in the Theia snow product. The snow depth measurements were converted to snow presence and absence using a threshold of $SD_0 = 0$ m. This eventually allowed us to compute a confusion matrix between a set of snow presence and absence data from in situ measurement and a set of simultaneous snow presence and absence data from the Theia Snow collection across the French Alps and Pyrenees. However, previous studies comparing MODIS binary snow

products with ground measurements showed that a value of 0 m may not be optimal (Klein and Barnett, 2003; Gascoin et al., 2015); therefore the sensitivity of the results to $SD_0$ was tested by recomputing the confusion matrix for 1 cm increments of $SD_0$ from 0 to 1 m.

### 4.1.2 Comparison with snow maps of higher spatial resolution

We used SPOT-6 and SPOT-7 images with a resolution of 1.5 m in panchromatic and 6 m in multispectral (blue, green, red, near-infrared) to evaluate the accuracy of the snow detection in Theia snow products. For this comparison, the LIS processor was run with its current operational configuration (i.e. the configuration used to routinely produce the Theia snow products since December 2017 Table 1). SPOT-6 and SPOT-7 sensors acquire images with a large radiometric depth coded in 16 bits, which enables us to identify surface features in alpine regions with dark shaded slopes and bright snow surfaces. We searched the catalogue of the Kalideos database for available SPOT images that could match a cloud-free (or nearly cloud-free) Theia snow map over the French Alps region (https://alpes.kalideos.fr, last access: 12 April 2019). We identified six pairs of images in 2016 and 2017 with a maximum time lag of 6 d (Table 4, Fig. 9). The SPOT images were obtained as orthorectified
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Map showing the location of the stations in the réseau nivo-météorologique (Météo-France snow observation network) used for this study (source: Météo-France) and the corresponding Sentinel-2 tiles.</figcaption></figure>

products from the Kalideos database. The SPOT images were used to generate reference snow maps by a pixel-based supervised classification. For each SPOT image, we manually drew about 15 polygons of homogeneous of snow and no-snow surfaces using the panchromatic image. The few cloud pixels in the SPOT images were also manually delineated with a large buffer to restrict the classification to strictly cloud-free areas. Colour composites made from the multispectral images were also used as a visual support to help the snow and cloud classification. These polygons were then used to extract training samples from the SPOT products. The samples were taken from both the panchromatic and the multispectral images. We tested a number of classification algorithms by splitting the polygons in calibration and validation data sets. The detail of this analysis is not shown here, but it allowed us to conclude that the random forest classifier was the best choice for our purpose, given its good accuracy and its numerical efficiency (Bouchet, 2018).

In addition, we generated the snow maps from the same Sentinel-2 level 1 products using ESA's Sen2Cor version 2.5 (Louis et al., 2016). The Sen2Cor output snow masks are expected to be nearly identical to the ones included in the distributed L2A products by ESA, since the same processor is used. We have generated the L2A products ourselves with Sen2Cor for this study, because the L2A products were not yet available when we made this analysis (ESA started operational production of Sentinel-2 data at level 2A in 2017). The Theia and Sen2Cor snow maps were compared to the reference snow maps from SPOT images using standard met-

<table><caption>**Table 4.** Pairs of Sentinel-2 and SPOT-6/7 products used for the evaluation of the snow detection accuracy (see also Fig. 9).</caption><thead><tr><th colspan="2">Sentinel-2 product</th><th colspan="2">Reference product</th></tr><tr><th>Tile</th><th>Date</th><th>Sensor</th><th>Date</th></tr></thead><tbody><tr><td>31TGK</td><td>13 Aug 2016</td><td>SPOT-6</td><td>8 Aug 2016</td></tr><tr><td>32TLS</td><td>12 Oct 2016</td><td>SPOT-7</td><td>12 Oct 2016</td></tr><tr><td>31TGK</td><td>11 Dec 2016</td><td>SPOT-7</td><td>17 Dec 2016</td></tr><tr><td>31TGL</td><td>1 Dec 2016</td><td>SPOT-7</td><td>3 Dec 2016</td></tr><tr><td>31TGK</td><td>11 Mar 2017</td><td>SPOT-6</td><td>11 Mar 2017</td></tr><tr><td>31TGL</td><td>11 Mar 2017</td><td>SPOT-6</td><td>11 Mar 2017</td></tr></tbody></table>

rics derived from the confusion matrix (accuracy, F1 score, kappa, false-detection rate and false-negative rate).

### 4.1.3 Evaluation metrics

From the confusion matrices of Sect. 4.1.1 and 4.1.2, we derived the following metrics: accuracy (the proportion of the total number of predictions that were correct), false-positive rate (FPR, i.e. the proportion of no-snow pixels that were incorrectly classified as snow), false-negative rate (FNR, i.e. the proportion of snow pixels that were incorrectly classified as no snow), F1 score (harmonic average of the precision and recall) and kappa coefficient (Cohen, 1960).
<table><caption>**Table 5.** Confusion matrix between the Theia snow products and in situ snow depth data (SD).</caption><thead><tr><th rowspan="2"></th><th colspan="2">Theia snow products</th></tr><tr><th>no snow</th><th>snow</th></tr></thead><tbody><tr><th>In situ snow depth</th><td></td><td></td></tr><tr><th>SD = 0</th><td>276</td><td>8</td></tr><tr><th>SD > 0</th><td>76</td><td>1054</td></tr></tbody></table>

### 4.1.4 Visual verification

The evaluation methods presented above undersample the actual resolution of the Theia snow collection, both in space (only 120 points of comparisons in Sect. 4.1.1) and time (only six dates of comparisons in Sect. 4.1.2). Given that we do not have a more extensive validation data set, we used a time series of 64 Theia snow products from 6 July 2015 to 2 July 2017 over tile 31TCH to control the consistency of the snow and cloud masks based on the visual inspection of the false-colour composites. This approach is efficient for evaluating the frequency of gross errors on the tile scale, i.e. large patches of false snow or false no-snow detection (Vidot et al., 2017).

## 4.2 Results

### 4.2.1 Comparison with in situ snow depth measurements

The confusion matrix between the Theia snow products and in situ snow depth measurements is given in Table 5. We find 1414 pairs of data for the study period (i.e. there are 1414 individual snow depth measurements for which a cloud-free retrieval can be found in the Theia Snow collection on the same day). From this confusion matrix, the accuracy (proportion of correct classifications) is 94% and the kappa coefficient is 0.83, which indicates excellent agreement according to Fleiss et al. (2013). The false-positive rate (2.8%) is lower than the false-negative rate (6.7%), which means that the Theia snow product is more likely to underestimate than to overestimate the snow detection at the station locations. The false-negative rate decreases if SD₀ is set to a higher value; however the false-positive rate also increases in such a way that an optimum is reached at SD₀ = 2 cm (Fig. 10).

The comparisons of each individual time series before matching the data by date is shown in Appendix B. These plots illustrate that the high revisit time of Sentinel-2 enables us to capture the seasonal cycle of the snow cover well. Even intermediate melt-out events at lower elevation stations can be identified over the course of the snow season.

### 4.2.2 Comparison with snow maps of higher spatial resolution

Table 6 shows that the MAJA-LIS workflow provides a better detection of the snow cover than the Sen2Cor processor in

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Location of the five SPOT-6/7 products used for the evaluation of the snow detection accuracy and the corresponding Sentinel-2 tiles (see also Table 4).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Sensitivity of the agreement between the in situ and Theia snow product to SD₀ in metres (threshold to convert the measured snow depth to snow presence or absence). The inset shows a close-up of the region of 0–0.2 m.</figcaption></figure>
<table><caption>**Table 6.** Results of the validation of Theia (LIS 1.2.1) and ESA (Sen2Cor 2.5) snow products using SPOT-6/7 reference snow maps.</caption><thead><tr><th colspan="2">Sentinel-2 product</th><th colspan="2">Accuracy</th><th colspan="2">F1</th><th colspan="2">Kappa</th><th colspan="2">FPR</th><th colspan="2">FNR</th></tr><tr><th>Tile</th><th>Date</th><th>Sen2Cor</th><th>LIS</th><th>Sen2Cor</th><th>LIS</th><th>Sen2Cor</th><th>LIS</th><th>Sen2Cor</th><th>LIS</th><th>Sen2Cor</th><th>LIS</th></tr></thead><tbody><tr><td>31TGK</td><td>11 Dec 2016</td><td>0.83</td><td>0.92</td><td>0.73</td><td>0.91</td><td>0.61</td><td>0.84</td><td>0.41</td><td>0.04</td><td>0.05</td><td>0.14</td></tr><tr><td>31TGL</td><td>11 Mar 2017</td><td>0.89</td><td>0.90</td><td>0.92</td><td>0.93</td><td>0.74</td><td>0.76</td><td>0.09</td><td>0.02</td><td>0.06</td><td>0.11</td></tr><tr><td>31TGL</td><td>1 Dec 2016</td><td>0.84</td><td>0.95</td><td>0.78</td><td>0.95</td><td>0.66</td><td>0.90</td><td>0.35</td><td>0.05</td><td>0.01</td><td>0.06</td></tr><tr><td>31TGK</td><td>11 Mar 2017</td><td>0.85</td><td>0.89</td><td>0.89</td><td>0.92</td><td>0.65</td><td>0.72</td><td>0.15</td><td>0.01</td><td>0.07</td><td>0.13</td></tr><tr><td>31TGK</td><td>13 Aug 2016</td><td>0.98</td><td>0.99</td><td>0.35</td><td>0.77</td><td>0.34</td><td>0.77</td><td>0.78</td><td>0.24</td><td>0.22</td><td>0.22</td></tr><tr><td>32TLS</td><td>12 Oct 2016</td><td>0.98</td><td>0.97</td><td>0.63</td><td>0.64</td><td>0.62</td><td>0.62</td><td>0.49</td><td>0.08</td><td>0.17</td><td>0.51</td></tr></tbody></table>

all the studied cases. Although the improvement in the accuracy coefficient is low in some cases (31TGL/2017-03-11, 31TGK/2016-08-13) or even slightly negative in one case (32TLS/2016-10-12), the F1 score and kappa coefficient of the Theia snow products are always greater than or equal to those of Sen2Cor. The differences are significant for four products among the six tested. The lower accuracy in the case of the 32TLS/2016-10-12 product is due to a higher false-negative rate. The false-positive rates in the Theia snow products are generally lower than the false-negative rates, which is consistent with the previous evaluation using in situ data. The spatial analysis of the errors indicates that the false-positives (i.e. overdetection of snow pixels) in the Theia products are mostly located near the snow cover edges (Fig. A1). This suggests that these errors might not be due to a false detection of snow (e.g. confusion with a lake or a cloud surface) but are rather due to the resolution discrepancy with the reference data. By contrast, the spatial comparison of the Sen2Cor products with the reference data shows large patches of false negatives (omission of snow pixels) in the snow-covered area, which are typically due to the omission of snow pixels in shaded areas (Fig. A2).

### 4.2.3 Visual verification

We found 4 products among 64 with significant patches of false-positive pixels (i.e. pixels falsely classified as snow). These false snow areas are exclusively located in cloud areas (Fig. 11). False-negative pixels (i.e. pixels falsely classified as snow-free) can be found by looking at higher-resolution data along the snow cover edges, but we did not find large areas of false-negative pixels.

## 5 Discussion

The evaluation of the Theia Snow collection using both in situ and remote-sensing reference data sets indicates an excellent accuracy of the snow detection, albeit with a tendency to snow underdetection. The omission of snow pixels can be due to the presence of trees or shadows. In the case of the in situ comparison, part of the errors may be explained by the uncertainty in the geolocalization of the Sentinel-2 data,

which is close to 11 m at 95% for both satellites according to the latest Sentinel-2 L1C data quality report (MPC Team, 2018). In addition, in situ measurements may not be representative of the snow conditions within the Sentinel-2 pixel. However, the sensitivity analysis of the accuracy and kappa coefficient to the SD₀ value suggests that the discrepancy between the scale of the in situ observations and the scale of the pixel observation is not significant, since the optimal SD₀ is close to 0 m. In contrast, a similar analysis with lower-resolution snow products from MODIS found an optimal SD₀ of 0.15 m Gascoin et al. (2015). This illustrates how high-resolution snow products from Sentinel-2 enable us to partly resolve the typical scale issue between in situ and remote-sensing products (Blöschl, 1999).

The visual inspection of the products reveals that the main issue in some Theia snow products is rather the occurrence of false positives due to the confusion with cloud surfaces. This problem can be due to two main factors:

- Cold clouds: cold clouds have a spectral signature that is close to the snow cover since they contain ice crystals (high NDSI). If these clouds are not accurately detected by the level 2A processor then the LIS algorithm will also classify them as snow after pass 1 (Sect. 2). As a result the snow line elevation z_s is not well estimated, which can generate erroneous snow patches within these clouds, even at low elevation (e.g. product of 20 July 2017 in Fig. 11).

- Shaded clouds: the three-dimensional structure of the clouds can form shadows within the cloud cover area (e.g. product of 4 March 2017 in Fig. 11). These shaded cloud areas have a lower reflectance in the visible wavelengths and therefore can be considered as dark clouds in the LIS algorithm, provided that the shaded cloud cover area is large enough to significantly reduce the reflectance in the red band, even after the downsampling step (Sect. 2).

Among 64 products, we find 4 cases with significant false snow cover detection in this tile (Fig. 11). Although it is not shown here, we have made this exercise over other tiles in the Alps and Atlas mountains, and we estimate that this proportion is similar in these areas. Regarding this issue we mainly
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Full time series of Theia snow products from 6 July 2015 to 2 July 2017 over tile 31TCH. Four products with visually evident false snow detections are framed in red and shown at higher resolution. Red arrows indicate the location of the false snow patches, while the green arrow indicates a true snow patch. The product of 16 April 2017 is shown as a reference to help the visualization of the errors.</figcaption></figure>

rely on our visual judgment using colour composites to assess the ability of the algorithm to discriminate the cloud and the snow cover. However, to further quantify the algorithm performance and eventually optimize its parameters, we would need a database of classified imagery with labelled regions of cloud and snow cover in many different cases.

Both Sen2Cor and LIS use the NDSI to detect the snow cover; however the algorithms differ in many aspects. The better performances of LIS in Sect. 4.2.2 can be related to the mono-date approach for atmospheric correction and cloud

detection in Sen2Cor, which can cause snow-cloud confusion in the L2A product, while LIS uses L2A products from the more accurate multi-temporal MAJA processor (Baetens et al., 2019). In addition, a unique feature of the LIS algorithm is the use of the snow line elevation concept to improve the robustness of the snow detection in mountain regions. In addition, LIS was designed to retrieve snow below thin clouds and to reclassify snow pixels which are frequently found near the snow cover edges in L2A products. Apart
from the MAJA processing, however, the effect of these al-
gorithmic differences could be less evident in flat areas.

6 Code and data availability

The Theia snow collection can be accessed
and cited using this digital object identifier:
https://doi.org/10.24400/329360/F7Q52MNK (Gascoin
et al., 2018). The let-it-snow source code is available
under the GNU Affero General Public License v3.0 in this
repository: https://gitlab.orfeo-toolbox.org/remote_modules/
let-it-snow/ (last access: 12 April 2019).

7 Conclusions

The Theia Snow collection is a free collection of snow products which indicate the presence or absence of snow on the land based on Sentinel-2 (20 m) and Landsat-8 (30 m) data. Most of the snow products are derived from Sentinel-2. However, in late September 2018, Landsat-8 snow products started being routinely added to the Theia Snow collection for the titles in France only. Previous Landsat-8 data were also reprocessed back to March 2017. At the time of writing (15 March 2019), for example in tile T31TCH in the Pyrenees, 20.5 % of all available products were derived from Landsat-8 data (63 among 244). In addition, the processing is done to facilitate the combination of Landsat-8 and Sentinel-2 snow products, since the data format and the tiling scheme are the same (only the pixel size differs). Although the co-registration of Sentinel-2 and Landsat-8 data can be still problematic (Storey et al., 2016), the combination of “analysis ready” Landsat-8 and Sentinel-2 is important for increasing the number of cloud-free observations.

The evaluation presented in this article indicate that the
Theia Snow product allows an accurate detection of the snow
cover (and of the absence of snow). However, it remains a
preliminary assessment that should be extended in the fu-
ture. For example, the evaluation was based only on Sentinel-
2 products, since the integration of Landsat-8 data started
during the writing of this paper and should be extended to
Landsat-8 products. However, the processing is done with
the same algorithm, which is based on the calculation of
the NDSI from slope-corrected surface reflectance (level 2A
data). In addition, the evaluation was limited to the French
Alps and French Pyrenees, while the collection covers other
mountain regions. Last, the spatial evaluation using higher-
resolution remote-sensing data was focused on the accuracy
of the snow detection using very high-resolution clear-sky
images, while, as discussed above, the main difficulty is not
to detect the snow when the sky is cloud-free but rather to
avoid false snow detection within clouds. The reduction of
the false-positive rate is the main challenge for the future de-
velopments in the LIS algorithm. In the meantime, we also
welcome feedback from other users.

The Theia Snow collection is based on optical observa-
tions; therefore it is not adapted to the detection of the snow
cover in dense forest areas where the ground is obstructed
by the canopy (Xin et al., 2012). This may typically occur in
evergreen conifer forests of alpine regions (e.g. Alps, Pyre-
nees). We have noted that the snow detection can be success-
ful even in alpine forests, but we lack data that can be used for
quantifying its accuracy. Therefore it is recommended to use
a land cover map to exclude these forest areas from the anal-
ysis. For the tiles in France, this can be done using the Theia
land cover map (Inglada et al., 2017), which is distributed at
the same spatial resolution. In the Pyrenees, the mean alti-
tude of zero-degree isotherm in winter is close to the treeline
elevation at 1600 m; therefore the impact of the forest cover
is limited (Gascoin et al., 2015). The LIS algorithm may also
fail to detect the snow cover on steep, shaded slopes if the
solar elevation is very low (typically below 20°). This can
occur in midlatitude areas in winter. In this case, the slope
correction in the L2A product is generally not applied. This
is indicated in the L2A mask but not in the snow product. A
future release of the Theia Snow collection should include
this information. However, the visual inspection of many im-
ages in regions of complex topography gives us the impres-
sion that the snow detection algorithm already performs well
on shaded slopes. We are planning to investigate this issue
further using alternative approaches (Sirguey et al., 2009).

In spite of all these limitations, the Theia snow products
have already been successfully used for the evaluation of the
MODIS snow products (Masson et al., 2018) and the assim-
ilation of snow-covered area data into a snowpack model in
the High Atlas (Baba et al., 2018). Given that the snow cover
is a key driver of many natural processes in mountains re-
gions, we envision various potential applications of the Theia
snow products, including the modelling of the distribution of
the permafrost in mountain regions, the validation and cal-
ibration of hydrological models in snow-dominated catch-
ments and the spatial modelling of biodiversity and produc-
tivity of ecosystems in mountain regions.
Appendix A: Spatial comparison between Theia and Sen2Cor snow products with SPOT-6/7 reference data

<figure><img src="image_14.png" /><figcaption>**Figure A1.** Theia snow products vs. SPOT-6/7 snow maps (green: true negative, white: true positive, red: false positive, yellow: false negative).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure A2.** Sen2cor snow product vs. SPOT-6/7 snow maps (green: true negative, white: true positive, red: false positive, yellow: false negative).</figcaption></figure>
Appendix B: Time series of Theia snow products and snow depth records at 120 snow observation stations from 1 September 2017 to 31 August 2018

<figure><img src="image_16.png" /><figcaption>Figure B1. Magenta dots: snow presence (a) or absence (b) in the Theia snow collection. Black line: snow depth time series.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure B2. Magenta dots: snow presence (a) or absence (b) in the Theia snow collection. Black line: snow depth time series.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure B3.** Magenta dots: snow presence (a) or absence (b) in the Theia snow collection. Black line: snow depth time series.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure B4.** Magenta dots: snow presence (a) or absence (b) in the Theia snow collection. Black line: snow depth time series.</figcaption></figure>
**Author contributions.** SG and OH defined the LIS algorithm. MG and GS worked on the implementation and development of the LIS processor. MG, MB and SG worked on the validation of the snow products. SG wrote the manuscript with inputs from all the co-authors.

**Competing interests.** The authors declare that they have no conflict of interest.

**Acknowledgements.** The development and production of the Theia Snow collection is supported by the CNES. We thank Tristan Klempka for his contribution to earlier developments of the LIS processor. We thank Marc Leroy, Arnaud Sellé and Philippe Pacholczyk for the coordination of the Theia project and Joelle Donadieu and Céline l'Helguen for the transfer to the production at MUSCATE. This manuscript was greatly improved thanks to the constructive comments of Gaia Piazzi, Javier Herrero, Jeff Dozier and an anonymous referee.

**Review statement.** This paper was edited by Alexander Kokhanovsky and reviewed by Gaia Piazzi, Javier Herrero, Jeff Dozier, and one anonymous referee.

## References

Baba, M. W., Gascoin, S., and Hanich, L.: Assimilation of Sentinel-2 Data into a Snowpack Model in the High Atlas of Morocco, Remote Sensing, 10, 1982, https://doi.org/10.3390/rs10121982, 2018.

Baetens, L., Desjardins, C., and Hagolle, O.: Validation of Copernicus Sentinel-2 Cloud Masks Obtained from MAJA, Sen2Cor, and FMask Processors Using Reference Cloud Masks Generated with a Supervised Active Learning Procedure, Remote Sensing, 11, 433, https://doi.org/10.3390/rs11040433, 2019.

Blöschl, G.: Scaling issues in snow hydrology, Hydrol. Process., 13, 2149–2175, 1999.

Bouchet, M.: Validation et amélioration des produits Surfaces Enneigées Sentinel-2, Zenodo, https://doi.org/10.5281/zenodo.1446460, 2018.

Cohen, J.: A coefficient of agreement for nominal scales, Educ. Psychol. Meas., 20, 37–46, 1960.

Commissariat général au développement durable: Plan d'applications satellitaires 2018, available at: https://www.ecologique-solidaire.gouv.fr/sites/default/files/Plan%20d%E2%80%99applications%20satellitaires%202018.pdf (last access: 12 April 2019), 2018.

Dozier, J.: Spectral signature of alpine snow cover from the Landsat Thematic Mapper, Remote Sens. Environ., 28, 9–22, 1989.

Drusch, M., Bello, U. D., Carlier, S., Colin, O., Fernandez, V., Gascon, F., Hoersch, B., Isola, C., Laberinti, P., Martimort, P., Meygret, A., Spoto, F., Sy, O., Marchese, F., and Bargellini, P.: Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services, Remote Sens. Environ., 120, 25–36, https://doi.org/10.1016/j.rse.2011.11.026, 2012.

Fleiss, J. L., Levin, B., and Paik, M. C.: Statistical methods for rates and proportions, John Wiley & Sons, 2013.

Frei, A., Tedesco, M., Lee, S., Foster, J., Hall, D., Kelly, R., and Robinson, D. A.: A review of global satellite-derived snow products, Adv. Space Res., 50, 1007–1029, https://doi.org/10.1016/j.asr.2011.12.021, 2012.

Gao, B.-C., Goetz, A. F., and Wiscombe, W. J.: Cirrus cloud detection from airborne imaging spectrometer data using the 1.38 µm water vapor band, Geophys. Res. Lett., 20, 301–304, 1993.

Gascoin, S., Hagolle, O., Huc, M., Jarlan, L., Dejoux, J.-F., Szczypta, C., Marti, R., and Sánchez, R.: A snow cover climatology for the Pyrenees from MODIS snow products, Hydrol. Earth Syst. Sci., 19, 2337–2351, https://doi.org/10.5194/hess-19-2337-2015, 2015.

Gascoin, S., Grizonnet, M., Klempka, T., and Salgues, G.: Algorithm theoretical basis documentation for an operational snow cover product from Sentinel-2 and Landsat-8 data (Let-it-snow), https://doi.org/10.5281/zenodo.1414452, 2018.

Gascoin, S., Grizonnet, M., Hagolle, O., and Salgues, G.: Theia Snow collection, CNES for THEIA Land data center, https://doi.org/10.24400/329360/f7q52mnk, 2018.

GDAL/OGR contributors: GDAL/OGR Geospatial Data Abstraction software Library, Open Source Geospatial Foundation, http://gdal.org (last access: 12 April 2019), 2018.

Grizonnet, M., Michel, J., Poughon, V., Inglada, J., Savinaud, M., and Cresson, R.: Orfeo ToolBox: open source processing of remote sensing images, Open Geospatial Data, Software and Standards, 2, 15, https://doi.org/10.1186/s40965-017-0031-6, 2017.

Hagolle, O., Huc, M., Pascual, D. V., and Dedieu, G.: A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENUS, LANDSAT and SENTINEL-2 images, Remote Sens. Environ., 114, 1747–1755, 2010.

Hagolle, O., Huc, M., Desjardins, C., Auer, S., and Richter, R.: MAJA Algorithm Theoretical Basis Document, https://doi.org/10.5281/zenodo.1209633, 2017.

Hall, D., Riggs, G., Salomonson, V., DiGirolamo, N., and Bayr, K.: MODIS snow-cover products, Remote Sens. Environ., 83, 181–194, 2002.

Inglada, J., Vincent, A., Arias, M., Tardy, B., Morin, D., and Rodes, I.: Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series, Remote Sensing, 9, 95, https://doi.org/10.3390/rs9010095, 2017.

Jarvis, A., Guevara, E., Reuter, H. I., and Nelson, A. D.: Hole-filled SRTM for the globe: version 4: data grid, CGIAR Consortium for Spatial Information, Washington D.C., USA, 2008.

Klein, A. and Barnett, A.: Validation of daily MODIS snow cover maps of the Upper Rio Grande River Basin for the 2000-2001 snow year, Remote Sens. Environ., 86, 162–176, https://doi.org/10.1016/S0034-4257(03)00097-X, 2003.

Krajčí, P., Holko, L., and Parajka, J.: Variability of snow line elevation, snow cover area and depletion in the main Slovak basins in winters 2001–2014, J. Hydrol. Hydromech., 64, 12–22, 2016.

Li, J. and Roy, D. P.: A Global Analysis of Sentinel-2A, Sentinel-2B and Landsat-8 Data Revisit Intervals and Implications for Terrestrial Monitoring, Remote Sensing, 9, 902, https://doi.org/10.3390/rs9090902, 2017.

Louis, J., Debaecker, V., Pflug, B., Main-Korn, M., Bieniarz, J., Mueller-Wilm, U., Cadau, E., and Gascon, F.: Sentinel-2 Sen2Cor: L2A Processor for Users, in: Proc. 'Living Planet Sym-
posium 2016', Prague, Czech Republic, 9–13 May 2016, ESA SP-740, 9–13, 2016.

Malnes, E., Buanes, A., Nagler, T., Bippus, G., Gustafsson, D., Schiller, C., Metsämäki, S., Pulliainen, J., Luojus, K., Larsen, H. E., Solberg, R., Diamandi, A., and Wiesmann, A.: User requirements for the snow and land ice services – CryoLand, The Cryosphere, 9, 1191–1202, https://doi.org/10.5194/tc-9-1191-2015, 2015.

Masson, T., Dumont, M., Mura, M., Sirguey, P., Gascoin, S., Dedieu, J.-P., and Chanussot, J.: An Assessment of Existing Methodologies to Retrieve Snow Cover Fraction from MODIS Data, Remote Sensing, 10, 619, https://doi.org/10.3390/rs10040619, 2018.

Metsämäki, S., Pulliainen, J., Salminen, M., Luojus, K., Wiesmann, A., Solberg, R., Böttcher, K., Hiltunen, M., and Ripper, E.: Introduction to GlobSnow Snow Extent products with considerations for accuracy assessment, Remote Sens. Environ., 156, 96–108, 2015.

MPC Team: Sentinel-2 L1C Data Quality Report, European Space Agency, Tech. Rep. 31, 80 pp., 2018.

Painter, T. H., Rittger, K., McKenzie, C., Slaughter, P., Davis, R. E., and Dozier, J.: Retrieval of subpixel snow covered area, grain size, and albedo from MODIS, Remote Sens. Environ., 113, 868–879, 2009.

Parajka, J. and Blöschl, G.: Spatio-temporal combination of MODIS images–potential for snow cover mapping, Water Resour. Res., 44, W03406, https://doi.org/10.1029/2007WR006204, 2008.

Ramsay, B. H.: The interactive multisensor snow and ice mapping system, Hydrol. Process., 12, 1537–1546, 1998.

Sirguey, P., Mathieu, R., Arnaud, Y., Khan, M., and Chanussot, J.: Improving MODIS Spatial Resolution for Snow Mapping Using Wavelet Fusion and ARSIS Concept, IEEE Geosci. Remote S., 5, 78–82, https://doi.org/10.1109/LGRS.2007.908884, 2008.

Sirguey, P., Mathieu, R., and Arnaud, Y.: Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment, Remote Sens. Environ., 113, 160–181, 2009.

Storey, J., Roy, D. P., Masek, J., Gascon, F., Dwyer, J., and Choate, M.: A note on the temporary misregistration of Landsat-8 Operational Land Imager (OLI) and Sentinel-2 Multi Spectral Instrument (MSI) imagery, Remote Sens. Environ., 186, 121–122, https://doi.org/10.1016/j.rse.2016.08.025, 2016.

U.S. Geological Survey Earth Resources Observation And Science Center: Landsat OLI Level-2 Surface Reflectance (SR) Science Product, https://doi.org/10.5066/f78s4mzj, 2014.

Vidot, J., Bellec, B., Dumont, M., and Brunel, P.: A daytime VIIRS RGB pseudo composite for snow detection, Remote Sens. Environ., 196, 134–139, https://doi.org/10.1016/j.rse.2017.04.028, 2017.

Xin, Q., Woodcock, C. E., Liu, J., Tan, B., Melloh, R. A., and Davis, R. E.: View angle effects on MODIS snow mapping in forests, Remote Sensing of Environment, 118, 50–59, https://doi.org/10.1016/j.rse.2011.10.029, 2012.

Zhu, Z., Wang, S., and Woodcock, C. E.: Improvement and expansion of the Fmask algorithm: Cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images, Remote Sens. Environ., 159, 269–277, 2015.