Development and Evaluation of a Cloud-Gap-Filled MODIS Normalized Difference Snow Index Product over High Mountain Asia

Abstract

Accurate snow cover data are critical for understanding the Earth’s climate system, and exploring hydrological processes and regional water resource management over High Mountain Asia (HMA). However, satellite-based remote sensing observations of snow cover have inevitable data gaps originating from cloud cover, sensor, orbital limitations and other factors. Here an effective cloud-gap-filled (CGF) method was developed to fully fill the data gaps in Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference snow index (NDSI) product. The CGF method combines the respective strengths of the cubic spline interpolation method and the spatio-temporal weighted method for generating the CGF Terra-Aqua MODIS NDSI product over HMA from 2000 to 2021. Based on the validation results of in situ snow-depth observations, the CGF NDSI product achieves a high range overall accuracy (OA) of 93.54–98.08%, a low range underestimation error (MU) of 0.15–3.49% and an acceptable range overestimation error (MO) of 0.84–5.77%. Based on the validation results of high-resolution Landsat images, this product achieves the OA of 88.52–92.40%, the omission error (OE) of 1.42–10.28% and the commission error (CE) of 5.97–17.58%. The CGF MODIS NDSI product can provide scientific support for eco-environment sustainable management in the high mountain region.

1. Introduction

Snow cover has been identified as an essential geophysical parameter for understanding the Earth’s climate system, covering about 40% to 50% of the Northern Hemisphere during winter [1,2,3]. Snow cover exerts a strong control on the surface energy budget, water cycle, primary productivity and surface gas exchange [4,5,6,7]. As a natural solid water reservoir within the cryosphere, snow cover holds a crucial role as the source of water supply, benefiting over 17% of the global population with fresh drinking water through seasonal snowmelt-driven runoff [8,9,10].

High Mountain Asia (HMA) is the largest expanse of snow cover outside of the polar regions and feeds several large rivers, often referred to as the “Earth’s third pole” and the “Water Tower of Asia” [11,12]. Since the 1970s, HMA has undergone drastic environmental changes, including the air temperature rate more than double the global average [13]. The warming climate may be accountable for accelerated melting of glaciers and snow cover over HMA and surrounding, which could have a significant impact on water supply and water security in this region. Given these concerns, it is crucial to monitor the variability of snow cover over HMA, while the long-term and highly accurate snow cover products are fundamental for snow cover research over HMA.

Compared to the limited field surveys and the sparse, uneven distribution of available snowfall observation data, satellite remote sensing data exhibit great advantages in depicting the long-term spatio-temporal patterns of snow cover in large-scale and rugged terrain regions (like HMA). The rapid development of remote sensing technology during the last decades has generated a series of snow cover products including optical-based and microwave-based (mainly derived from passive sensors) products. Microwave-based snow products (e.g., snow depth (SD) and snow water equivalent) derived from passive sensors have been proven to be useful for monitoring snow cover [14,15,16,17,18]. However, their coarse (∼25 km) footprints and the saturation of observations in deep snow (>0.8 m depth) greatly limit the applicability for mountainous areas [19,20,21]. In comparison, optical-based sensors can also provide snow cover data, in terms of more appropriate temporal and spatial resolutions by making full use of snow’s unique reflective properties in the visible and shortwave infrared bands [18,22,23,24,25]. Currently, a variety of long-time series optical-based snow cover products with different temporal and spatial resolutions have been widely used, such as Advanced Very High Resolution Radiometer (AVHRR) Global Area Coverage (GAC) daily snow cover fraction products [26], the daily Visible Infrared Imaging Radiometer Suite (VIIRS) snow cover products [27] and the Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover products [28]. Among them, the daily MODIS snow cover products from the Terra (MOD10A1) and Aqua (MYD10A1) satellites are widely used to depict the spatio-temporal patterns of snow cover, due to the advantages of a long-time series (available since 2000), high spatio-temporal resolution (i.e., 500 m and daily), global coverage and being freely available. Previous studies have demonstrated that MODIS snow cover products perform well in extracting snow cover information, exhibiting an overall accuracy exceeding 90% under clear sky conditions, despite some remaining uncertainties in forest coverage and mountainous areas [25,29,30,31,32,33,34]. Nevertheless, cloud contamination in MODIS snow cover products often causes numerous data gaps, which greatly limits their applications [35,36,37].

To reduce or eliminate the effects of cloud data gaps in the MODIS snow cover products, various methods have been developed over the past decade. Traditional cloud removal algorithms can be categorized into four types: temporal, spatial, multi-source fusion and spatio-temporal combination methods [18]. In addition, these methods have mainly been aimed at binary snow cover (BSC) or fractional snow cover (FSC) data from the MODIS snow cover product Collection 5 (C5) [18]. Temporal methods include the Terra and Aqua combination (TAC) [38], adjacent temporal deduction [39,40], multi-day combination [38,41] (e.g., the 8-day composite products of MOD/MYD10A2 product), season filter [39,42] and temporal interpolation using a mathematical function [31,43,44,45]. These temporal methods utilize the instability of cloud cover and the temporal correlation of snow cover to effectively reduce cloud cover, either partially or completely, with high overall accuracy. However, they have the problems of sacrificing temporal resolution, failing to reduce the cloud cover completely and having low accuracy when cloud coverage is continuous. Spatial methods include the spatial neighborhood filter [42], the snowline mapping approach [46] and spatial interpolation based on a regression function [47]. These spatial methods utilize the spatial distribution characteristics of snow cover for reclassifying cloud pixels. However, they have the problems of high computational complexity, and having low accuracy in a region with high cloud coverage for all spatial methods. Multi-source fusion methods, including the optical/microwave observations fusion, optical/meteorological station observations fusion and optical/microwave/meteorological station observations fusion [16,48,49,50]. These multi-source fusion methods have the capability to eliminate all cloud pixels, while the accuracy depends on the complementarity and precision of the input data [18]. However, they often sacrifice spatial resolution and result in varying degrees of uncertainty. Considering all the appearing problems from the above three methods, the spatio-temporal combination methods can achieve satisfactory cloud removal effectiveness and accuracy. Spatio-temporal combination methods rely on the correlations of snow cover in space and time, typically taking two basic forms [1]. One is to take advantage of multiple spatial methods and temporal methods step by step [38,40,42,47,51]. The other is to utilize the spatial and temporal information simultaneously [52,53,54,55]. Between the two forms of methods, the latter has shown more promise in recent years because it uses spatio-temporal information from the snow cover and removes all cloud contamination in one step. For example, Li et al. [54] developed an adaptive spatio-temporal weighted method to estimate cloud pixels by considering the probability of snow cover; Hou et al. [52] introduced a gap-filling method that utilizes a non-local spatio-temporal filter method to eliminate cloud contamination for daily MODIS FSC products; Huang et al. [53] proposed a spatio-temporal model based on hidden Markov random fields that integrates spatio-temporal-spectral information along with environmental relationships. However, these methods suffer from issues related to high computational complexity and excessively high time costs.

Currently the version of MODIS snow cover product has moved from C5 to Collection 6 (C6), because the forward processing of MODIS product C5 has been discontinued in 2016. In addition, since the NDSI from C6 is a more accurate description index of the snow detection as compared to the FSC from C5 [56], more studies are shifting to the use of MODIS C6 products. Research works have demonstrated that the MODIS snow cover product C6 has high accuracy in the Tibetan Plateau, with Terra product C6 being comparable to Terra product C5, and Aqua product C6 truly having better accuracy than Aqua product C5 [57,58,59]. However, the data gaps from cloud contamination in C6 still exist, and thus a cloud removal algorithm for MODIS NDSI products is necessary. Given the advantages and disadvantages of the above cloud removal methods, we propose a cloud-gap-filled (CGF) method to retrieve the missing NDSI information beneath the cloud gaps in the MODIS NDSI product. The CGF method can effectively combine the temporal interpolation method with the spatio-temporal weighted method, leveraging their individual strengths. The ultimate goal is to develop a long-term daily cloud-free MODIS NDSI product for HMA with high precision. This product will facilitate climate and glacio-hydrological modelling and understanding of the present dynamics of the cryosphere in the region. The abbreviations for key terminology mentioned in this paper are shown in

Table A1. Summary table of abbreviations of key terminology in this paper.

No.	Abbreviation	Full Terminology
1	BSC	Binary snow cover
2	CGF	Cloud-gap-filled
3	CPD	Cloud persistence days
4	CE	Commission error
5	CSI	Cubic spline interpolation
6	HMA	High Mountain Asia
7	LI	Linear interpolation
8	MAE	Mean absolute error
9	MODIS	Moderate Resolution Imaging Spectroradiometer
10	NDSI	Normalized difference snow index
11	OE	Omission error
12	OA	Overall accuracy
13	MO	Overestimation error
14	QI	Quadratic interpolation
15	RMSE	Root–mean–square–error
16	SCD	Snow-covered days
17	SCE	Snow-covered extent
18	SD	Snow depth
19	STW	Spatio-temporal weighted
20	TAC	Terra and Aqua combination
21	TI	Temporal interpolation
22	MU	Underestimation error

.

2. Study Area and Data

2.1. Study Area

The HMA refers to the vast high-altitude geographical area in Central Asia (spanning 22–47°N and 64–107°E), encompassing the Asian mountain ranges (e.g., Tien Shan, Pamir, Karakoram, Kunlun Mountain, and the Himalayas) surrounding the Tibetan Plateau (Figure 1). Fifteen subregions of HMA according to the Randolph Glacier Inventory version 6.0 (RGI v6.0) are shown in Figure 1 [60]. With an average elevation exceeding 4000 m, the HMA constitutes the largest expanse of snow cover outside of the Earth’s polar regions. Many rivers (e.g., Yangtze, Yellow, Indus, Brahmaputra and Syr Darya) are recharged by the abundant snow/glacier melt water in HMA and provide a major headwater for almost 2 billion people [12,61,62,63,64,65]. In addition, the HMA has typical alpine vegetation cover, mainly including alpine grassland (1.52 × 10⁶ km²), alpine meadow, temperate grassland, deciduous broad-leaved forest and forest [66]. The changes in snow cover triggered by climate change have a significant impact on terrestrial ecosystems, notably affecting the phenology of alpine vegetation [67,68].

2.2. Data

2.2.1. MODIS Snow Cover Products

The version six MODIS Terra (MOD10A1 C6) and MODIS Aqua (MYD10A1 C6) daily snow cover products from 2000 to 2021 from the National Aeronautics and Space Administration (NASA) Earthdata website (https://search.earthdata.nasa.gov/search, accessed on 10 April 2022) are used in this study, covering the whole study area with nine tiles (Figure 1). The C6 product has implemented several algorithm improvements of snow cover and replaces the BSC and FSC layers in the C5 products with two new layers: raw NDSI and NDSI snow cover. The original MODIS NDSI layer is mosaicked, resampled and converted to a geographic projection (WGS84 coordinate system, 0.005 degree resolution). After that, the pixel value in the original MYD10A1 and MOD10A1 products are reclassified into two categories: (1) valid NDSI_Snow_Cover values ranging from 0 to 100 (with values for inland water and ocean set to 0); (2) invalid cloud-gap values of 250, encompassing the original classifications for cloud, missing data, detector saturated, night and no decision.

2.2.2. Ground Snow Depth (SD) Measurements

The SD observations at 99 meteorological stations (Figure 1) from 2001 to 2013 (https://data.tpdc.ac.cn/en/data/72d6dadf-8e1c-458b-b24e-91539042dfe6/, accessed on 5 August 2022) are utilized for accuracy validation [69,70]. The released SD data have performed quality control and homogenization procedures to be a high-quality data resource.

2.2.3. Landsat OLI Satellite Images

The Landsat-8 OLI satellite images (45 scenes, Figure 1) are used to derive high-resolution BSC maps for comparison and validation of the CGF MODIS NDSI product in this study. The images with less than 30% cloud coverage are selected through the ‘CLOUD_COVER’ under clear sky and visual identification. The images are provided and processed by the Google Earth Engine platform (https://earthengine.google.com/, accessed on 12 August 2022). The detailed selected images information is shown in

Table A2. Statistical description of Landsat-8 OLI images used for comparison and validation of the CGF MODIS NDSI product in this study.

Image Pair No.	Path/Row	Date of Acquisition	Cloud Cover (%)
1	134/37	2017/04/18	4.28
2	134/40	2017/06/05	11.78
3	134/40	2020/02/06	4.82
4	134/40	2021/03/28	4.44
5	135/33	2017/11/19	2.72
6	135/33	2019/05/01	3.12
7	135/33	2021/12/16	2.61
8	140/36	2018/04/15	11.66
9	141/36	2021/11/08	2.92
10	141/40	2016/12/28	2.78
11	141/40	2018/05/08	13.14
12	141/40	2019/02/04	3.49
13	141/40	2021/12/10	2.95
14	142/37	2016/01/18	7.05
15	142/37	2021/12/17	3.23
16	143/36	2017/06/04	12.29
17	143/36	2021/11/22	0.29
18	144/36	2021/11/13	2.04
19	145/30	2016/03/27	26.43
20	145/30	2018/10/27	4.9
21	145/30	2019/09/12	6.61
22	145/39	2017/06/02	12.39
23	145/39	2018/04/18	6.15
24	145/39	2019/12/01	3.27
25	145/39	2020/02/03	2.31
26	145/39	2021/04/26	2.75
27	146/35	2016/07/24	3.83
28	146/35	2018/05/11	3.37
29	147/37	2017/01/23	3.26
30	147/37	2017/07/02	13.45
31	147/37	2018/06/03	2.45
32	147/37	2019/10/28	1.49
33	148/35	2016/12/29	28.13
34	148/35	2018/05/09	12.37
35	148/35	2019/06/29	13.19
36	149/34	2018/01/24	2.99
37	149/34	2021/07/11	2.26
38	150/34	2016/07/20	2.46
39	150/34	2016/11/09	4.46
40	151/33	2017/04/09	0.3
41	151/33	2019/05/01	4.46
42	151/33	2020/08/23	2.37
43	151/33	2021/07/09	1.59
44	151/35	2017/04/09	0.43
45	151/35	2019/05/01	9.88

. All data processing is on the platform.

2.2.4. Digital Elevation Model (DEM) Data

The Shuttle Radar Topography Mission (SRTM) 90 m gridded digital elevation model (DEM) version 4.1 data are available from CGIAR Consortium for Spatial Information (CGIAR-CSI, http://srtm.csi.cgiar.org/, accessed on 13 April 2022). These data are resampled to 500 m for matching the MODIS snow cover product and used for extracting elevation information during the procedure of CGF method.

3. Methods

In this study, the implementation of the CGF method included three substeps: (1) the Terra and Aqua daily combination is executed to remove partial cloud gaps; (2) three temporal interpolation methods (TI) and a spatio-temporal weighted (STW) method are experimented for reconstructing NDSI information of cloud pixel separately; (3) the “Cloud Assumption” approach and cloud persistence days (CPD) are introduced to assess the performance of these four methods, and then the final CGF method is determined by how these four methods are combined (Section 4.1). Figure 2 presents a detailed flowchart of the proposed CGF method for MODIS NDSI product. Finally, daily CGF MODIS NDSI product is generated, and its accuracy is evaluated based on in situ SD observations and high-resolution snow cover maps derived from Landsat images.

3.1. Cloud-Gap-Filled (CGF) Method

3.1.1. Terra and Aqua Daily Combination (TAC)

The TAC refers to merge MOD10A1 (Terra) and MYD10A1 (Aqua) products on the same day on a pixel-by-pixel basis, efficiently reducing the cloud gaps with negligible precision sacrifice [71]. Thus, this method is introduced as a preprocessing to reduce cloud coverage preliminarily. The complete combination rule is: (1) when a pixel is cloud-free in both MOD10A1 and MYD10A1, the combination pixel value is set to the MOD10A1 NDSI value; (2) when a pixel is cloud-free only in one of the products, the combination pixel value is set to the cloud-free NDSI value; (3) for the remaining pixels, they are identified as cloud pixels; (4) for the first 2 years only Terra was in orbit, thus the period before 4 July 2002 is only based on MOD10A1 product.

3.1.2. Temporal Interpolation (TI) Method

After TAC, there are still a large number of clouds in the MODIS NDSI products. The NDSI can offer an accurate depiction of the gradual changes in snow cover over time, due to NDSI’s ability to accurately depict the gradual changes in snow cover over time, the clouds in MODIS NDSI data can be removed through temporal interpolation filtering. Tang et al. [71] developed a cubic spline interpolation (CSI) for MODIS NDSI products, and this method effectively filled in the missing data under the cloud cover. In this study, three temporal interpolation methods: linear interpolation (LI), quadratic interpolation (QI) and CSI are used for experimental test of temporal interpolation filtering. All interpolation processes are carried out in the Interactive Data Language platform (IDL, version 8.5).

3.1.3. Spatio-Temporal Weighted (STW) Method

However, the temporal interpolation algorithm may cause some error in the case of long cloud persistence days [71]. Inspired by the work of Jing et al. [24] and Li et al. [54], this study proposes STW as the comparison with the above three temporal interpolation methods, which simultaneously takes advantage of the spatial and temporal correlations of snow cover. In the proposed approach, the NDSI information of candidate cloud-free pixels within a certain spatio-temporal window are used to reconstruct the missing NDSI information of the target cloud pixel. The implementation of STW included three substeps: (1) the spatio-temporal selection of candidate cloud-free pixels; (2) the calculation of spatio-temporal weights for candidate cloud-free pixels; (3) the calculation of NDSI values for target cloud pixel based on the spatio-temporal weights (Figure 2).

In a target TAC image, N candidate cloud-free pixels are extracted within a 3 × 3 × t (7 ≤ t ≤ 15) space-time cube, which denotes a spatial window of dimensions 3 × 3 and a time window of t days with the center of the time window corresponding to the day of data gap. A space-time cube that is too small may not provide sufficient candidate cloud-free pixels, while a space-time cube that is too large can result in increased computational and time costs during the cloud removal process [72]. Thus, the value of t depends on whether it is satisfied such that N accounts for at least 30% (this threshold has been tested to be appropriate) of the total number of pixels in the space-time cube. In addition, many studies have demonstrated that snow fraction exhibits a positive linear correlation with elevation [73,74,75]. Thus, an elevation control condition (i.e., the absolute value of the difference in elevation between the candidate cloud-free pixels and the target cloud pixel should be less than or equal to 500) is added in this study to constrain the selection of candidate cloud-free pixels.

The reconstructed NDSI information of the target cloud pixel can be expressed by Equations (1)–(5). The weight W_i decides the contribution of the candidate cloud-free pixels to reconstructing NDSI information of the target cloud pixel and is calculated by the inverse distance weighted (IDW) method. It is determined by the temporal distance, geographic distance and elevational distance between the target cloud pixel and the candidate cloud-free pixels. Shorter time intervals, smaller distance and closer elevation of the candidate cloud-free pixel to the target cloud pixel produce a higher weight (i.e., greater contribution) for the candidate cloud-free pixel.

$$NDS{I}_k^T={\displaystyle {\sum }_i^N{W}_i}\times NDS{I}_i^C\cdots \begin{array}{cc}\mathrm{where}& W_i=(1/D_i)/{\displaystyle {\sum }_i^N(1/D_i)}\end{array}$$

(1)

$$D_i=\sqrt{{(d_i^t)}^2+{(d_i^g)}^2+{(d_i^e)}^2}$$

(2)

where ${NDSI}_k^T$ denotes the final weighted NDSI of the target cloud pixel k. ${NDSI}_i^C$ is the preprocessed (i.e., TAC) NDSI of the candidate cloud-free pixel i. N is the number of the candidate cloud-free pixels. W_i denotes the spatio-temporal weight. D_i denotes the synthetic spatio-temporal distance.

The temporal distance $d_i^t$ between the target cloud pixel k and the candidate cloud-free pixel i can be calculated according to Equation (3). The unit of time is the day.

$$d_i^t=1+abs(tim{e}_i-tim{e}_k)/tw$$

(3)

where tw is the width of searching temporal window that is used to normalize the temporal distance.

The geographic distance $d_i^g$ between the target cloud pixel ($x_k$, $y_k$) and the candidate cloud-free pixel ($x_i$, $y_j$) can be calculated according to Equation (4). The unit of horizontal or vertical distance is the number of pixels between the two pixels in horizontal or vertical direction.

$$d_i^g=1+\sqrt{{(x_i-x_k)}^2+{(y_i-y_k)}^2}/gw$$

(4)

where gw is the width of searching geographic window (gw = 1) that is used to normalize the geographic distance.

The elevational distance $d_i^e$ between the target cloud pixel k and the candidate cloud-free pixel i can be calculated according to Equation (5). The unit of elevation is meter.

$$d_i^e=1+abs(elevatio{n}_i-elevatio{n}_k)/ew$$

(5)

where ew is the width of searching elevational window (ew = 500 m) that is used to normalize the elevational distance.

3.2. Validation Method

3.2.1. Accuracy Assessment Based on “Cloud Assumption”

Utilizing the original NDSI value of the MODIS snow cover product as the validation data provides the most direct way to assess the accuracy of different TI methods and STW method in filling in cloud gaps over the entire study area. Therefore, the strategy we perform is to assume that the images on multiple random dates are cloud-covered. In addition, then we use different TI methods and STW method to reconstruct NDSI values (i.e., simulated NDSI values) for these images and subsequently compare them with the original images assumed before (i.e., reference true NDSI value). The TAC NDSI image of twenty dates (i.e., 12 February, 12 May, 13 October and 27 October of 2017; 5 February, 7 April, 31 July and 4 November of 2018; 16 February, 12 April, 12 August and 8 December of 2019; 27 March, 23 May, 17 September and 10 November of 2020; 20 March, 25 June, 26 August, 11 November of 2021) from 2017 to 2021 are selected. The performance of each method is evaluated using two metrics: mean absolute error (MAE) and root mean square error (RMSE), calculated using Equation (6) and Equation (7), respectively.

$$\mathrm{MAE}=\frac{{\displaystyle \sum _{i=1}^{\mathrm{n}}\left|c_i-s_i\right|}}{n}$$

(6)

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^{\mathrm{n}}{\left(c_i-s_i\right)}^2}}{n}}$$

(7)

where c_i and s_i are the simulated NDSI value and reference true NDSI value, respectively; n is the total number of reference true NDSI pixels.

In addition, we introduce the cloud persistence days (CPD) combined with MAE and RMSE to evaluate the accuracy of different reconstruction methods and determine the final CGF method. The CPD represents the number of consecutive days of cloud observed for a pixel from the last cloud-free observation to the next cloud-free observation [71]. The CPD for each cloud pixel are calculated using the daily TAC images from 2017 to 2021.

3.2.2. Validation Based on In-Situ SD Observations

The CGF MODIS NDSI product are compared with the daily SD observed by ninety-nine meteorological stations from 2001 to2013. The confusion matrix widely used [76], presented in

Table 1. Accuracy evaluation of confusion matrix based on in situ SD observations. The metrics derived from confusion matrix are: overall accuracy (OA), underestimation error (MU) and overestimation error (MO). The ε₁ and ε₂ indicate the defined SD threshold value and MODIS NDSI threshold value, respectively.

		Observed SD
		Snow Cover (≥ε1 cm)	No Snow (<ε1 cm)
MODIS NDSI	Snow cover (≥ε2)	a	b
	No Snow (<ε2)	c	d
	Cloud	e	f
OA = ((a + d)/(a + b + c + d + e + f)) × 100%
MU = (c/(a + b + c + d)) × 100%   MO = (b/(a + b + c + d)) × 100% Number of available pixels = a + b + c + d

. The ε₁= 1, 2, 3 and 5 cm are the defined threshold values (referring to different measurements of SD) for the in situ SD observations. In addition, the ε₂= 0.10 [57], 0.29 [58] and 0.40 [37] are the defined threshold values for the MODIS NDSI to determine whether a pixel is covered by snow.

3.2.3. Validation Based on Landsat-8 OLI Images Derived BSC Maps

The BSC maps with 30 m resolution retrieved from Landsat-8 OLI images based on the Google Earth Engine platform are used to verify the CGF MODIS NDSI product in this study. The method of BSC map generation consists of two parts: the first part is the SNOMAP algorithm [2,77], which defined the NDSI of a snow pixel as greater than or equal to 0.40. Furthermore, to eliminate the influence of water bodies and other land cover, the reflectance of NIR greater than 0.11 and green band greater than 0.10 are adopted in this algorithm [78]. The second part is proposed by Wang et al. [79] and the normalized difference vegetation index (NDVI), which are introduced in this study to enhance the identification of the forest snow cover. The specific strategy is that a pixel with NDSI < 0.4, NDFSI > 0.4 and NDVI < 0.6 is identified as the snow pixel in the forest [80]. The NDSI, NDFSI and NDVI are defined as follows Equations (8)–(10):

$$\mathrm{NDSI}=\frac{{\rho }_{\mathrm{green}}-{\rho }_{\mathrm{swir}}}{{\rho }_{\mathrm{green}}+{\rho }_{\mathrm{swir}}}$$

(8)

$$\mathrm{NDFSI}=\frac{{\rho }_{\mathrm{nir}}-{\rho }_{\mathrm{swir}}}{{\rho }_{\mathrm{nir}}+{\rho }_{\mathrm{swir}}}$$

(9)

$$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{nir}}-{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{nir}}+{\rho }_{\mathrm{red}}}$$

(10)

where ρ_green, ρ_swir, ρ_nir and ρ_red correspond to the reflectance of the green, SWIR, NIR and red bands measured in Landsat imagery.

In addition, the 500 m BSC maps are aggregated from the Landsat 30 m BSC maps through 50%-pixel aggregation to match the spatial resolution of the MODIS NDSI product. The specific strategy is that if an aggregated pixel (500 m) contains more than 50% snow cover pixels (30 m) and lower than 5% cloud pixels (30 m), it is assigned to “snow covered”; otherwise, it is assigned “snow free” or “cloud” (those with >5% of the 30 m pixels flagged as “cloud” or “cloud shadow”) [58,78]. Both cloud and cloud shadow are identified using the quality assessment (QA) band of Landsat-8 OLI image. Then, the CGF MODIS NDSI derived BSC maps are evaluated with the Landsat BSC maps based on the confusion matrix including the metrics of OA, commission error (CE) and omission error (OE) (

Table 2. Confusion matrix comparing MODIS NDSI derived BSC maps with Landsat BSC maps. The ε₂ indicates the defined MODIS NDSI threshold value.

		Landsat OLI
		Snow Cover	No Snow
MODIS NDSI	Snow cover (≥ε2)	SS	NS
	No now (<ε2)	SN	NN
OA = ((SS + NN)/(SS + NS + SN + NN)) × 100%
OE = (SN/(SS + SN)) × 100%   CE = (NS/(NN + NS)) × 100%

).

4. Results

4.1. The Determination and Effectiveness of the CGF Method

As shown in Figure 3a, the frequency of the CPD reduces gradually as the increasing of CPD value, and the frequency of CPD ≤ 5 d and CPD ≤ 15 d for the cloud-covered pixels reach 70.92% and 94.07%, respectively. It can be clearly seen that the accuracy of these four methods is related to CPD with the MAE and RMSE increasing as the rising of CPD (Figure 3b,c). The accuracy of LI and QI is obviously lower than the CSI and STW due to higher values of MAE and RMSE. For the CSI method and STW method, when the CPD < 8 d, the MAE and RMSE curves of the two methods are approximately the same, i.e., a close accuracy; in the case of the CPD ≥ 8 d, the STW has the slightly higher accuracy than CSI due to the lower values of MAE and RMSE. The overall MAE/RMSE of the methods of LI, QI, CSI and STW using the frequency of CPD as the weight are 0.080/0.137, 0.088/0.144, 0.066/0.101 and 0.064/0.099, respectively.

To directly show the cloud removal results, Figure 4 shows the spatial comparison in three regions (R1, R2 and R3, Figure 1), including the TAC NDSI image, real cloud and cloud assumption image, CPD image, CSI NDSI image, STW NDSI image and the CSI-STW NDSI image. The CSI-STW image is the combination of the CSI NDSI image (when the CPD < 8 d, and the frequency of CPD is 80.53%) and STW NDSI image (when the CPD ≥ 8 d, and the frequency of CPD is 19.47%) based on the CPD NDSI image. It can be clearly seen that both CSI NDSI and STW NDSI images can generally maintain well the spatial distribution of snow cover compared with the true TAC NDSI image, and the distribution of reconstructed pixels have spatial continuity. However, there are some abrupt changes in the reconstructed pixels of the CSI NDSI image with high CPD. In addition, the reconstructed result of STW NDSI image tends to have some overestimations in low CPD region, but performs better in high CPD region than the CSI NDSI image. Therefore, the combination of the CSI and STW (i.e., CSI-STW NDSI image) is determined as the final CGF method, which can take full advantage of the respective precision advantages and spatial accuracies of CSI and STW NDSI (Figure 4, column 6). Moreover, this combination ensures high efficiency, due to higher time cost (high computational complexity) for STW. The CPD threshold (CPD = 8 d) for the combination is taken from the accuracy results in Figure 3.

4.2. The Accuracy of the CGF MODIS NDSI Product

4.2.1. Validation Based on In Situ SD Observations

The in situ SD observations are used as the ground truth to verify the snow classification accuracy of the TAC NDSI and CSI-STW NDSI datasets (

Table 3. Validation results for TAC NDSI and CSI-STW NDSI based on in situ SD observations in the Tibetan Plateau from 1 January 2001 to 31 December 2013. The definitions and calculations of the indices in this table are shown in Table 1.

Threshold Values for Snow Cover	TAC NDSI			CSI-STW NDSI
	OA (%)	MU (%)	MO (%)	OA (%)	MU (%)	MO (%)
ε1 = 1 cm, ε2 = 0.1	59.21	1.70	3.23	93.54	2.32	4.14
ε1 = 1 cm, ε2 = 0.29	60.26	2.07	1.16	95.39	3.15	1.46
ε1 = 1 cm, ε2 = 0.4	60.43	2.26	0.70	95.67	3.49	0.84
ε1 = 2 cm, ε2 = 0.1	59.61	0.71	3.57	94.16	1.19	4.65
ε1 = 2 cm, ε2 = 0.29	60.86	0.93	1.34	96.71	1.56	1.73
ε1 = 2 cm, ε2 = 0.4	61.10	1.07	0.83	97.18	1.81	1.01
ε1 = 3 cm, ε2 = 0.1	59.66	0.32	3.87	94.30	0.61	5.10
ε1 = 3 cm, ε2 = 0.29	61.03	0.45	1.55	97.15	0.83	2.02
ε1 = 3 cm, ε2 = 0.4	61.33	0.54	0.99	97.78	0.99	1.22
ε1 = 5 cm, ε2 = 0.1	59.52	0.08	4.35	94.08	0.15	5.77
ε1 = 5 cm, ε2 = 0.29	61.00	0.11	1.94	97.26	0.20	2.54
ε1 = 5 cm, ε2 = 0.4	61.36	0.15	1.33	98.08	0.27	1.65
Average value	60.45	0.87	2.07	95.94	1.39	2.68
Number of available pixels	289,595			465,010

). For different threshold combinations of ε₁ and ε₂, the OAs, MUs and MOs obtained by the CSI-STW NDSI compared with the TAC NDSI show satisfactory performance. The CSI-STW NDSI efficiently captures the snow cover referring to the in situ measurements, with an average OA of more than 95%. However, the TAC NDSI is insufficient to accurately detect the snow cover, with an average OA of only 60.45%. Note: The reason for the significant difference in OA between TAC NDSI and CSI-STW NDSI is that the calculation of OA is considering the case where TAC is cloud-covered when station records SD (e and f in Table 1). Although the MU and MO of the CSI-STW NDSI are slightly higher than the TAC NDSI, they are still within an acceptable range. To sum up, the CGF method proposed in this study can fill all of the data-gap pixels, capture more snow events and restore the snow cover information with a high reliability.

4.2.2. Validation Based on High-Resolution BSC Maps

In the absence of meteorological station data, a good way to evaluate the accuracy of the CGF MODIS NDSI product is to compare it with snow cover maps derived from higher-resolution sensors. In this study, a total of 45 Landsat-8 images (Figure 1, Table A2) with different snow coverages from January 2017 to December 2021 are selected for this accuracy validation. As shown in

Table 4. Validation results for TAC NDSI and CSI-STW NDSI based on Landsat-8 BSC maps under clear-sky conditions. The definitions and calculations of the indices in this table are shown in Table 2.

Metrics (%)	NDSI Threshold: 0.1		NDSI Threshold: 0.29		NDSI Threshold: 0.40
	TAC NDSI	CSI-STW NDSI	TAC NDSI	CSI-STW NDSI	TAC NDSI	CSI-STW NDSI
SS	1,765,524	2,255,700	1,688,278	2,158,014	1,608,492	2,053,084
NS	575,201	662,668	290,431	340,524	190,413	225,094
SN	21,234	32,577	98,480	130,263	178,266	235,193
NN	2,990,539	3,107,687	3,275,309	3,429,831	3,375,327	3,545,261
Total	5,352,498	6,058,632	5,352,498	6,058,632	5,352,498	6,058,632
OE	1.19	1.42	5.51	5.69	9.98	10.28
CE	16.13	17.58	8.15	9.03	5.34	5.97
OA	88.86	88.52	92.73	92.23	93.11	92.40

, the OA of TAC NDSI is 88.86–93.11% for different MODIS NDSI thresholds. Based on the CGF method, the data gaps are filled very effectively, with the OA of CSI-STW NDSI (88.52–92.40%) slightly lower than TAC NDSI (with less loss of accuracy).

Table 5. Validation results of CSI-STW NDSI for filling the data-gap pixels in the TAC NDSI dataset, based on the Landsat-8 BSC maps under clear-sky conditions. The definitions and calculations of the indices in this table are shown in Table 2.

Metrics (%)	NDSI Threshold: 0.1	NDSI Threshold: 0.29	NDSI Threshold: 0.40
SS	489,451	469,736	444,592
NS	85,163	50,093	34,681
SN	12,068	31,783	56,927
NN	119,452	154,522	169,934
Total	706,134	706,134	706,134
OE	2.41	6.34	11.35
CE	41.62	24.48	16.95
OA	86.23	88.41	87.03

further presents the validation results of CSI-STW NDSI for filling the data-gap pixels in the TAC NDSI dataset, based on the Landsat-8 BSC maps under clear-sky conditions. The OA of CSI-STW NDSI for data-gap pixels in TAC NDSI is 86.23–88.41% for different MODIS NDSI thresholds, which is 2.63–6.08% slightly lower than TAC NDSI (gap-free pixels) in Table 4. The results from Table 4 and Table 5 demonstrate that the CGF method proposed in this study can develop more reliable snow cover products, with high OA and acceptable OE and CE. In addition, the MODIS NDSI threshold is an important factor affecting the accuracy of MODIS snow cover mapping. The OA tends to increase, the CE tends to decrease and the OE tends to increase as NDSI threshold increases. However, the NDSI threshold of 0.29 notably improves the detection of snow cover and achieve outstanding OE and CE, compared to the other two thresholds. It can be concluded that the NDSI threshold of 0.29 is the best threshold to describe the spatial pattern of snow cover of HMA when using MODIS snow-cover data, which is coherent with the conclusion of Zhang et al. [58]. Figure 5 shows the snow cover mapping results, i.e., TAC BSC, Landsat BSC, Landsat BSC aggregation and CSI-STW BSC, correspond to four Landsat scenes to facilitate the analysis of their spatial consistency. In addition, the MODIS NDSI threshold is set to 0.29. It can be seen that the spatial distribution of CSI-STW BSC is fairly consistent with Landsat BSC, with notably improvement in the detection of snow cover information compared to the TAC BSC. Overall, the CGF MODIS NDSI product can accurately reflect the spatio-temporal pattern of snow cover in HMA and reliably serve the research field of hydrology and climate change.

4.3. Spatio-Temporal Patterns of Snow Cover over HMA

The CGF MODIS NDSI product developed in this study can be utilized to detect the spatio-temporal patterns of snow cover over HMA. The snow cover over HMA exhibit large spatial heterogeneity due to the complex topographic and climatic conditions, and a characteristic of a marked reduction in space during the snow melt period (Figure 6). Using 0.29 of the MODIS NDSI threshold, the snow-covered days (SCD) and snow-covered extent (SCE) can be calculated [71]. As seen in Figure 7, the areas with high SCD (SCD > 120 d) are mainly distributed in high-altitude mountain ranges, such as Tien Shan, Karakoram, Kunlun, Qilian Shan, Hengduan Shan and Himalayas, accounting for 22.25% of HMA (elevation greater than 1500 m). The areas with low SCD (SCD < 20 d) are primarily distributed in the interior of the Tibetan Plateau and low altitude areas around it, accounting for 34.45%. Figure 8 displays the interannual variations in the daily snow-covered extent (SCE), monthly SCE and yearly SCE for HMA (elevation greater than 1500 m) from 2000 to 2021. The daily SCE and monthly SCE presents strong interannual volatility in winter months, while there is steady fluctuation in summer months. As for the yearly SCE, the mean SCE and minimum SCE show smooth interannual fluctuations with an insignificant trend. However, the maximum SCE shows significant interannual fluctuation, which can reflect some extreme weather events to some degree, such as the extreme snow event in 2008.

5. Discussion

In recent years, several researchers have successively developed good long-term cloud-free snow cover products for different regional scales (

Table 6. Typical long-term cloud-free snow cover products for different regional scales.

Product Type	References and Dataset DOI	Spatial Coverage	Temporal Coverage	Temporal Resolution	Spatial Resolution
SCE	Huang [83] https://doi.org/10.12072/ncdc.CCI.db0044.2020, accessed on 21 March 2022	Northern Hemisphere	2000–2015	Daily	∼1 km
SCE	Hao et al. [84] https://doi.org/10.11888/Snow.tpdc.271381, accessed on 21 March 2022	China	1981–2019	Daily	∼5 km
SCE	Hao et al. [85] https://doi.org/10.12072/ncdc.I-SNOW.db0001.2020, accessed on 15 May 2022	China	2000–2020	Daily	∼500 m
SCE	Muhammad and Thapa [86] https://doi.org/10.1594/PANGAEA.918198, accessed on 21 March 2022	HMA	2002–2019	Daily	∼500 m
SCE	Huang et al. [87] https://doi.org/10.11888/Cryos.tpdc.272204, accessed on 30 September 2022	Tibetan Plateau	2002–2021	Daily	∼500 m
SCE	Li et al. [88] https://doi.org/10.57760/sciencedb.j00076.00112, accessed on 27 September 2022	HMA	1982–2019	Daily	∼5 km
FSC	Qiu et al. [89](<10% Cloud coverage) https://doi.org/10.11922/sciencedb.457, accessed on 21 March 2022	HMA	2002–2018	Daily	∼500 m
NDSI	Hall and Riggs [81] https://doi.org/10.5067/MODIS/MOD10A1F.061, accessed on 21 March 2022	Global coverage	2000–present	Daily	∼500 m
NDSI	Han et al. [90] https://doi.org/10.12072/ncdc.I-SNOW.db0024.2021, accessed on 21 March 2022	Northeast China	2000–2020 (snow season)	Daily	∼500 m
NDSI	Tang et al. [71] https://doi.org/10.11888/Cryos.tpdc.272836, accessed on 29 September 2022	HMA	2000–2021	Daily	∼500 m
NDSI	Jing et al. [59] https://doi.org/10.5281/zenodo.5644386, accessed on 12 July 2022	China	2001–2020	Daily	∼500 m

). Overall, cloud-free SCE products produced by composite algorithms are frequently released, while only a few NDSI products have been produced. The global cloud-gap-filled MODIS NDSI dataset (MOD10A1F) were generated by retaining clear-sky views of the surface from previous days in MOD10A1 NDSI product to fill the cloud-covered pixel [81,82]. However, this product performs poorly in China, where periodic and transient snow is dominant [59], and so does HMA. Jing et al. [59] developed the Spatio-Temporal Adaptive fusion method with erroR correction to generate the cloud-free STAR NDSI collection. This product provides a detailed snow cover dataset with high accuracy for China and have excellent application prospects. Tang et al. [71] developed a cloud removal method based on CSI for MODIS NDSI products. Even though this method effectively filled in the missing data under the cloud cover, it may cause some error in the case of long cloud persistence days. Our study is an inheritance and a great improvement of our previous work [71] in a sense. In this study, three temporal interpolation methods (LI, QI and CSI) and a STW method are synthesized and compared for determining the final CGF method. Accuracy evaluations based on in situ SD observations and high-resolution Landsat-8 OLI images derived snow cover maps verified the reliability of the CGF MODIS NDSI product in terms of accuracy and consistency. In addition, the CGF MODIS NDSI product effectively improves the accuracy of NDSI reconstructions in areas with long cloud persistence days. Additionally based on this high accuracy product, the binary snow cover (BSC) product could be derived using a certain MODIS NDSI threshold (such as 0.29 of Zhang et al. [58]) or the fractional snow cover (FSC) product could be derived using the empirical relationship of Zhang et al. [59] between FSC and MODIS NDSI, which in turn can be engaged in other applications.

However, limitations of the CGF MODIS NDSI product (including its accuracy validation) may come from the following: (1) the original MODIS NDSI products are not perfect due to the difficulty of detecting snow in the mountainous area with complex terrain and land cover, as well as snow/cloud confusion errors in the cloud-masking algorithm for MODIS snow data; (2) the slope and aspect of the terrain, and forest canopy obstruction may affect the reconstruction of NDSI or the accuracy of NDSI from MODIS original products [33,34]; (3) the NDSI approach cannot always accurately distinguish between snow and cloud, which leads to snow misclassification errors [91,92]; (4) the in situ SD observations are sorely scarce, especially in high-altitude areas, which creates challenges in verifying snow category classification [58]; when there is high cloud coverage in Landsat images, aggregating from the Landsat snow cover map with 30 m spatial resolution to the snow cover map with 500 m MODIS spatial resolution may introduce some uncertainties; and (5) there are also spatial inconsistencies between the CGF MODIS NDSI product and Landsat snow cover, possible reasons are as follows: firstly there is a difference between the Landsat OLI and MODIS sensors and their snow cover mapping algorithms; secondly the snow cover information obtained from the two images is inconsistent, even on the same day, because the snow can change rapidly affected by wind speed, sunshine duration, air temperature and precipitation [78]; lastly due to scaling effects, Landsat is superior to MODIS in portraying snow cover in areas of complex terrain. In addition, as the temporal and spatial resolution of satellite remote sensing increases, different strategies need to be adopted when reconstructing snow cover information, taking into account different snow conditions, microtopographic factors (e.g., as slope gradient and aspect), land-cover types and vegetation coverage in forests.

6. Conclusions

The MODIS NDSI snow cover product is susceptible to the influence of the cloud contamination, leading to numerous data gaps. In this study, an effective CGF method was developed to fully fill the data gaps in MODIS NDSI snow cover product. The CGF MODIS NDSI product generally has the following strengths:

(1) The proposed CGF method combines the respective strengths of the CSI method (high accuracy and computational efficiency in the case of short cloud persistence days) and the STW method (comprehensively considering spatial and temporal correlations of the snow cover).

(2) An accuracy evaluation based on in situ SD observations verified the reliability of the CGF MODIS NDSI product in terms of accuracy and consistency. The CGF MODIS NDSI product achieves a high-range OA of 93.54–98.08%, a low-range MU error of 0.15–3.49% and an acceptable-range MO error of 0.84–5.77% for different combinations of NDSI and SD observed thresholds. Compared with the high-resolution Landsat-8 OLI images derived snow cover maps, the CGF MODIS NDSI product largely corresponds to the Landsat snow cover maps, even if there are scaling effects. In addition, the CGF MODIS NDSI product effectively improves the accuracy of NDSI reconstructions in areas with long cloud persistence days.

(3) Overall, the CGF MODIS NDSI product performs well and is able to provide a set of long-term, spatiotemporally continuous and highly accurate snow cover dataset for HMA, and thereby provide a valuable input dataset for hydrological and climate modeling, snow cover dynamics and other water-related studies.