Open AccessArticle

An Adaptive Method for the Estimation of Snow-Covered Fraction with Error Propagation for Applications from Local to Global Scales

Lars Keuris

Markus Hetzenecker

Thomas Nagler

Nico Mölg

and

Gabriele Schwaizer

ENVEO IT GmbH, Fürstenweg 176, 6020 Innsbruck, Austria

Author to whom correspondence should be addressed.

Remote Sens. 2023, 15(5), 1231; https://doi.org/10.3390/rs15051231

Submission received: 10 January 2023 / Revised: 15 February 2023 / Accepted: 20 February 2023 / Published: 23 February 2023

(This article belongs to the Special Issue Remote Sensing of the Cryosphere)

Download

Browse Figures

Versions Notes

Abstract

Snow can cover over 50% of the landmass in the Northern Hemisphere and has been labelled as an Essential Climate Variable by the World Meteorological Organisation. Currently, continental and global snow cover extent is primarily monitored by optical satellite sensors. There are, however, no large-scale demonstrations for methods that (1) use all the spectral information that is measured by the satellite sensor, (2) estimate fractional snow and (3) provide a pixel-wise quantitative uncertainty estimate. This paper proposes a locally adaptive method for estimating the snow-covered fraction (SCF) per pixel from all the spectral reflective bands available at spaceborne sensors. In addition, a comprehensive procedure for root-mean-square error (RMSE) estimation through error propagation is given. The method adapts the SCF estimates for shaded areas from variable solar illumination conditions and accounts for different snow-free and snow-covered surfaces. To test and evaluate the algorithm, SCF maps were generated from Sentinel-2 MSI and Landsat 8 OLI data covering various mountain regions around the world. Subsequently, the SCF maps were validated with coincidentally acquired very-high-resolution satellite data from WorldView-2/3. This validation revealed a bias of 0.2% and an RMSE of 14.3%. The proposed method was additionally tested with Sentinel-3 SLSTR/OLCI, Suomi NPP VIIRS and Terra MODIS data. The SCF estimations from these satellite data are consistent (bias less than 2.2% SCF) despite their different spatial resolutions.

Keywords:

snow cover; multi-spectral; unmixing; Sentinel; Landsat; error estimation; mountain

1. Introduction

1.1. Background

Snow is a pivotal element within the ensemble of cryospheric climate components and impacts the Earth at various scales. On average, snow cover reaches up to approximately 44% (47.4 ± 1.5 million km

^{2}

) of the land area in the Northern Hemisphere during the winter period [1]. The vast extent and its high albedo cause a drop in the Earth’s net energy balance as snow reflects a large portion of the incoming solar radiation. The World Meteorological Organisation (WMO) has listed snow-covered area as an Essential Climate Variable (ECV) for its indisputable role in the Earth’s climate system. More than one billion people depend on seasonal snow and glaciers for their water supply [2]. Moreover, snow melt runoff is an important driver for the generation of electricity through hydropower stations [3] and is especially relevant during the current energy transition towards greener sources. At the same time, snow melt runoff can take place rapidly—potentially flooding populated areas [4,5]. In essence, the location, state and evolution of snow is a highly relevant topic of today for many world citizens, scientists and policymakers alike.

Due to the enormous area that can be covered by seasonal snow each year, satellite remote sensing is currently the most feasible tool for mapping the world’s spatially largest cryospheric component. Currently, there are numerous continental-to-global-scale snow-cover products available (e.g., [6,7,8,9,10,11,12,13,14]). The majority of these are binary snow products with a medium-to-low ground resolution (0.005 to 0.25

^{°}

corresponding to about 0.5 to 25 km). Most of these products use only a few spectral bands to discriminate snow-covered from snow-free areas. Moreover, a quantitative uncertainty estimate for individual pixels does not exist for binary snow products.

1.2. Previous Work

Three popular groups for the estimation of snow-covered fraction (SCF) have been identified and are discussed hereafter. Initially, SCF mapping efforts have widely used the normalized difference snow index (NDSI) [15]. The NDSI takes advantage of the high reflectance of snow in the visible part of the spectrum (VIS) and the low reflectance of snow in the shortwave infrared part of the spectrum (SWIR). Accordingly, the NDSI is computed from a VIS and SWIR top-of-atmosphere (TOA) reflectance band,

r_{VIS}

and

r_{SWIR}

, respectively, in Equation (1):

NDSI = \frac{r_{VIS} - r_{SWIR}}{r_{VIS} + r_{SWIR}}

(1)

The NDSI ranges from −1 to 1 and positively correlates with the SCF. A popular approach is to scale the NDSI linearly (e.g., [16,17]) or non-linearly (e.g., [18]) to directly estimate the SCF. An advantage of this approach is that it provides a simple and robust solution with only two reflective bands that are available for most optical satellite sensors. Furthermore, a normalized difference metric reduces atmospheric transfer effects.

Correspondingly, NDSI maps from the MODerate resolution Imaging Spectroradiometer (MODIS) onboard Terra (MOD10_L2) and Aqua (MYD10_L2) are computed and disseminated on a daily basis [19]. Together with the long time series of data acquired by this sensor (>20 years), this has led to a widespread use of NDSI for both small and large scale snow mapping applications (e.g., [17,20,21]). Furthermore, the Pan-European High-Resolution Snow and Ice (HR-I&S) monitoring system as part of the Copernicus Land Monitoring Service (CLMS), which is coordinated by the European Environmental Agency (EEA), applies a non-linear scaling of the NDSI from Sentinel-2 MSI data to estimate the SCF [22].

An approach that can fully utilize the spectral information available is multispectral unmixing. Multispectral unmixing separates an observed spectrum into fractions of known reference spectra (called endmembers). Endmembers can be extracted from a pre-defined spectral library or, alternatively, a model (or a combination with a model) is introduced to further improve the representation of the endmembers [23,24,25]. Even though the reported validations vary, multispectral unmixing has demonstrated an overall improvement in SCF estimates compared to NDSI-derived SCF estimates [26,27,28,29,30].

A third approach uses machine learning, which commonly models the input bands as well as auxiliary information (e.g., normalized difference indices, topography and land cover type) to SCF in a linear or non-linear fashion (e.g., [31,32,33,34]). Although machine-learning methods are more complex than NDSI-based methods, they can incorporate auxiliary information more easily. Ref. [34] successfully demonstrated that a machine-learning method, such as Support Vector Machine (SVM), can provide reasonable SCF estimates in complex terrain (i.e., in alpine terrain with a low sun elevation and fractional snow coverage).

1.3. Limitations of Existing Methods

Although NDSI-based methods are popular due to their simplicity, they have a significant weakness in shaded regions. In fact, shade amplifies NDSI estimates without any snow present, which commonly leads to overestimated SCFs in the shade for NDSI-based algorithms. Shade can be cast by clouds and topography and is prominent in high-resolution optical remote-sensing imagery. In alpine terrain, north-facing slopes (in the Northern Hemisphere) are commonly poorly illuminated during the winter. In these shaded areas, the reduced observed radiance is dominated by diffuse light instead of (direct) specular light, such that atmospheric scatterers play a relatively large role [35,36].

In the Earth’s atmosphere, the scattering of solar light is especially dominant for shorter wavelengths (within the visible and IR part of the spectrum; Rayleigh scattering), such that the relative difference in reflectance between shorter and longer wavelengths increases in shaded regions. Hence, the NDSI is amplified, and the SCF is potentially overestimated. It must also be noted that scaling down well-illuminated snow spectra by a single constant to model their shaded counterpart is, from this perspective, incorrect and can lead to inaccurate SCF estimates. An example of the behaviour of the NDSI in illuminated and shaded areas is shown in Figure 1.

Figure 1 shows that the NDSI is nearly as sensitive to shade as it is to snow. The methods by [17], i.e., FRA6T and FRA7U, linearly scale the NDSI to estimate the SCF and were originally developed for MODIS onboard Aqua and Terra, respectively. The required bands for calculating the NDSI onboard MODIS Aqua and Terra have a native spatial resolution of approximately 500 m. At such resolutions, shade plays a smaller role than at high resolutions. Figure 1 shows that the NDSI is unsuitable for the detection of snow-free and snow-covered areas in shaded terrain, which is particularly prominent in mountainous areas.

Alternatively, multispectral unmixing and machine-learning methods are more versatile in terms of data input and have demonstrated promising results [3,24,25,32,33,37]. Current prominent multispectral unmixing implementations use endmember spectra from a spectral library, which hinders a straightforward implementation for global-scale applications. One of the main reasons for this is that snow-free and snow-covered signatures vary widely both locally and globally. This is demonstrated in Figure 2 where the spectral signatures of three different snow-free and fully snow-covered pixels are displayed from (1) well-illuminated pixels on the left-hand side and (2) shaded pixels on the right-hand side. These were extracted from the marked pixel locations in Figure 1, where the round and diamond markers indicate the well-illuminated and shaded pixels, respectively. Evidently, more challenging spectral signatures can be expected for continental and global-scale applications.

1.4. Objective

There have been promising methodological developments for the estimation of the snow-covered fraction (SCF in percent per pixel) from remote-sensing data, among which linear multispectral unmixing approaches (e.g., [3,23,24,25,30,36,38]) and machine-learning methods (e.g., [31,32,33,39]) are the most prominent. These advanced methods have been tested and applied primarily to high-resolution satellite data from Landsat 8 (OLI) and Sentinel-2 (MSI) and with regional coverage to Terra (MODIS) data (see Table 1 for the sensor abbreviations).

In this paper, we propose a new method for estimating the SCF from multispectral optical satellite data from high-to-medium spatial resolution (10 m to 1 km) sensors, i.e., those onboard Sentinel-2 (MSI), Landsat 8/9 (OLI(-2)), Sentinel-3 (OLCI and SLSTR), Suomi National Polar-orbiting Partnership (NPP) (VIIRS) and Terra (MODIS). As the proposed method adapts to local variations of the solar illumination and spatial variations of snow-free and snow-covered land surfaces, it is applicable to a wide range of environments. Additionally, a pixel-wise estimate of the SCF uncertainty, based on error propagation, is included.

2. Data and Methodology

2.1. Input Data

The proposed method uses the TOA reflective bands as input to estimate the visible SCF. The reflective bands from a variety of sensors are used independently to demonstrate the multi-resolution applicability of the method. The high-resolution SCF maps are calculated from data from Sentinel-2A/B MSI and Landsat 8/9 OLI(-2). Four more medium-resolution sensors are used in this paper for demonstration: OLCI and SLSTR onboard Sentinel-3A/B, VIIRS onboard the Suomi NPP and MODIS onboard Terra.

The reflective bands used in this paper are shown in Figure 3, and their band names can be found in Table 1. The reflective bands are aggregated or resampled to create a raster pyramid (25 m (MSI, OLI), 250 m (OLCI, SLSTR), 500 m (VIIRS) and 1000 m (MODIS)), which is used for intercomparison over the European Alps.

The European Space Agency (ESA) Climate Change Initiative (CCI) land-cover map from 2015 [40] and the Copernicus Global Surface Water dataset (maximum extent) [41] were used to mask water bodies in medium- and high-resolution maps, respectively. Five very-high-resolution (2 m × 2 m and higher) WorldView-2/3 reference scenes were available for validation. The WorldView-2/3 data covers different climate zones, different times of the year and different terrains (see Table A2 for more details), such that the impact on the final validation results of these external factors is minimized. These scenes were manually checked for minimal snow fall or melt that could have taken place between the WorldView-2/3 acquisition (3 days or less) and the Sentinel-2 or Landsat 8 acquisition that is used for comparison. The Arctic DEM (Digital Elevation Model) ([42], version 7) and the EU-DEM ([43], version 1.1) were used for orthorectification of the WorldView-2/3 images.

2.2. Methodology

The proposed viewable SCF estimation methodology was split up into four steps: (1) pre-processing, (2) local endmember selection, (3) Locally Adaptive MultiSpectral Unmixing (LAMSU) and (4) post-processing. A detailed outline of the workflow is given in Figure 4.

2.2.1. Pre-Processing

In the first step, the data are converted from digital numbers (DN) to TOA reflectance values. In the next step, all pixels that are not used for further evaluation in the processing line are masked out. For instance, pixels affected by clouds in alpine areas are masked out. The Simple Cloud Detection Algorithm (SCDA v2.3) (ESA Snow Climate Change Initiative, Algorithm Theoretical Basis Document, in review), which is based on a hierarchical decision tree, is used for cloud detection if applicable. Furthermore, open-water areas are masked out with the Global ESA CCI land-cover map ([40], v2.0.7) or the Copernicus Global Surface Water dataset [41].

2.2.2. Local Endmember Selection

The selection of endmembers for snow-free and snow-covered pixels is a crucial step for estimating the SCF with multispectral unmixing. The large number of endmembers required for an adequate representation of the surface classes around the world is a major drawback for the application of multispectral unmixing with a spectral library. Instead, a local and adaptive solution is proposed (advancing the work by [38,44]) to accommodate the highly variable nature of surface class endmembers on a global scale. Our approach is to create two groups of endmembers, i.e., snow-free and snow-covered, and to directly select these from the pixels within the image. Snow-free endmembers capture various surface types (grass, sand, rock, etc.) that are typical for a scene and vary for different environments and regions.

Snow-covered endmembers capture various snow spectra, e.g., with varying grain sizes. The SCF of the selected snow-free and snow-covered pixels is directly estimated at 0 and 100% SCF, respectively. Additionally, different illumination conditions are considered, focusing on two classes: shaded areas (only diffuse solar light) and illuminated areas. Shaded areas are particularly prominent in mountainous regions and often coincide with (fractional) snow cover. The SCFs of the pixels that are not selected are estimated using linear spectral unmixing with the identified snow-free and snow-covered endmembers from within the scene.

In short, each pixel potentially belongs to one of four endmember classes: (1) illuminated and snow-free, (2) illuminated and fully snow-covered, (3) shaded and snow-free or (4) shaded and fully snow-covered as illustrated in Figure 5.

The identification of the four endmember classes takes place with one hierarchical decision tree for shaded and illuminated pixels and two hierarchical decision trees for snow-free and snow-covered pixels (see Figure 5). The hierarchical decision trees have been manually derived from 19 independent Landsat-8 scenes from different climate zones, during different seasons and involved repeated and careful visual inspection.

The different rules in the hierarchical decision tree use five spectral bands at wavelengths of 560, 650, 860, 1610 and 2200 nm, which are commonly available from optical sensors onboard Sentinel-2 and Landsat 8/9 as well as medium-resolution sensors, e.g., SLSTR, VIIRS and MODIS. Hence, the methodological set-up allows for a broad (multi-sensor) applicability of the rules for identifying endmembers. The NDSI [15], normalized difference vegetation index (NDVI) [45], two-band enhanced vegetation index (EVI) [46] and eight more normalized spectral indices were calculated from those five spectral bands. The general characteristics for each of the spectral indices are described in Table 2. Snow is a strong reflector in the visible and is nearly opaque in the shortwave infrared part of the spectrum (see Figure 1). Shade decreases the TOA reflectance of the observed spectrum—in particular, for longer wavelengths.

Multiple spectral indices with similar characteristics but with different band combinations are used to increase the robustness of the selection procedure. Moreover, not all snow-free and snow-covered pixels must be found, only those with high certainty. Hence, conservative thresholds are applied to various combinations of these spectral indices. Additionally, the shaded pixels are compared to nearby illuminated pixels, which are also considered shaded if their vector norm is of a similar magnitude (up to 125%).

Subsequently, the 5th through the 95th percentiles (in steps of 5 percent) of the spectral norm are selected from each of the four endmember classes. Finally, the spectral information divergence (SID) [47] is computed between the remaining unclassified pixels and the percentiles of each of the four endmember classes. Unclassified pixels that are similar enough to one of the four endmember classes are then added to the respective endmember class. SID is a probabilistic metric of dissimilarity, where a SID of 0 “nats” indicates equal spectra, and is formulated as in Equation (2):

SID ({\vec{r}}_{a}, {\vec{r}}_{b}) = \sum_{i = 0}^{M} r_{a, i} \cdot log (\frac{r_{a, i}}{r_{b, i}}) + \sum_{i = 0}^{M} r_{b, i} \cdot log (\frac{r_{b, i}}{r_{a, i}})

(2)

where

{\vec{r}}_{a}

and

{\vec{r}}_{b}

are normalized probability vectors of length M (equal to the number of bands) derived from spectral reflectances for pixel a and for pixel b, respectively.

r_{a, i}

and

r_{b, i}

are the ith elements of probability vector

{\vec{r}}_{a}

and

{\vec{r}}_{b}

, respectively. Ref. [47] demonstrated that SID has an average 10-times better relative spectral discriminability power (indicating class discriminability) with respect to the spectral similarity used by the Spectral Angle Mapper (SAM).

A pixel is added to the corresponding endmember class if the SID is below 0.0006, and the NDSI is (1) above 0.75 (for illuminated snow), (2) above 0.85 (for shaded snow), (3) below 0.15 (for illuminated snow-free) and (4) below 0.90 (for shaded snow-free). Shaded snow-covered endmembers are discarded as endmembers if they are near water pixels (three pixels in distance) to avoid false positives. Finally, snow-free and snow-covered endmembers that are adjacent to each other are discarded as endmembers. The newly identified snow-free and snow-covered pixels are assigned a 0 and 100% SCF, respectively.

A separate Voronoi diagram is created from the endmembers of each class, which allows for rapid nearest neighbour querying in a topological space. Endmembers are only included for the creation of the Voronoi diagrams if they are not directly adjacent to four endmembers of the same endmember class. This further speeds up the Voronoi diagram creation and querying. The Voronoi diagrams are created with Fortune’s algorithm [48] using Mathias Westerdahl’s implementation (https://github.com/JCash/voronoi, accessed on 30 June 2022).

2.2.3. Locally Adaptive MultiSpectral Unmixing

The SCF of each remaining unclassified pixel is derived from a set of local endmembers with linear multi-spectral unmixing. The Locally Adaptive MultiSpectral Unmixing (LAMSU) method assumes that an observed spectrum can be made up out of two endmembers: snow-free and snow-covered. This concept can be written as in Equation (3) and describes the forward problem.

{\vec{R}}_{o b s} = \sum_{i = 1}^{N} {\vec{R}}_{i} f_{i} + \vec{ϵ}

(3)

where

{\vec{R}}_{o b s}

is the reflectance column vector with reflectances of the observed spectrum for different wavelengths on each row,

f_{i}

is the fraction of endmember column vector

{\vec{R}}_{i}

that is present in

{\vec{R}}_{o b s}

\vec{ϵ}

is the residual error vector and N is equal to the number of endmembers (2 in this case). Hence, only one (shaded or illuminated) snow-free (0 % SCF) and one (shaded or illuminated) snow (100 % SCF) endmember are used. The forward model has the form

\vec{y} = A \vec{x} + \vec{e}

, such that the least-squares solution can be achieved with Equation (4) under the assumption that the noise in

\vec{y}

is uncorrelated (Gauss–Markov theorem).

\hat{\vec{x}} = {(A^{T} A)}^{- 1} A^{T} \vec{y}

(4)

where

\hat{\vec{x}}

is the vector with unknown fractions for the endmembers in design matrix

A

with one snow-free and one snow-covered endmember column and

\vec{y}

is the observation vector that contains the reflectances of the unclassified pixel. A sum-to-one constraint is added to Equation (4) to ensure that the estimated fractions in

\hat{\vec{x}}

add up to one. This can be realised by introducing a row with ones to

\vec{y}

and

A

. The default lower and upper limit for the fractions in

\hat{\vec{x}}

are zero and one, respectively. The lower and upper limit are accomplished with an active/passive set least-squares algorithm [49]. The unbiased model variance, also known as misfit, (

Var (A \hat{\vec{x}})

) can be derived with Equation (5).

Var (A \hat{\vec{x}}) = \frac{{(\vec{y} - A \hat{\vec{x}})}^{T} (\vec{y} - A \hat{\vec{x}})}{M + 1 - N} = \frac{{\vec{r}}^{T} \vec{r}}{M + 1 - N}

(5)

where

\vec{r}

is the residual vector, M is the number of spectral bands and N is the number of estimated parameters (2 in this case). The estimated fractions (and misfits) are no longer unbiased because of the value constraints placed on the outcome of

\hat{\vec{x}}

. This bias can be calculated from the residuals with Equation (6):

bias (A \hat{\vec{x}}) = \frac{1}{M + 1} \sum_{i}^{M + 1} (y_{i} - A_{i} \hat{\vec{x}}) = \frac{1}{M + 1} \sum_{i}^{M + 1} r_{i}

(6)

where i selects the ith row of

\vec{y}

A

and

\vec{r}

. The misfit mean square error (MSE

_{misfit}

) is then formulated as in Equation (7):

{MSE}_{misfit} = Var (A \hat{\vec{x}}) + bias {(A \hat{\vec{x}})}^{2}

(7)

The MSE

_{misfit}

can be propagated to the estimated fractions in

\hat{\vec{x}}

with Equation (8) under the same assumption of uncorrelated observational errors.

Q_{\hat{\vec{x}} \hat{\vec{x}}} = {MSE}_{misfit} \cdot {(A^{T} A)}^{- 1}

(8)

where

Q_{\hat{\vec{x}} \hat{\vec{x}}}

is the variance–covariance matrix of

\hat{\vec{x}}

. The diagonal elements of

Q_{\hat{\vec{x}} \hat{\vec{x}}}

represent the variances of the estimated snow-free and snow-covered fraction, i.e.,

Var (1 - SCF)

and

Var (SCF)

, respectively.

Var (SCF)

does not account for the uncertainty of the design matrix set-up, which is dominated by the uncertainty of selected endmembers. This uncertainty (

{MSE}_{model}

) is approximated in Section 3.1 and added to

Var (SCF)

for the final uncertainty description with Equation (9).

{MSE}_{tot} = Var (SCF) + {MSE}_{model}

(9)

where MSE

_{tot}

is the total MSE for the SCF estimate.

The snow-free and snow-covered endmembers in

A

are taken from the vicinity of the unclassified pixel using the Voronoi diagrams. If the unclassified pixel is located in the shade, then the five closest snow-free and snow-covered endmembers from the shaded group are used. Otherwise, the nearest five are taken from the illuminated snow-free and snow-covered classes. Additionally, another five snow-free and snow-covered endmembers are looked up in the average opposite direction (with respect to the query point) of the first five closest endmembers. This additional step partially accounts for the fact that the closest endmember may not always be the most suitable endmember. The inclusion of ten snow-free and ten snow-covered endmembers improves the robustness against invalid endmembers of the LAMSU method and will be further improved in the next step.

The individual vector norms of the discovered snow-free and snow-covered endmembers are normalized to the 50th percentile vector norm from their corresponding endmember class if the distance to the query point is 50 pixels or more. Endmembers adjacent to the query point are normalized to the average vector norm of the corresponding set of nearest endmembers. The weighting is linearly scaled for distances in between. Subsequently, up to 100 (10 snow-free times 10 snow-covered endmember combinations) SCF solutions are calculated. The MSE for each combination is calculated with Equation (9). The impact of incorrect endmembers is further reduced by discarding the upper quantile of all the MSEs. The remainder is used as weighting for deriving the weighted average SCF (

\bar{SCF}

) with Equation (10) and the weighted average SCF variance (

Var (\bar{SCF})

) with Equation (11).

\bar{SCF} = \sum_{k} \frac{{SCF}_{k}}{{MSE}_{k}} / \sum_{k} \frac{1}{{MSE}_{k}}

(10)

and

Var (\bar{SCF}) = \sum_{k} \frac{Var ({SCF}_{k})}{{MSE}_{k}} / \sum_{k} \frac{1}{{MSE}_{k}}

(11)

are the equations for weighting the single pixel SCF estimates (

{SCF}_{k}

) and SCF variances (

Var ({SCF}_{k})

), respectively, using the individual linear multispectral unmixing uncertainty estimates

{MSE}_{k}

(Equation (7)). The SCF estimate and the square root of the corresponding uncertainty estimate in Equations (10) and (11), respectively, are the final estimates for a single pixel.

2.2.4. Post-Processing

The range of observed reflectances in shaded regions is narrower with respect to the range of reflectances observed in illuminated regions. The signal-to-noise ratio (SNR) is significantly lower in shaded regions and the calculation of accurate SCF estimations becomes challenging. A direct consequence of the low SNR is the false detection of low SCFs in the output. Groups of adjacent snow pixels in shaded regions with low SCF (less than 5% SCF average, less than 12% SCF maximum) are, therefore, set to 0% SCF.

Classification errors in the static water mask can lead to false detections of low SCFs near water bodies. A pixel with a low SCF is, therefore, set to 0% if all SCFs within a seven pixel radius (including the centre) are 5% SCF or lower and a water pixel (from the water mask) is present. Finally, the transition zones from illuminated to shaded regions are difficult to resolve. These transition zones are smoothed using a 5 × 5 1-

σ

scaled Gaussian and majority filter window as given by Equation (12):

\begin{matrix} {SCF}_{s} = \frac{1}{\sum_{x} \sum_{y} f_{x, y}} \sum_{x} \sum_{y} f_{x, y} {SCF}_{x, y} \end{matrix}

(12)

\begin{matrix} f_{x, y} = \{\begin{matrix} 0, & if \sqrt{x^{2} + y^{2}} > 2 \\ \frac{0.25}{\sqrt{2 π σ}} e^{- \frac{x^{2} + y^{2}}{2 σ^{2}}}, & if (x, y) is shaded \\ \frac{1}{\sqrt{2 π σ}} e^{- \frac{x^{2} + y^{2}}{2 σ^{2}}}, & if (x, y) is illuminated \end{matrix} \end{matrix}

(13)

where

{SCF}_{s}

is the smoothed SCF value,

f_{x, y}

is the weight assigned to

{SCF}_{x, y}

, where x and y both range from −2 to 2 and

σ

is set to one and is the standard deviation of the Gaussian window. Equation (13) gives higher weight to illuminated SCF estimates than to shaded SCF estimates. If more than half of the individual

{SCF}_{x, y}

are snow-free or fully snow-covered, then

{SCF}_{s}

is also set to snow-free or fully snow-covered, respectively. The squared SCF difference between the former and updated SCF is added to the MSE of the relevant pixel for each of the adjustments made.

2.3. Validation of LAMSU SCF Estimates

For validation, snow products from very-high-resolution data were used. The approach from [25] was followed, using WorldView-2/3 data for validation of the LAMSU SCF estimates. The WorldView-2/3 images were orthorectified using the enclosed Rational Polynomial Coefficient (RPC) sensor model and a DEM. The orthorectification was manually inspected and unresolved regions were manually masked out. Subsequently, a subset of snow-free and snow-covered pixels was manually selected and extracted as training data. Images were trained independently with the gradient tree boosting classification method [50] from the XGBoost Python module (version 1.5.2, https://xgboost.readthedocs.io, accessed on 2 July 2022).

Gradient tree boosting is a state-of-the-art machine-learning technique for classification tasks and generally outperforms random forest classifications [51]. The required training data were adjusted until convergence, and a snow-free and snow classification map was produced. Water bodies, clouds and the edges of the WorldView scenes were manually delineated and masked out. Finally, the very high-resolution binary snow classification was averaged to the resolution of the intersecting LAMSU SCF map derived from Landsat 8 or Sentinel-2 data such that pixel-wise comparison became possible. Additionally, the SCF was also computed with FRA6T [17], an NDSI-based approach, which serves as a baseline.

The WorldView-2/3 scenes’ unequal distribution of snow-free and snow-covered pixels almost certainly leads to biased statistics if untreated. Hence, random sampling was used to resolve the unequal distribution. The average bias (sample minus reference) and the average root-mean-square error (RMSE) were calculated from one million random realisations. Each single realisation was independent and randomly drawn and had 50% snow-free and 50% snow-covered pixels. The set size was equal to 95% of the number of pixels from the smallest class (snow-free or snow) within the WorldView-2/3 scene and was the same for each iteration. The average bias and average RMSE were also calculated for all scenes combined in the same manner.

Finally, the consistency of the SCF estimates was also tested through intercomparison of medium-resolution LAMSU SCF estimates from Sentinel-3 SLSTR/OLCI, Suomi NPP VIIRS and Terra MODIS. An SCF map was computed over the European Alps on 2 April 2020 for each of the medium-resolution sensors. The medium-resolution maps were then aggregated and matched to assess the consistency of the algorithm at various spatial scales. They were also compared to selected aggregated Sentinel-2 MSI and Landsat 8 OLI SCF maps.

The consistency between each sensor pair was quantified using the bias and RMSE and was called the multi-sensor bias and multi-sensor RMSE. Seven suitable Sentinel-2 MSI and two Landsat 8 OLI scenes were found to overlap the area within one days’ time. The same sampling strategy was applied as for the WorldView-2/3 data to accommodate unequal class sizes. Water bodies and clouds were masked out with the Global ESA CCI Land cover classification map (v2.0.7) and a manual delineation, respectively.

3. Results

In the first part of this section, LAMSU SCF estimates are compared for validation to SCF estimates derived from very-high-resolution Worldview-2/3 data using a supervised machine-learning technique. Subsequently, its performance is illustrated through the application to high and medium-resolution data. Additional attention is paid to shaded regions and consistency.

3.1. Validation with WorldView-2/3 Imagery

The validation results of the Landsat 8 (30 m) and Sentinel-2 (20 m) SCF product is shown in Table 3. WorldView-2/3 reference scenes 2, 3 and 5 are strongly affected by shade and have the highest RMSEs, reaching up to 21% SCF. The increased bias and RMSE for reference scene 2 can be explained by the glaciated regions within the scene. LAMSU estimates high SCFs for ice-covered pixels, whereas the WorldView-2 reference scene displays low SCFs for the same pixels. The individual biases are equal to or smaller than 8.1% SCF and are approximately 0.2% SCF on average.

In all cases, the LAMSU SCF estimates have noticeably smaller biases (7.1% smaller overall) and RMSEs (9.2% smaller overall) than the FRA6T SCF estimates. The NDSI-based method also shows a positive bias on all accounts, which are a consequence of the overestimated SCF in shaded regions.

The bias and RMSE are assumed to mainly originate from the LAMSU method and its input data quality. Although minor, they can also come from (1) different detection limits between the OLI and MSI sensors on one hand and the WorldView-2/3 sensor on the other, (2) orthorectification and geocoding errors and (3) temporal changes (3 days for reference scenes 3 and 5). The upper limit of the MSE

_{tot}

(Equation (9)) can be approximated from the RMSE results in Table 3.

Two assumptions can be made to quantify the MSE

_{model}

: (1) shade should receive a higher design model uncertainty because the endmember selection is more difficult in shade and (2) the model misfit is small compared to the design model uncertainty. Hence, conservative zeroth-order approximations of 10% and 15% SCF LAMSU design model uncertainty (square root of MSE, Section 2.2.3) are attributed to each SCF estimation for illuminated and shaded regions, respectively, (see Equation (9)).

3.2. High-Resolution Sensor Application

Figure 6 illustrates the performance of the LAMSU method applied to Sentinel-2 in a mountainous area. LAMSU can separate snow-free and snow-covered areas in the shade, where reflectances often drop below 0.1. SCF estimates are less accurate in the shade than in the illuminated regions, which is reflected by the RMSE. The SCF results show a stark contrast to the SCFs that can be expected from the NDSI-based method in Figure 1.

Figure 7 shows the LAMSU result over Mount Shashta in North-California, U.S. The snow is nearly radially symmetrically distributed above the tree line. However, the illumination conditions are not the same around Mount Shashta. Nevertheless, the SCF map primarily indicates fully snow-covered pixels above the tree line; below the tree line, SCFs drop rapidly as the forest canopy obscures the snow underneath. Large spectral changes at short distances (at a scale of a single pixel) can be challenging for algorithms that utilize information from nearby pixels. Nevertheless, the SCFs also reveal the clear-cuts within the scene. Furthermore, LAMSU is relatively insensitive to the clouds in this scene, which are classified as snow-free. The RMSEs for these cloud-covered pixels are generally larger than the RMSEs of cloud-free pixels in the same scene.

Figure 8 displays a challenging scene with thin and broken snow cover on the Tibetan plateau. The broken snow introduces many mixed pixels, which reduces the number of potential endmembers. Most of the broken snow cover is captured by LAMSU. However, several very thinly covered patches are underestimated and detected as snow-free. The limited sensitivity originates from the endmember selection, whose snow-free endmembers are, in this case, very similar to these underestimated patches. Although, the overall RMSE estimates are comparable to the RMSE estimates in Figure 6 and Figure 7, the South-East corner of the image has relatively high RMSEs.

After further investigation, it was found that many snow-covered endmembers come from (near) the frozen river sections. These endmembers have spectral signatures that do not necessarily resemble the spectra in the snow-covered region in the South-East corner. Nevertheless, these endmembers still fit significantly better than the snow-free endmembers in the vicinity. As a consequence, the SCF is high; however, the RMSE is also high due to the propagated misfit.

Additional LAMSU results are given in the Appendix A in Figure A6, which display the region around the Zeravshan glacier and the Darvoz mountain range in Tajikistan.

3.3. Medium-Resolution Sensor Application and Intercomparison

There is virtually no limit to the number of TOA reflective bands that the LAMSU method can include for the estimation of the SCF. The same procedure (Figure 4) is applied to data from both optical sensors (SLSTR and OLCI) onboard of the Sentinel-3 A/B satellites combined, as well as from VIIRS and MODIS individually, utilizing 27, 5 and 7 reflective bands, respectively. The selected scenes are largely cloud-free and during the comparison, the clouds are manually masked out to avoid false classifications that could otherwise be introduced by the cloud screening algorithm. The result is given in Figure 9.

The cumulus clouds present in the bottom left of Figure 9a are detected as snow-free by LAMSU, again displaying the limited sensitivity to the presence of (cumulus) clouds. Accurate SCF estimates become increasingly relevant for lower ground resolutions. The consistency of the medium-resolution products is, therefore, evaluated through a pixel-wise comparison.

In addition to the intercomparison of the full medium-resolution SCF maps, the results from overlapping high-resolution scenes have also been added. The extents of these high-resolution scenes are illustrated with dashed boxes in Figure 9 and their details are listed in Table A1. The intercomparison results are shown in Table 4. The multi-sensor bias increases as the spatial resolution becomes lower but remains less than 2.2% for all comparisons. This is a consequence of the increasingly lower detection limit. Hence, low resolution SCF maps tend to estimate less snow than fine resolution SCF maps. The multi-sensor RMSEs follow the same trend, except for the comparisons against Sentinel-2 MSI. The comparisons of the individual Sentinel-2 and Landsat 8 SCF maps with the medium-resolution SCF maps are shown in Table 5 to further investigate the increased multi-sensor RMSEs for Sentinel-2 MSI.

Table 5 shows that the SCF estimates from individual Sentinel-2 scenes on 1 April 2020 display more snow than the medium-resolution SCF maps from 2 April 2020. However, SCF estimates from Sentinel-2 scenes on 3 April 2020, 2 days later, show less snow. The scenes are taken during the melt season, causing snow cover to retreat to higher elevations where it is colder. The crossover from a positive multi-sensor bias to a negative multi-sensor bias lies between 1 April 2020 and 3 April 2020, which is the date at which the medium-resolution maps were recorded. This is most obvious for entries 1 and 6 and 2 and 7, which cover the same UTM granule. The 0% bias crossover takes place between the two dates, i.e., on 2 April 2020, which indicates consistent behaviour that was not apparent from Table 4. A close-up of the European Alps is given in Figure A7 in the Appendix A, for which the spatial extent is outlined by the red box in Figure 9a.

The ability to adapt to local surfaces and environments makes the LAMSU method also applicable to global snow mapping. An example of a global SCF map is presented for 8 March 2020 from Sentinel-3A/B SLSTR and OLCI data in Figure 10. The SCF estimates over the Greenland and the Antarctic ice sheet are (near) 100%. High SCFs can also be seen at high latitudes in the Northern Hemisphere. The boreal forest belt is also visible, which partially covers the snow present under its canopy. In the Southern Hemisphere and at lower latitudes, snow is only prominent in mountain ranges, such as the Caucasus mountains and the Himalayas.

4. Discussion

LAMSU and its current endmember selection procedure is based on the principle, similar to the NDSI, that snow reflects strongly in the visible and weakly in the shortwave-infrared part of the spectrum. In addition, it was shown that shade decreases the TOA reflectance, especially for longer wavelengths. These characteristics are captured by the spectral indices in Table 2 and were subject to extensive visual inspection of spectra under snow-free and snow-covered conditions. Endmember redundancy as well as a conservative selection and weighting significantly reduce the impact of erroneously selected endmembers on the final SCF estimates.

Nevertheless, for the SCF uncertainty estimate, a zeroth-order approximation of the endmember selection uncertainty is included and propagated. High SCF uncertainties were found in shaded regions and when the selected endmembers did not fit the query spectra. The thresholds for the individual rules for endmember identification were manually selected using TOA reflectance data. Therefore, the thresholds need to be adjusted if surface reflectance data are used as input.

In regard to providing endmembers manually, we suggest to use endmembers from an earlier date or adjacent scenes with similar surface conditions. This step becomes particularly relevant when processing large regions because the Earth’s surface has vast areas without snow-free or snow-covered endmembers. Additionally, water bodies and shaded areas can have similar spectral signatures. Strong similarities can also occur for certain cloud types and snow. A cloud and water mask (e.g., [40]) is, therefore, recommended.

The overall bias (0.2% SCF) and RMSE (14.3% SCF) of the LAMSU method were derived by comparison against very-high-resolution SCF maps from WorldView-2/3 data. The highest RMSEs (up to 21.4%) were found in scenes with complex and shaded terrain. Additionally, it was found that LAMSU detects icy surfaces as (near) fully snow-covered, whereas the reference data suggests lower SCFs. In comparison to the standard FRA6T algorithm [17], the LAMSU achieves a lower bias and RMSE in all cases. This finding is in line with the behavior of NDSI-based methods in complex terrain (see Section 1.3).

The aforementioned overall bias and RMSE are of similar magnitude as validation results from other recently developed SCF retrieval methods; the evaluation of SCF estimates generated with a machine learning technique from Landsat TM/ETM+ data with three reference snow maps from Ikonos data [32] resulted in a bias less than 2.0% SCF and RMSEs in the range of 9.0% to 14.0% SCF. Ref. [30] (pre-print) found a median RMSE of 12.4% SCF for Landsat 8 OLI based SCF maps, calculated with the spectral unmixing method called Snow Property Inversion From Remote Sensing (SPIReS) [25], using snow depth measurements from airborne laser scanning surveys as ground truth information. The SCF maps from Sentinel-2 MSI data generated with the NDSI-based method of [18] showed an average bias and RMSE of −5.6% and 26.8% SCF, respectively, compared to reference snow maps from SPOT-6/-7 scenes.

The distribution of the snow-free and snow-covered validation data was not reported in these studies, although this may negatively or positively affect validation results. Future intercomparisons of the performance of SCF methods will most certainly benefit from supplementary documentation on the used sampling strategy and distribution. The SPIReS method and the LAMSU method are both based on multispectral unmixing. SPIReS additionally includes a canopy correction, grain radius, albedo and dust concentration estimate but does not, to date, provide corresponding uncertainties and is not locally adaptive.

The SPIReS method has demonstrated high accuracy at regional scales for retrievals based on Terra MODIS and Landsat 8 OLI, with lower fractional detection limits (cut-off) of 7% and 10% SCF, respectively. Overall, both the LAMSU method and the SPIReS method show improved results when a cloud and water mask is included.

The LAMSU method also automatically addresses multiple sensors in a consistent manner. LAMSU was applied to Sentinel-2A/B MSI, Landsat 8/9 OLI(-2), Sentinel-3A/B SLSTR/OLCI, Suomi NPP VIIRS and Terra MODIS TOA reflectance data. SCF estimates computed from medium-resolution data tend to estimate less snow than SCF estimates from high resolution data. However, the multi-sensor bias and RMSE among these different sensors do not exceed 2.2% and 9.7% SCF, respectively. The highest multi-sensor RMSEs occurred for a scene that was recorded one day earlier and later whilst snow melt was present. These intercomparisons demonstrated a consistency across multiple sensors that is well within the SCF RMSE from the validation.

The implementation of the software was done entirely in C [52] with OpenMP [53] for parallelization. The computation time of the Alpine scene from 27 Sentinel-3 bands (Figure 9d) was 329 ± 1 s (average of 10 runs) on a laptop (AMD Ryzen 7 PRO 5850U, 27 GB RAM).

5. Conclusions

We demonstrated that multispectral unmixing can be made locally adaptive, thereby, circumventing a spectral library. The presented methodology computes SCF estimates under varying illumination conditions and surface types, whilst maintaining reasonable processing times. LAMSU, therefore, incorporates the fact that the spectra of snow-free and snow-covered surfaces could widely vary. Furthermore, the presented methodology performed consistently for five different high- and medium-resolution satellite sensors. It is expected that the aforementioned will allow for tractable large-scale multi-sensor processing to estimate the SCF globally in near-real-time with a quantitative uncertainty.

We will further investigate and strengthen the performance of LAMSU for obtaining consistent snow products from different sensors, accounting for specific features of different regimes of the global snow cover. A particular focus will also be placed on the further development of a multi-sensor snow product, combining data from sensors that capture the Earth’s surface at different resolutions.

Furthermore, the algorithm will be expanded and tested with hyperspectral data from the German EnMAP (Environmental Mapping and Analysis Program) and the Italian PRISMA (PRecursore IperSpettrale della Missione Applicativa) mission as a contribution to the scientific preparation for the proposed Copernicus Expansion Mission CHIME (Copernicus Hyperspectral Imaging Mission for the Environment).

Author Contributions

Conceptualization, L.K. and T.N.; methodology, L.K., N.M. and G.S.; software, L.K. and M.H.; validation, L.K.; formal analysis, L.K.; investigation, L.K.; data curation, L.K.; writing—original draft preparation, L.K.; writing—review and editing, all; visualization, L.K. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the ESA Snow CCI (Contract No. 4000124098/18/I-NB) project and the ESA EXPRO+ AlpSnow - Alps Regional Initiative project (Contract No. 4000132770).

Data Availability Statement

The Sentinel-2 MSI and Sentinel-3 SLSTR/OLCI data are available on the Copernicus Open Access Hub (https://scihub.copernicus.eu/, accessed on 14 July 2022). The Landsat 8 OLI data are available on the USGS Earth Explorer (https://www.usgs.gov/tools/earthexplorer, accessed on 4 July 2022). The Terra MODIS and Suomi-NPP VIIRS data are available on (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on July 2022). The EU-DEM is available on the Copernicus Land Monitoring Service (https://land.copernicus.eu/imagery-insitu/eu-dem/eu-dem-v1.1, accessed on 1 July 2022) and the ArcticDEM is available on the Polar Geospatial Center of the University of Minnesota (https://www.pgc.umn.edu/data/arcticdem/, accessed on 6 July 2022). The ESA Land Cover CCI is available on (https://www.esa-landcover-cci.org/?q=node/164, accessed on 24 June 2022). The Copernicus Global Surface Water dataset is available on https://global-surface-water.appspot.com/download, accessed on 16 January 2023. The WorldView data are available on the WorldView ESA archive (https://earth.esa.int/eogateway/catalog/worldview-esa-archive, accessed on 21 June 2022) after the Fast Approval Registration. The LAMSU codebase is written in C and is closed source. The Voronoi diagram implementation from Mathias Westerdahl is available on GitHub (https://github.com/JCash/voronoi, accessed on 8 June 2022). The XGBoost Python module is available on GitHub (https://github.com/dmlc/xgboost/, accessed on 12 July 2022).

Acknowledgments

I wish to thank my colleagues at ENVEO IT GmbH for their insightful comments and support, in particular Helmut Rott, Ursula Fasching and Stefan Scheiblauer.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Dataset Details

Appendix A.1. LAMSU Input, Validation and Intercomparison Datasets

Table A1. List of Sentinel-2 MSI, Landsat 8 OLI, Sentinel-3 SLSTR/OLCI, Suomi-NPP VIIRS and Terra MODIS data, which were used for the inter-comparison of the LAMSU SCF estimates for the estimation of the multi-sensor bias and RMSE.

№	Scene ID	Date
1	S2A_MSIL1C_20200401T102021_N0209_R065_T32TPS_20200401T123732	1 April 2020
2	S2A_MSIL1C_20200401T102021_N0209_R065_T32TPT_20200401T123732	1 April 2020
3	S2B_MSIL1C_20200403T100549_N0209_R022_T32TNR_20200403T130157	3 April 2020
4	S2B_MSIL1C_20200403T100549_N0209_R022_T32TNS_20200403T130157	3 April 2020
5	S2B_MSIL1C_20200403T100549_N0209_R022_T32TPR_20200403T130157	3 April 2020
6	S2B_MSIL1C_20200403T100549_N0209_R022_T32TPS_20200403T130157	3 April 2020
7	S2B_MSIL1C_20200403T100549_N0209_R022_T32TPT_20200403T130157	3 April 2020
8	LC08_L1TP_190027_20200402_20200822_02_T1	2 April 2020
9	LC08_L1TP_190028_20200402_20200822_02_T1	2 April 2020
10	S3B_OL_1_EFR_20200402T094445_20200402T094745_20200403T130637_0180_037_193_2160_LN1_O_NT_002	2 April 2020
	S3B_SL_1_RBT_20200402T094445_20200402T094745_20200403T145347_0180_037_193_2160_LN2_O_NT_004
11	SVI01_npp_d20200402_t1138030_e1143434_b43685_ c20200402154344482690_nobc_ops	2 April 2020
	SVI02_npp_d20200402_t1138030_e1143434_b43685_ c20200402154344489621_nobc_ops
	SVM04_npp_d20200402_t1138030_e1143434_b43685_ c20200402154344475375_nobc_ops
	SVM10_npp_d20200402_t1138030_e1143434_b43685_ c20200402154344428499_nobc_ops
	SVM11_npp_d20200402_t1138030_e1143434_b43685_ c20200402154344434443_nobc_ops
13	MOD021KM.A2020093.1050.061.2020093191243	2 April 2020

Table A2. List of WorldView data that were used for validation of the LAMSU method.

Ref. №	WorldView Scene ID
	Date	Bands	Resolution	Location	Center Longitude, Latitude	Solar Elevation Angle
	Comparison № Landsat 8/Sentinel 2 Scene ID Snow-Free Snow
1	WV2_OPER_WV1_4B__2A_20150626T130722_N65-304_W018-337_0001_v0100
	26-06-2015	4	2 m × 2 m	Eyjafjarðarsveit, Iceland	−18.3362, 65.3035	48.2 $^{°}$
	1 LC08_L1TP_220014_20150626_20200909_02_T 174.8% 25.2%
2	WV2_OPER_WV1_4B__2A_20200827T101318_N45-537_E007-279_0001_v0100
	27-08-2020	3	0.5 m × 0.5 m	Aosta Valley, Italy	7.2800, 45.5375	50.6 $^{°}$
	1 LC08_L1TP_195028_20200827_20200905_02_T 149.0% 51.0% 2 S2B_MSIL1C_20200827T102559_N0209_R108_T32TLR 32.6% 67.4%
3	WV2_OPER_WV1_4B__2A_20210806T103439_N45-541_E007-314_0001_v0100
	06-08-2021	3	0.4 m × 0.4 m	Aosta Valley, Italy	7.3138, 45.5414	58.4 $^{°}$
	1 S2B_MSIL1C_20210809T101559_N0301_R065_T32TLR 75.5% 24.5%
4	WV3_OPER_WV1_4B__2A_20200623T104458_N45-538_E007-317_0001_v0100
	23-06-2020	3	0.3 m × 0.3 m	Aosta Valley, Italy	7.3171, 45.5388	66.0 $^{°}$
	1 LC08_L1TP_195028_20200624_20200824_02_T1 12.6% 87.4% 2 S2A_MSIL1C_20200623T103031_N0209_R108_T32TLR 12.4% 87.6%
5	WV3_OPER_WV1_8B__2A_20201122T104312_N45-409_E007-036_0001_v0100
	22-11-2020	8	1.2 m × 1.2 m	Val-d’Isère, France	7.0381, 45.4096	24.0 $^{°}$
	1 S2B_MSIL1C_20201125T103349_N0209_R108_T32TLR 39.0% 61.0%

Appendix A.2. Worldview-2/3 Classification

Figure A1. Left: WorldView-2 true colour composite taken on 26 June 2015 (reference №1 in Table A2). Right: Binary snow classification of the left-side image.

Figure A2. (Left): WorldView-2 true colour composite taken on 27 August 2020 (reference №2 in Table A2). (Right): Binary snow classification of the left-side image.

Figure A3. (Left): WorldView-2 true colour composite taken on 6 August 2021 (reference №3 in Table A2). (Right): Binary snow classification of the left-side image.

Figure A4. (Left): WorldView-3 true colour composite taken on 23 June 2020 (reference №4 in Table A2). (Right): Binary snow classification of the left-side image.

Figure A5. (Left): WorldView-3 true colour composite taken on 22 November 2020 (reference №5 in Table A2). (Right): Binary snow classification of the left-side image.

Appendix B. LAMSU SCF Examples

Figure A6. The left column shows the Zeravshan glacier between the Turkestan and Zarafshan mountain range near the Tajik–Kyrgyz border, and the right column shows the Darvoz mountain range in Tajikistan. First row: subsections from Landsat 8 OLI false colour composite (bands 6, 5 and 3) taken on 14 May 2013 (path 153, row 033). Second row: LAMSU SCF derived from Landsat 8 OLI data at 30 m resolution. Third row: LAMSU SCF RMSE at 30 m resolution with 90% transparency over the false colour composites.

Figure A7. (Top left): Sentinel-3B SLSTR false colour composite (bands 5, 3 and 1) taken on 2 April 2020. (Top right): LAMSU SCF from Terra MODIS data at 1000 m ground resolution. (Center left): LAMSU SCF from Suomi-NPP VIIRS data at 500 m ground resolution. (Center right): LAMSU SCF from Sentinel-3B SLSTR and OLCI data at 250 m ground resolution. (Bottom left): Sentinel-2A MSI false colour composite (bands 11, 8A and 3) taken on 1 April 2020. (Bottom right): LAMSU SCF from Sentinel-2A MSI data at 20 m ground resolution. The top and middle row are from 2 April 2020. The bottom row is from 1 April 2020.

References

Estilow, T.W.; Young, A.H.; Robinson, D.A. A long-term Northern Hemisphere snow cover extent data record for climate studies and monitoring. Earth Syst. Sci. Data 2015, 7, 137–142. [Google Scholar] [CrossRef] [Green Version]
Barnett, T.; Adam, J.; Lettenmaier, D. Potential Impacts of a Warming Climate on Water Availability in Snow-Dominated Regions. Nature 2005, 438, 303–309. [Google Scholar] [CrossRef]
Vikhamar, D.; Solberg, R. Snow-cover mapping in forests by constrained linear spectral unmixing of MODIS data. Remote Sens. Environ. 2003, 88, 309–323. [Google Scholar] [CrossRef]
Dankers, R.; Feyen, L. Flood hazard in Europe in an ensemble of regional climate scenarios. J. Geophys. Res. Atmos. 2009, 114. [Google Scholar] [CrossRef]
Krøgli, I.K.; Devoli, G.; Colleuille, H.; Boje, S.; Sund, M.; Engen, I.K. The Norwegian forecasting and warning service for rainfall- and snowmelt-induced landslides. Nat. Hazards Earth Syst. Sci. 2018, 18, 1427–1450. [Google Scholar] [CrossRef] [Green Version]
Robinson, D.A.; Estilow, T.W. Rutgers Northern Hemisphere 24 km Weekly Snow Cover Extent, September 1980 Onward, Version 1. 2021. Available online: https://nsidc.org/data/g10035/versions/1 (accessed on 12 July 2022).
U.S. National Ice Center. IMS Daily Northern Hemisphere Snow and Ice Analysis at 1 km, 4 km, and 24 km Resolutions, Version 1. 2008. Available online: https://nsidc.org/data/g02156/versions/1 (accessed on 11 July 2022).
Hori, M.; Sugiura, K.; Kobayashi, K.; Aoki, T.; Tanikawa, T.; Kuchiki, K.; Niwano, M.; Enomoto, H. A 38-year (1978–2015) Northern Hemisphere daily snow cover extent product derived using consistent objective criteria from satellite-borne optical sensors. Remote Sens. Environ. 2017, 191, 402–418. [Google Scholar] [CrossRef]
Nagler, T.; Schwaizer, G.; Mölg, N.; Keuris, L.; Hetzenecker, M.; Metsämäki, S. ESA Snow Climate Change Initiative (Snow_cci): Daily global Snow Cover Fraction - viewable snow (SCFV) from MODIS (2000–2020). 2022. Available online: https://catalogue.ceda.ac.uk/uuid/ebe625b6f77945a68bda0ab7c78dd76b (accessed on 1 July 2022).
Romanov, P. Global Multisensor Automated satellite-based Snow and Ice Mapping System (GMASI) for cryosphere monitoring. Remote Sens. Environ. 2017, 196, 42–55. [Google Scholar] [CrossRef]
Riggs, G.A.; Hall, D.K.; Tschudi, M.; NASA VIIRS Land SIPS. VIIRS/NPP Snow Cover Daily L3 Global 375m SIN Grid, Version 1. 2019. Available online: https://nsidc.org/data/vnp10a1/versions/1 (accessed on 30 June 2022).
Hall, D.K.; Riggs, G.A. MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid. 2016. Available online: https://nsidc.org/data/mod10a1/versions/6 (accessed on 27 June 2022).
Solberg, R.; Trier, ø.; Breivik, L.A.; Godøy, ø.; Killie, M.; Andreassen, L.M.; Hausberg, J.; Olsen, O. CryoClim—A New System for Cryospheric Climate Monitoring. In Proceedings of the 33rd International Symposium on Remote Sensing of Environment, ISRSE 2009, Stresa, Italy, 4–8 May 2009. [Google Scholar]
Siljamo, N.; Hyvärinen, O.; Riihelä, A.; Suomalainen, M. MetOp/AVHRR Snow Detection Method for Meteorological Applications. J. Appl. Meteorol. Climatol. 2020, 59, 2001–2019. [Google Scholar] [CrossRef]
Dozier, J. Spectral signature of alpine snow cover from the landsat thematic mapper. Remote Sens. Environ. 1989, 28, 9–22. [Google Scholar] [CrossRef]
Salomonson, V.; Appel, I. Estimating fractional snow cover from MODIS using the normalized difference snow index. Remote Sens. Environ. 2004, 89, 351–360. [Google Scholar] [CrossRef]
Salomonson, V.; Appel, I. Development of the Aqua MODIS NDSI fractional snow cover algorithm and validation results. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1747–1756. [Google Scholar] [CrossRef]
Gascoin, S.; Barrou Dumont, Z.; Deschamps-Berger, C.; Marti, F.; Salgues, G.; López-Moreno, J.I.; Revuelto, J.; Michon, T.; Schattan, P.; Hagolle, O. Estimating Fractional Snow Cover in Open Terrain from Sentinel-2 Using the Normalized Difference Snow Index. Remote Sens. 2020, 12, 2904. [Google Scholar] [CrossRef]
Riggs, G.; Hall, D.; Román, M.O. MODIS Snow Products Collection 6.1 User Guide. 2019. Available online: https://nsidc.org/sites/default/files/c61_modis_snow_user_guide.pdf (accessed on 24 August 2022).
Riggs, G.; Hall, D.; Vuyovich, C.; DiGirolamo, N. Development of Snow Cover Frequency Maps from MODIS Snow Cover Products. Remote Sens. 2022, 14, 5661. [Google Scholar] [CrossRef]
Jing, Y.; Li, X.; Shen, H. STAR NDSI collection: A cloud-free MODIS NDSI dataset (2001–2020) for China. Earth Syst. Sci. Data 2022, 14, 3137–3156. [Google Scholar] [CrossRef]
HR-S&I Consortium. Algorithm Theoretical Basis Document for Snow Products Based on Sentinel-2 (COSIMS-DT-062-MAG_ATBD_SNOW). 2022. Available online: https://land.copernicus.eu/user-corner/technical-library/hrsi-snow-atbd (accessed on 28 June 2022).
Painter, T.H.; Dozier, J.; Roberts, D.A.; Davis, R.E.; Green, R.O. Retrieval of subpixel snow-covered area and grain size from imaging spectrometer data. Remote Sens. Environ. 2003, 85, 64–77. [Google Scholar] [CrossRef] [Green Version]
Painter, T.H.; Rittger, K.; McKenzie, C.; Slaughter, P.; Davis, R.E.; Dozier, J. Retrieval of subpixel snow covered area, grain size, and albedo from MODIS. Remote Sens. Environ. 2009, 113, 868–879. [Google Scholar] [CrossRef] [Green Version]
Bair, E.H.; Stillinger, T.; Dozier, J. Snow Property Inversion From Remote Sensing (SPIReS): A Generalized Multispectral Unmixing Approach With Examples From MODIS and Landsat 8 OLI. IEEE Trans. Geosci. Remote Sens. 2021, 59, 7270–7284. [Google Scholar] [CrossRef]
Dietz, A.J.; Kuenzer, C.; Gessner, U.; Dech, S. Remote sensing of snow—A review of available methods. Int. J. Remote Sens. 2012, 33, 4094–4134. [Google Scholar] [CrossRef]
Rittger, K.; Painter, T.H.; Dozier, J. Assessment of methods for mapping snow cover from MODIS. Adv. Water Resour. 2013, 51, 367–380. [Google Scholar] [CrossRef]
Masson, T.; Dumont, M.; Mura, M.D.; Sirguey, P.; Gascoin, S.; Dedieu, J.P.; Chanussot, J. An Assessment of Existing Methodologies to Retrieve Snow Cover Fraction from MODIS Data. Remote Sens. 2018, 10, 619. [Google Scholar] [CrossRef] [Green Version]
Aalstad, K.; Westermann, S.; Bertino, L. Evaluating satellite retrieved fractional snow-covered area at a high-Arctic site using terrestrial photography. Remote Sens. Environ. 2020, 239, 111618. [Google Scholar] [CrossRef]
Stillinger, T.; Rittger, K.; Raleigh, M.S.; Michell, A.; Davis, R.E.; Bair, E.H. Landsat, MODIS, and VIIRS snow cover mapping algorithm performance as validated by airborne lidar datasets. Cryosphere Discuss. 2022, 2022. [Google Scholar] [CrossRef]
Dobreva, I.D.; Klein, A.G. Fractional snow cover mapping through artificial neural network analysis of MODIS surface reflectance. Remote Sens. Environ. 2011, 115, 3355–3366. [Google Scholar] [CrossRef] [Green Version]
Czyzowska-Wisniewski, E.H.; van Leeuwen, W.J.; Hirschboeck, K.; Marsh, S.E.; Wisniewski, W.T. Fractional snow cover estimation in complex alpine-forested environments using an artificial neural network. Remote Sens. Environ. 2015, 156, 403–417. [Google Scholar] [CrossRef]
Hofmeister, F.; Arias-Rodriguez, L.; Premier, V.; Marin, C.; Notarnicola, C.; Disse, M.; Chiogna, G. Intercomparison of Sentinel-2 and modelled snow cover maps in a high-elevation Alpine catchment. J. Hydrol. X 2022, 15, 100123. [Google Scholar] [CrossRef]
Barella, R.; Marin, C.; Gianinetto, M.; Notarnicola, C. A Novel Approach to High Resolution Snow Cover Fraction Retrieval in Mountainous Regions. In Proceedings of the IGARSS 2022 - 2022 IEEE International Geoscience and Remote Sensing Symposium, Kuala Lumpur, Malaysia, 17–22 July 2022; pp. 3856–3859. [Google Scholar] [CrossRef]
Chen, Y.; Hall, A.; Liou, K.N. Application of three-dimensional solar radiative transfer to mountains. J. Geophys. Res. Atmos. 2006, 111. [Google Scholar] [CrossRef]
Sirguey, P.; Mathieu, R.; 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. 2009, 113, 160–181. [Google Scholar] [CrossRef]
Kuter, S.; Bolat, K.; Akyurek, Z. A machine learning-based accuracy enhancement on EUMETSAT H-SAF H35 effective snow-covered area product. Remote Sens. Environ. 2022, 272, 112947. [Google Scholar] [CrossRef]
Nolin, A.W.; Dozier, J.; Mertes, L.A.K. Mapping alpine snow using a spectral mixture modeling technique. Ann. Glaciol. 1993, 17, 121–124. [Google Scholar] [CrossRef] [Green Version]
Cannistra, A.F.; Shean, D.E.; Cristea, N.C. High-resolution CubeSat imagery and machine learning for detailed snow-covered area. Remote Sens. Environ. 2021, 258, 112399. [Google Scholar] [CrossRef]
Defourny, P.; Lamarche, C.; Bontemps, S.; De Maet, T.; Van Bogaert, E.; Moreau, I.; Brockmann, C.; Boettcher, M.; Kirches, G.; Wevers, J.; et al. Land Cover CCI Product User Guide Version 2. Technical Report. 2017. Available online: https://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf (accessed on 18 April 2022).
Pekel, J.F.; Cottam, A.; Gorelick, N.; Belward, A.S. High-resolution mapping of global surface water and its long-term changes. Nature 2016, 540, 418–422. [Google Scholar] [CrossRef] [PubMed]
Porter, C.; Morin, P.; Howat, I.; Noh, M.J.; Bates, B.; Peterman, K.; Keesey, S.; Schlenk, M.; Gardiner, J.; Tomko, K.; et al. ArcticDEM. 2018. Available online: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/OHHUKH (accessed on 18 April 2022).
González, J.C.G.; Redondo, J.A.; Garzón, A. EU-DEM Upgrade: Documentation EEA User Manual (C4EO17 1/2). 2015. Available online: https://land.copernicus.eu/user-corner/technical-library/eu-dem-v1-1-user-guide (accessed on 9 May 2022).
Nagler, T.; Ripper, E.; Rott, H. Preparations for Snow Monitoring using Sentinel-3 SLSTR and OLCI. 2012. Available online: https://www.researchgate.net/publication/258776433_Preparation_for_Snow_Cover_Monitoring_Using_Sentinel-1_and_Sentinel-3_Data (accessed on 3 May 2022).
Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring Vegetation Systems in the Great Plains with ERTS. 1973. Available online: https://ntrs.nasa.gov/citations/19740022614 (accessed on 30 August 2022).
Jiang, Z.; Huete, A.R.; Didan, K.; Miura, T. Development of a two-band enhanced vegetation index without a blue band. Remote Sens. Environ. 2008, 112, 3833–3845. [Google Scholar] [CrossRef]
Chang, C.I. Spectral information divergence for hyperspectral image analysis. In Proceedings of the IEEE 1999 International Geoscience and Remote Sensing Symposium. IGARSS’99, Hamburg, Germany, 8 June–2 July 1999; Volume 1, pp. 509–511. [Google Scholar] [CrossRef]
Fortune, S. A sweepline algorithm for Voronoi diagrams. Algorithmica 1987, 2. [Google Scholar] [CrossRef]
Stark, P.; Parker, R. Bounded-variable least-squares: An algorithm and applications. Comput. Stat. 1995, 129–141. [Google Scholar]
Chen, T.; Guestrin, C. XGBoost. In Proceedings of the Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, 13–17 August 2016. [Google Scholar] [CrossRef] [Green Version]
Hastie, T.; Tibshirani, R.; Friedman, J.H. Boosting and Additive Trees. In The Elements of Statistical Learning, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2009; pp. 337–384. [Google Scholar]
Kernighan, B.W.; Ritchie, D.M. The C Programming Language; Prentice Hall: Englewood Cliffs, NJ, USA, 2006. [Google Scholar]
OpenMP Architecture Review Board. OpenMP Application Program Interface Version 4.5. 2015. Available online: https://www.openmp.org/wp-content/uploads/openmp-4.5.pdf (accessed on 20 January 2022).

Figure 1. (Left): Sentinel-2B false colour composite (bands 11, 8A and 3) taken on 25 November 2020 at 10:33 AM UTC in the Chamonix valley in the European Alps (32TLR) with the RGB colours scaled from 0 to 0.25, 0 to 0.15 and 0 to 0.15, respectively. The round and diamond markers are placed in illuminated and shaded areas, respectively. The solar elevation angle is 26.5

^{°}

. (Right): NDSI value along the profile line in the left image. The line colours indicate the snow-free (green) and snow-covered (pink) pixels along the profile. The shaded region is indicated with a grey background and the NDSI values for 80% and 100% SCF using FRA6T [17] are indicated with dashed and solid horizontal lines, respectively.

^{°}

Figure 2. Both figures show three green curves for the snow-free spectral signatures and three purple curves for the snow-covered spectral signatures. The values are the Sentinel-2B MSI TOA reflectances from Figure 1. The (left) hand figure displays the reflectances for the well-illuminated spectra (round markers in Figure 1), and the (right) hand figure displays the reflectance for the shaded spectra (diamond markers in Figure 1). Please note the different scales on the y-axis.

Figure 3. The reflective bands within the electromagnetic spectrum listed per sensor that were used within this paper as input data by LAMSU. The total number of reflective bands is given on the far right. The bands’ approximate ground resolution that is available for download is indicated by the second y-axis.

Figure 4. Proposed novel adaptive SCF estimation method with LAMSU at the core of the processing scheme.

Figure 5. Workflow for the selection of local endmembers within a multispectral image. The spectral indices are described in Table 2.

Figure 6. (Left): Sentinel-2B MSI false colour composite (bands 11, 8A and 3) taken on 25 November 2020 in the European Alps near Chamonix (32TLR) with the RGB colours scaled from 0 to 0.25, 0 to 0.15 and 0 to 0.15, respectively. (Center): LAMSU SCF derived from Sentinel-2B MSI data at 20 m resolution. Right: LAMSU SCF RMSE at 20 m resolution. Water mask source: EC JRC/Google.

Figure 7. (Left): Landsat 8 OLI false colour composite (band 6, 5 and 3) over Mount Shasta in North-California in the U.S. taken on 2 April 2017 (path 045, row 031). (Center): LAMSU SCF derived from Landsat 8 OLI data at 30 m resolution. (Right): LAMSU SCF RMSE at 30 m resolution. Water mask source: EC JRC/Google.

Figure 8. (Left): Landsat 9 OLI-2 false colour composite (bands 6, 5 and 3) near the Northern boundary of the Tibetan plateau (South of the Kunlun mountains) taken on 18 January 2022 (path 140, row 034). (Center): LAMSU SCF derived from Landsat 9 OLI-2 data at 30 m resolution. (Right): LAMSU SCF RMSE at 30 m resolution. Water mask source: EC JRC/Google.

Figure 9. Example over the Alps on 2 April 2020 from (a) Sentinel-3B SLSTR false colour composite (bands 5, 3 and 2). (b) LAMSU SCF derived from Terra MODIS data at 1000 m. (c) LAMSU SCF derived from Suomi NPP VIIRS data at 500 m. (d) LAMSU SCF derived from Sentinel-3B SLSTR and OLCI data combined at 250 m. An enlarged section is presented in Figure A7, which covers the red box in subfigure (a). The yellow and purple dashed boxes in subfigure (a) illustrate the extent of the comparison with high-resolution snow maps from Table 4 and Table 5.

Figure 10. LAMSU SCF derived from Sentinel-3A/B SLSTR and OLCI data at 250 m resolution for 8 March 2020.

Table 1. Summary of band details of the input data with their native resolutions used in this study. The sensor abbreviations are MultiSpectral Instrument (MSI), Operational Land Imager (OLI), Ocean and Land Colour Instrument (OLCI), Sea and Land Surface Temperature Radiometer (SLSTR), Visible Infrared Imaging Radiometer Suite (VIIRS), MODerate Resolution Imaging Spectroradiometer (MODIS) and WorldView-110 (WV110).

Satellite	Sensor	Bands
Sentinel-2 A/B	MSI	2, 3, 4, 8 (10 m) 5, 6, 7, 8A, 11, 12 (20 m) 1, 9, 10 (60 m)
Landsat 8/9	OLI(-2)	8 (15 m) 1, 2, 3, 4, 5, 6, 7, 9 (30 m)
Sentinel-3 A/B	OLCI	Oa01, Oa02, Oa03, Oa04, Oa05, Oa06, Oa07, Oa08, Oa09, Oa10, Oa11, Oa12, Oa13, Oa14, Oa15, Oa16, Oa17, Oa18, Oa19, Oa20, Oa21 (300 m)
Sentinel-3 A/B	SLSTR	S1, S2, S3, S4, S5, S6 (500 m)
Suomi NPP	VIIRS	i01, i02 (375 m) m04, m10, m11 (750 m)
Terra	MODIS	1, 2 (250 m) 3, 4, 5, 6, 7 (500 m) (all bands are also available as aggregated 1 km products)
WorldView-2	WV110	panchromatic (1) (0.41 m) multispectral (8) (1.64 m)
WorldView-3	WV110	panchromatic (1) (0.31 m) multispectral (8) (1.24 m)

Table 2. The spectral indices that are used for the local endmember selection. Reflectances of five spectral bands are used:

r_{560}

r_{650}

r_{860}

r_{1610}

and

r_{2200}

, which represent the reflectance of green, red, NIR and two types of SWIR light. The subscript indicates their wavelengths in nanometers.

Table 2. The spectral indices that are used for the local endmember selection. Reflectances of five spectral bands are used:

r_{560}

r_{650}

r_{860}

r_{1610}

and

r_{2200}

, which represent the reflectance of green, red, NIR and two types of SWIR light. The subscript indicates their wavelengths in nanometers.

№	Spectral Index	Characteristic
1	$\frac{r_{560} - r_{1610}}{r_{560} + r_{1610}}$	Increases with snow presence and shade. This index is also known as the NDSI [15].
2	$\frac{r_{860} - r_{650}}{r_{860} + r_{650}}$	Increases with vegetation presence and decreases within shade. This index is also known as the NDVI [45].
3	$\frac{2.5 \cdot (r_{860} - r_{650})}{r_{860} + 2.4 \cdot r_{650} + 1}$	Increases with vegetation presence and decreases within shade. This index is also known as the two-band EVI [46].
4	$2 \frac{r_{560} + r_{660} - r_{860} - r_{1610}}{r_{560} + r_{660} + r_{860} + r_{1610}}$	Increases with snow presence but is more robust and sensitive to snow in the shade than NDSI.
5	$2 \frac{2 \cdot r_{560} - r_{660} - r_{1610}}{2 \cdot r_{560} + r_{660} + r_{1610}}$	Same as (4) but is less sensitive to vegetation.
6	$2 \frac{2 \cdot r_{560} - r_{650} - r_{850}}{2 \cdot r_{560} + r_{650} + r_{850}}$	Same as (4) but is more sensitive to shade than snow.
7	$2 \frac{2 \cdot r_{650} - r_{860} - r_{1610}}{2 \cdot r_{650} + r_{860} + r_{1610}}$	Same as (4) but is less sensitive to snow.
8	$\frac{r_{560} - r_{650}}{r_{560} + r_{650}}$	Increases with shade and is close to zero for snow.
9	$\frac{r_{650} - r_{1610}}{r_{650} + r_{1610}}$	Increases with snow presence and shade. Slightly lower sensitivity to snow presence and shade than (1).
10	$\frac{r_{860} - r_{1610}}{r_{860} + r_{1610}}$	Increases with snow presence and shade. Lower sensitivity to snow presence and shade (in particular in illuminated areas) than (1).
11	$\frac{r_{860} - r_{2200}}{r_{860} + r_{2200}}$	Increases with snow presence and shade. Lower sensitivity to snow presence and shade (in particular in illuminated areas) than (1).

Table 3. LAMSU SCF validation results. Bias is the average SCF bias (LAMSU SCF or FRA6T [17] SCF minus reference SCF) and RMSE is the average SCF root-mean-square error (both) from one million realisations. N is the number of valid pixel-wise comparisons from which the subsets are drawn. More details on the reference data can be found in Table A2 in the Appendix A. L8 and S2 refer to the Landsat 8 OLI (30 m) and Sentinel-2 MSI (20 m) SCF products, respectively.

		LAMSU		FRA6T [17]
Reference №	L8/S2	Bias (%)	RMSE (%)	Bias (%)	RMSE (%)	N (#)
1	L8	−3.30	10.80	4.12	13.54	29652
2	L8	8.10	21.45	17.37	35.20	18203
	S2	3.85	18.56	7.92	25.09	12327
3	S2	0.62	16.79	8.36	26.21	102602
4	L8	0.42	9.30	3.37	13.18	13115
	S2	0.87	10.10	4.33	14.38	45732
5	S2	−1.95	14.17	9.87	27.59	84644
Overall		0.15	14.28	7.21	23.48	306275

Table 4. Comparison of SCF maps using the LAMSU method for different medium-resolution sensors. The lower and upper triangle show the multi-sensor SCF bias (column minus row) and multi-sensor RMSE in percent SCF, respectively. Further details can be found in Table A1. SCF maps from Sentinel-2 MSI and Landsat 8 OLI are calculated at 25 m, Sentinel-3 SLSTR/OLCI at 250 m, Suomi NPP VIIRS at 500 m and Terra MODIS at 1000 m.

	Sentinel-2 MSI	Landsat 8 OLI	Sentinel-3 SLSTR/OLCI	Suomi NPP VIIRS	Terra MODIS
Bias	Sentinel-2 MSI	Landsat 8 OLI	Sentinel-3 SLSTR/OLCI	Suomi NPP VIIRS	Terra MODIS
Sentinel-2 MSI	-	-	9.68	7.35	6.99
Landsat 8 OLI	-	-	4.92	4.33	5.32
Sentinel-3 SLSTR/OLCI	−1.46	0.76	-	8.19	8.93
Suomi NPP VIIRS	−0.34	1.28	0.65	-	8.86
Terra MODIS	0.47	1.71	2.16	2.15	-

Table 5. Comparison of SCF maps using the LAMSU method for different medium-resolution sensors against individual SCF maps from Sentinel-2 MSI data using the LAMSU method. The two columns for each sensor show the multi-sensor SCF bias (row minus column) and multi-sensor RMSE in percent SCF, respectively. The high-resolution granules are identified by their date and UTM tile for Sentinel-2 MSI and path/row for Landsat 8 OLI. More details can be found in Table A1.

№	Date	Tile	Bias	RMSE	Bias	RMSE	Bias	RMSE
Sentinel-2 MSI			Sentinel-3 SLSTR/OLCI		Suomi NPP VIIRS		Terra MODIS
1	1 April 2020	32TPS	0.59	7.15	0.68	6.72	1.19	7.59
2	1 April 2020	32TPT	2.21	8.63	1.97	7.73	3.36	8.35
3	3 April 2020	32TNR	−6.00	14.83	−1.97	7.82	−1.61	6.23
4	3 April 2020	32TNS	−1.93	10.06	−1.34	7.95	−0.29	8.37
5	3 April 2020	32TPR	−3.07	10.60	0.08	6.83	0.42	5.68
6	3 April 2020	32TPS	−1.53	8.62	−0.86	6.98	−0.10	7.67
7	3 April 2020	32TPT	−0.47	7.88	−1.16	6.57	0.30	5.07
Landsat 8 OLI			Sentinel-3 SLSTR/OLCI		Suomi NPP VIIRS		Terra MODIS
1	2 April 2020	190027	0.42	5.24	0.98	4.60	1.68	5.81
2	2 April 2020	190028	1.09	4.60	1.58	4.06	1.73	4.82

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Keuris, L.; Hetzenecker, M.; Nagler, T.; Mölg, N.; Schwaizer, G. An Adaptive Method for the Estimation of Snow-Covered Fraction with Error Propagation for Applications from Local to Global Scales. Remote Sens. 2023, 15, 1231. https://doi.org/10.3390/rs15051231

AMA Style

Keuris L, Hetzenecker M, Nagler T, Mölg N, Schwaizer G. An Adaptive Method for the Estimation of Snow-Covered Fraction with Error Propagation for Applications from Local to Global Scales. Remote Sensing. 2023; 15(5):1231. https://doi.org/10.3390/rs15051231

Chicago/Turabian Style

Keuris, Lars, Markus Hetzenecker, Thomas Nagler, Nico Mölg, and Gabriele Schwaizer. 2023. "An Adaptive Method for the Estimation of Snow-Covered Fraction with Error Propagation for Applications from Local to Global Scales" Remote Sensing 15, no. 5: 1231. https://doi.org/10.3390/rs15051231

APA Style

Keuris, L., Hetzenecker, M., Nagler, T., Mölg, N., & Schwaizer, G. (2023). An Adaptive Method for the Estimation of Snow-Covered Fraction with Error Propagation for Applications from Local to Global Scales. Remote Sensing, 15(5), 1231. https://doi.org/10.3390/rs15051231

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

An Adaptive Method for the Estimation of Snow-Covered Fraction with Error Propagation for Applications from Local to Global Scales

Abstract

1. Introduction

1.1. Background

1.2. Previous Work

1.3. Limitations of Existing Methods

1.4. Objective

2. Data and Methodology

2.1. Input Data

2.2. Methodology

2.2.1. Pre-Processing

2.2.2. Local Endmember Selection

2.2.3. Locally Adaptive MultiSpectral Unmixing

2.2.4. Post-Processing

2.3. Validation of LAMSU SCF Estimates

3. Results

3.1. Validation with WorldView-2/3 Imagery

3.2. High-Resolution Sensor Application

3.3. Medium-Resolution Sensor Application and Intercomparison

4. Discussion

5. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Appendix A. Dataset Details

Appendix A.1. LAMSU Input, Validation and Intercomparison Datasets

Appendix A.2. Worldview-2/3 Classification

Appendix B. LAMSU SCF Examples

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI