A Conditional Probability Interpolation Method Based on a Space-Time Cube for MODIS Snow Cover Products Gap Filling

Abstract

Seasonal snow cover is closely related to regional climate and hydrological processes. In this study, Moderate Resolution Imaging Spectroradiometer (MODIS) daily snow cover products from 2001 to 2018 were applied to analyze the snow cover variation in northern Xinjiang, China. As cloud obscuration causes significant spatiotemporal discontinuities in the binary snow cover extent (SCE), we propose a conditional probability interpolation method based on a space-time cube (STCPI) to remove clouds completely after combining Terra and Aqua data. First, the conditional probability that the central pixel and every neighboring pixel in a space-time cube of 5 × 5 × 5 with the same snow condition is counted. Then the snow probability of the cloud pixels reclassified as snow is calculated based on the space-time cube. Finally, the snow condition of the cloud pixels can be recovered by snow probability. The validation experiments with the cloud assumption indicate that STCPI can remove clouds completely and achieve an overall accuracy of 97.44% under different cloud fractions. The generated daily cloud-free MODIS SCE products have a high agreement with the Landsat–8 OLI image, for which the overall accuracy is 90.34%. The snow cover variation in northern Xinjiang, China, from 2001 to 2018 was investigated based on the snow cover area (SCA) and snow cover days (SCD). The results show that the interannual change of SCA gradually decreases as the elevation increases, and the SCD and elevation have a positive correlation. Furthermore, the interannual SCD variation shows that the area of increase is higher than that of decrease during the 18 years.

1. Introduction

Snow cover has a high albedo, high emissivity, and an important impact on regional climate and energy variation [1,2,3]. Additionally, snow is considered to be one of the sensitive factors to reflect climate change, especially regarding precipitation and temperature [4,5,6,7]. To understand climate change interactions and predicting climate warming tendencies, it is imperative to monitor accurately the snow cover extent (SCE) [8]. Moreover, snow cover is also a critical water source and significantly influences the hydrological and biochemical processes in surrounding lowland areas [9,10,11]. Hence, investigation of the spatiotemporal characteristics of snow is meaningful for effective management of limited water resources.

Snow is a main water resource in northern Xinjiang, China, and there is long-term snow cover in winter [4]. Sufficient snowfall is beneficial for agricultural irrigation in spring and crop wintering, but it also brings about frequent snow disasters [12]. Due to the influence of topography and geomorphology, the distribution of snow cover is uneven in different regions and seasons. Therefore, it is necessary to generate cloud-free SCE products and analyze the spatiotemporal change in snow cover in northern Xinjiang, China.

Remote sensing has shown great advantages in monitoring snow cover in large-scale and rugged terrain regions. Depending on the data source, research on snow variations includes in situ-based analysis, microwave-based analysis, and optical-based analysis [13,14]. The in situ data, which provide accurate and long record point measurements, are suitable for the analysis of climate tendency [15,16]. However, in situ snow monitoring is too rare to adequately quantify snow variations [17,18]. Passive microwave remote sensing data are effective in monitoring the spatiotemporal variations in snow cover by the snow depth and the snow water equivalent. However, their spatial resolution is usually low (>1 km), and detailed snow information is not acquired from them [19,20].

In optical remote sensing data, the Moderate Resolution Imaging Spectroradiometer (MODIS), which is onboard the Terra and Aqua satellites, is widely used to detect snow cover change because it has a high temporal resolution, high spatial resolution (500 m), and high accuracy of binary snow cover products [21,22,23,24,25]. Yang et al. [8] compared the accuracy of the MOD/MYD10A1 snow products and the Interactive Multisensor Snow and Ice Mapping System (IMS) products, and showed that the daily MODIS snow cover products have better performance. However, the daily MODIS snow cover products have many spatial and temporal gaps that are mainly due to cloud contamination.

To remove MODIS cloud pixels, several methods have been proposed in recent decades. MOD10A2 and MYD10A2 generated an 8-day maximum SCE, which reduced the temporal resolution of the original data [26]. Multisource fusion methods make full use of the high spatial resolution of MODIS and the good cloud-penetrating capacity of microwave products to remove cloud; e.g., MODIS and the advanced microwave scanning radiometer for earth observing system (25 km) [27,28,29] and MODIS and IMS (4 km) [30] have been fused. These methods resample the coarse spatial resolution of microwave products to 500 m to replace the cloud pixels in the MODIS products, and the fusion reduces the spatial resolution and leads to some uncertainty due to the reclassification accuracy of the cloud pixel depending on the microwave products [31,32,33]. Meanwhile, other spatial filter and temporal filter methods have also been developed to recover cloud gaps. Parajka and Blöschl [34] first presented a cloud removal method with three steps: Terra and Aqua image combination (TAC), spatial filtering using the nearest neighbors, and temporal filtering. This method has been improved by incorporating the snow line, air temperature, elevation information, and so on [35,36,37]. However, all these cloud removal methods use only spatial and temporal information step by step and ignore some important neighborhood information at small spatiotemporal distances.

Dozier et al. [38] first proposed the concept of taking the snow cover variation as a space-time cube to fill cloud gaps in fractional snow cover. Huang et al. [39] established a hidden Markov random field framework to fill SCE gaps, which were also based on a space-time cube and the snow mapping accuracy reached 88.0% for gap-filled areas. A space-time cube can make full use of neighborhood information to recover the ground characteristics. Additionally, the conditional probability interpolation method [40,41] can effectively calculate the conditional probability of snow pixels being covered by clouds and meteorological data to remove clouds, but this method has limited capacity for removing clouds in areas with few in situ observations. Therefore, we propose a conditional probability interpolation method based on a space-time cube (STCPI) that combines the advantages of a space-time cube and conditional probability interpolation method to produce daily cloud-free MODIS SCE products. This method takes the conditional probability as the weight of the space-time neighborhood pixels to calculate the snow probability of the cloud pixels, and then the snow condition of the cloud pixels can be recovered by the snow probability. The process can fill all of the cloud gaps in the MODIS SCE products.

The generated daily cloud-free MODIS snow products can effectively quantify the spatial and temporal changes of snow cover in northern Xinjiang, China. Yang et al. [42] analyzed the change of snow cover in northern Xinjiang with site data, showing that snow cover increased from 1961 to 2017. Chen et al. [43] used MODIS 8-day synthetic product to analyze the correlation between snow cover and elevation in Xinjiang, China, and there are some other studies about the snow cover variation in Tianshan Mountains [44,45,46]. However, systematic analysis of the spatial and temporal variation of snow cover in northern Xinjiang, China, is still rare, and current analyses are mainly based on the site data and MODIS version 5 snow products. In MODIS version 5 snow products, snow is detected from clouds and other land cover types with the normalized difference snow index (NDSI) threshold of 0.4 [47,48]. Compared with version 5, the new MODIS version 6 products released in 2016 better handles atmospheric correction, cloud and snow misclassification, recovers Aqua band 6, and provides NDSI products. Many researchers assessed the MODIS version 6 snow products exhibited that an NDSI threshold of 0.1 provides the highest accuracy [49,50]. Therefore, it is necessary and urgent to carry out more accurate spatiotemporal analysis of snow cover in northern Xinjiang, China based on higher precision version 6 snow products.

As a result, we first propose a new cloud remove method to generate daily cloud-free MODIS SCE products. The validation procedure of this method is based on the cloud assumption method and Landsat–8 OLI images. Then, we investigate the snow cover spatiotemporal variation by the snow cover area (SCA) and snow cover duration (SCD) in northern Xinjiang, China, from 2001 to 2018.

2. Study Area and Data

2.1. Study Area

Northern Xinjiang, China, is located on the Eurasian continent in the range of 42°–50°N latitude and 79°–92°E longitude. This region covers a total area of 0.39 million km² with an average elevation of 1256 m above sea level (Figure 1). The topographical features are the Junggar Basin situated in the middle of the region, and Altai Mountains and Tianshan Mountains are distributed on both sides of north and south. Due to a considerable amount of snowfall in this region, snow resources play a vital role in the water supply and ecosystem.

2.2. Data

Version 6 of the NDSI snow cover products in MOD/MYD10A1 was employed in this study. To conveniently fill in gaps, we reclassified the NDSI snow cover into three types–NDSI $\ge$ 0.1 was reclassified as snow; NDSI $<$ 0.1, inland water and ocean were reclassified as snow-free; the other classes (cloud, missing, no decision, night, detector saturated, and fill) were combined into a cloud class. The data were acquired from the Google Earth Engine (https://code.earthengine.google.com); there were 6556 images of MOD10A1 and 5959 images of MYD10A1 from 1 January 2001, to 31 December 2018.

Landsat–8 OLI data, which has a high spatial resolution of 30 m, were derived from the U.S. Geological Survey (USGS) (https://Earthexplorer.usgs.gov) to validate the produced daily cloud-free MODIS SCE products. Figure 1 shows the locations of the Landsat–8 OLI images.

Digital elevation data were derived from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) version 2 with a 30 m spatial resolution. To match the MODIS snow products, the GDEM was resampled to 500 m grids by calculating the average elevation. These data were mainly used to analyze snow cover variation in different altitudinal zones.

3. Methodology

3.1. Gap-Filling Method

We propose a new gap-filling method that aims to produce daily cloud-free MODIS SCE products. First, TAC is used to remove partial cloud pixels, and then using STCPI method removed the remaining cloud completely. The cloud removal procedure is shown in Figure 2.

3.1.1. Terra and Aqua Daily Combination (TAC)

TAC has been used in many kinds of research to remove clouds [34,35,36,37,40]. As the Aqua and Terra satellites have different transit times but similar designs and performance, it is feasible to combine Terra (MOD10A1) and Aqua (MYD10A1) data to remove clouds. The combination rule is that first, when a pixel is snow in both MOD10A1 and MYD10A1, snow is taken as daily combination pixels; after that, when a pixel is snow-free in only one of the products, snow-free is taken as daily combination pixels; finally, the remainder pixels are taken as cloud for daily combination pixels [36]. Through defining snow > snow-free > cloud, this rule can be expressed as:

$$S_{\left(x,y\right)}^C=max(S_{\left(x,y\right)}^T,S_{\left(x,y\right)}^A)$$

(1)

where x represents the row location, y represents the column location, and S^T, S^A, and S^C represent the Terra, Aqua, and daily combination pixels, respectively. In this work, we refer to the daily combination products as MOYD.

3.1.2. Conditional Probability Interpolation Based on a Space-Time Cube (STCPI)

After TAC, there are still a large number of clouds in the MODIS SCE products. To fill the cloud gaps and ensure accuracy, STCPI is applied to recover snow conditions from the surrounding cloud-free pixels. This method takes the conditional probability as the weight of the space-time neighborhood pixels to calculate the snow probability of the cloud pixels, and then the snow condition of the cloud pixels can be recovered.

First, the conditional probability of each pixel in a clear sky is calculated, given that each pixel is in a space-time cube. The conditional probabilities are calculated using the following formula:

$$p(C_{x,y}|C_n)=\frac{{\displaystyle \sum }\left(1-ABS\left(C_{x,y,t}-C_{n,t\prime }\right)\right)}{N_{x,y}}$$

(2)

where $p\left(C_{x,y}|C_n\right)$ is the conditional probability of the central pixel with coordinates x,y and the n-th neighboring pixel in the space-time cube having the same snow condition. $C_{x,y,t}$ and $C_{n,t\prime }$ represent snow (C = 1) or snow-free (C = 0) conditions in the MODIS SCE products for days t and t´, respectively. Nx,y is the number of days in which the central pixel and n-th neighboring pixel are cloud-free over the study period 2001–2018. $p\left(C_{x,y}|C_n\right)$ ranges from 0 to 1. The closer $p\left(C_{x,y}|C_n\right)$ is to 1, the more consistent the changes in the n-th neighboring pixel and central pixel. That is, the n-th neighboring pixel and central pixel have a higher probability of having the same snow condition. Conversely, the closer $p\left(C_{x,y}|C_n\right)$ is to 0, the greater the probability that the n-th neighboring pixel and central pixel have opposite snow conditions.

Figure 3 shows the conditional probability for the snow conditions of the central pixel with coordinates (x,y,t) and neighboring pixel with coordinates (x+2,y+2,t+2). The conditional probability in the mountains is lower than that on flat ground, especially in the Tianshan Mountains, because of the greater heterogeneity of the land surface in mountainous areas. To ensure that every central pixel has corresponding high conditional probability neighborhood pixels, and exclude the influence of low conditional probability neighborhood pixels to the final results, only the conditional probability $p\left(C_{x,y}|C_n\right)>0.8$ is used to calculate the snow probability of the central pixel in this study. In general, the smaller the space-time distances between pixels, the greater the probability that they have the same snow condition. Therefore, it is necessary to select the smallest space-time cube with sufficient cloud-free neighborhood pixels. A too small space-time cube does not ensure sufficient cloud-free neighboring pixels to calculate the snow probability of central pixel, and a too big space-time cube will lead to more computational and time cost in the process of cloud removal. To maintain a small computational cost and ensure sufficient cloud-free pixels in the neighborhood space-time cube, the space-time cube of a 5 × 5 × 5 window (a spatial window of dimensions 5 × 5 and a time window of 5 days) was appropriate in this experiment.

Based on the conditional probability calculation in the space-time cube, we can use the current snow condition of the space-time cube to calculate the final snow probability ($P\left(x0,y0,t0\right))$ of the cloud gaps:

$$P\left(x0,y0,t0\right)=\frac{{\displaystyle \sum }P(C_{x0,y0}|C_n)\times Sn }{{\displaystyle \sum }P(C_{x0,y0}|C_n) \times Vn},$$

(3)

where

$$\begin{gathered} S_n =\{\begin{array}{c}1, if n­th pixel has snow\hfill \\ 0, otherwise\hfill \end{array}, \\ V_n =\{\begin{array}{c}1, if n­th pixel is cloud­free\hfill \\ 0, otherwise\hfill \end{array} \end{gathered}$$

As the range of snow probability is [0,1], and the real category of cloud pixels is only snow or snow-free, reconstructing cloud pixel information is a binary classification problem of judging cloud pixels to be snow or snow-free, which is an equal possibility event. Therefore, taking the median value of 0.5 as the threshold to judge the category of cloud is appropriate, and the snow condition $C_{\left(x0,y0,t0\right)}$ can be determined by Equation (4). Then, the cloud pixels can be removed completely by performing STCPI iteration.

$$C_{\left(x0,y0,t0\right)} =\{\begin{array}{c}snow, P\ge 0.5\hfill \\ snow­free, P<0.5\end{array}$$

(4)

3.2. Validation Methodology

3.2.1. Validation Based on the Cloud Assumption

The truth snow condition in cloud gaps is unknown. Here, we apply a cloud assumption method to verify the gap-filling method [35,51,52]. The idea of the cloud assumption is to use clouds to mask the known SCE and then compare the recovered SCE with the known SCE to assess the performance of the cloud removal results.

3.2.2. Validation Based on Landsat–8 OLI Binary Snow Mapping

Binary snow maps retrieved from Landsat–8 OLI were used to verify the produced cloud-free MODIS SCE products. The SNOMAP algorithm defined a snow pixel as having NDSI ≥ 0.4 and the reflectance of near infrared (NIR) band ≥ 0.11 [47,53]. Due to the forest canopy obstructing the snow information from being transmitted to the sensor, the traditional method of snow cover mapping has low accuracy in a snow-covered forest [54,55]. To perform high-accuracy snow mapping in forest areas, a multi-index method focuses on snow-covered forests with NDSI < 0.4, making use of the normalized difference forest snow index (NDFSI) > 0.4 and the normalized difference vegetation index (NDVI) < 0.6 to extract forest snow; this method has an accuracy of more than 90% in northern Xinjiang, China [56]. Therefore, we used this method to extract the SCE from Landsat–8 OLI images and defined snow pixels as having NDSI ≥ 0.4 and NIR band≥0.11 or NDSI < 0.4, NDFSI > 0.4 and NDVI < 0.6. The normalized difference indexes are defined as follows:

$$\mathrm{NDSI}=\frac{{\rho }_{Green}-{\rho }_{SWIR}}{{\rho }_{Green}+{\rho }_{SWIR}}$$

(5)

$$\mathrm{NDFSI}=\frac{{\rho }_{NIR}-{\rho }_{SWIR}}{{\rho }_{NIR}+{\rho }_{SWIR}}$$

(6)

$$\mathrm{NDVI}=\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$$

(7)

where ${\rho }_{SWIR}$, ${\rho }_{NIR}$, ${\rho }_{Green}\mathrm{and}$ ${\rho }_{Red}$ represent the reflectance of the shortwave infrared (SWIR), near infrared (NIR), green, and red bands in the Landsat–8 OLI image.

After that, we aggregate the Landsat–8 binary snow cover from 30 m to MODIS 500 m resolution to produce a fractional snow cover product, those pixels with fractional snow cover $\ge$ 15% are defined as snow, and others are snow-free [25,57,58,59]. Then, we validated the accuracy of the produced daily cloud-free MODIS SCE products.

3.2.3. Accuracy Assessment Metrics

The basic evaluation metrics based on the confusion matrix (

Table 1. Confusion matrix.

	Results: Snow	Results: Snow-Free
Truth: snow	a	b
Truth: snow-free	c	d

) were applied to evaluate the results of the produced daily cloud-free MODIS SCE products; the metrics included the overestimation error (OE), underestimation error (UE), overall accuracy (OA), and F-score (Fs) [52,60]. OE is the probability that the recovered pixel is snow but the corresponding pixel in the true image is snow-free; UE is the probability that the recovered pixel is snow-free but the true pixel is snow; OA is defined by the number of correctly recovered pixels divided by the total number of pixels; Fs penalizes both the overestimation and underestimation of snow without the influence of extensive snow-free areas. The formulas for OE, UE, OA, and Fs are defined as follows:

$$\mathrm{OE}=\frac{\mathrm{c}}{\mathrm{a}+\mathrm{b}+\mathrm{c}+\mathrm{d}}$$

(8)

$$\mathrm{UE}=\frac{\mathrm{b}}{\mathrm{a}+\mathrm{b}+\mathrm{c}+\mathrm{d}}$$

(9)

$$\mathrm{OA}=\frac{\mathrm{a}+\mathrm{d}}{\mathrm{a}+\mathrm{b}+\mathrm{c}+\mathrm{d}}$$

(10)

$$\mathrm{Fs}=\frac{2\times \mathrm{a}}{2\times \mathrm{a}+\mathrm{b}+\mathrm{c}}$$

(11)

3.3. Analysis of the Snow Cover Variation

Snow cover days (SCD) were calculated through the generated cloud-free MODIS SCE products, and the SCD is defined as the total number of days on which each pixel has snow in a given year. To capture the possible positive and negative trends of the snow cover variation and significance level, the Mann–Kendall and Theil–Sen median methods were conducted on each SCD series.

The Mann–Kendall method is a nonparametric statistical test method that is widely used to judge the significance of snow cover variation. This method does not need data that obey normal distribution, and it is not affected by a few outliers [6,61,62]. For a series Xi = (X₁, X₂,⋯ ⋯, X_n) with n samples, the specific formula is expressed as follows:

$$\mathrm{Z}=\{\begin{array}{c}\frac{\mathrm{S}-1}{\sqrt{\mathrm{Var}\left(\mathrm{S}\right)}}, if S>0\\ 0, if S=0\\ \frac{\mathrm{S}+1}{\sqrt{\mathrm{Var}\left(\mathrm{S}\right)}}, if S<0\end{array},$$

where

$$\begin{gathered} \mathrm{S}={{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=i+1}^{n}sgn(X_j-X_i), \\ sgn\left(X_j-X_i\right) =\{\begin{array}{c}+1, if\left(X_j-X_i\right)>0\\ 0, if\left(X_j-X_i\right)=0\\ -1, if\left(X_j-X_i\right)<0\end{array} \\ \mathrm{Var}(\mathrm{S})=\frac{\mathrm{n}\left(\mathrm{n}-1\right)\left(2\mathrm{n}+5\right)-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{t}}_{\mathrm{i}}({\mathrm{t}}_{\mathrm{i}}-1)\left(2{\mathrm{t}}_{\mathrm{i}}+5\right)}{18} \end{gathered}$$

(12)

Additionally, n, m, and t_i are the number of years, tied groups, and data values in the p-th group, respectively. The test statistic Z was used for the two-sided trend test when n > 10. Additionally, Z > 0 represents an upward trend, while Z < 0 represents a downward trend. As for significance level α, |Z| > Z_1−α/2 represents a significant trend in the time series, otherwise it is nonsignificant; Z_1−α/2 is acquired from the standard normal cumulative distribution table.

Theil–Sen median trend analysis is an effective nonparametric method for statistical trend variation, and the result is not greatly affected by abnormal values [63]. The calculation process of the Theil–Sen median slope for a series of length n is defined as follows:

$$\mathsf{\beta }=\mathrm{Median}\left(\frac{X_j-X_i}{\mathrm{j}-\mathrm{i}}\right), j>i$$

(13)

where β > 0 and β < 0 represent an increasing and decreasing trend, respectively.

4. Results

4.1. Evaluation of the Gap-Filling Methods

Based on the gap-filling method proposed in this paper, daily cloud-free MODIS SCE products for 2001–2018 in northern Xinjiang, China can be produced. Figure 4 shows the monthly variation in the cloud fraction (CF) for Terra, Aqua, and MOYD. For the original MODIS data, the CF of the Aqua products is similar to that of the Terra products. After TAC, the CF of MOYD is reduced by 8.99% on average compared to the Terra products. However, the remainder of the average CF is still 38.67%. Following STCPI, all the residual clouds are removed completely, and this method is effective in filling large cloud gaps. TAC has been extensively proven to have high accuracy, so we only evaluate the accuracy of STCPI for the generated daily cloud-free MODIS SCE products. The cloud assumption and Landsat–8 validation are implemented in the following.

4.1.1. Validation Based on the Cloud Assumption

To evaluate the effectiveness of STCPI, the MOYD images with the lowest CF were selected in every month in 2018 as the true images, and a CF close to 25%, 50%, and 75% for every month in 2017 is applied to mask the true images. Figure 5 shows the STCPI prediction accuracy of the selected images. The average values of overestimation error, underestimation error, overall accuracy, and F-score are 1.19%, 1.37%, 97.44%, and 86.76%, respectively.

From Figure 5, we can draw the following four conclusions: (1) STCPI can effectively fill the gaps in binary snow cover and maintain high accuracy because OA maintains a high value all year. (2) The average OE is close to UE, which shows that the STCPI method does not lead to an overall overestimation or underestimation of the SCE. (3) Fs in summer is lower than that in winter, this is because the snow fraction of northern Xinjiang, China, is very small in summer, such as the snow fraction of only 2.35% on June 06, 2018. When there is little snow cover in summer and the image is covered by clouds, there is no effective neighborhood information that can be used to recover the snow information under the cloud coverage. This small part of the snow information recovery error does not affect the overall snow recovery results, which still have a high OA value. (4) A greater amount of CF shows better accuracy, which is related to the different spatial distribution of cloud-free pixels and the total number of cloud pixels for calculating accuracy. The recovered the snow condition of the pixel masked by cloud based on the neighborhood spatiotemporal information is more likely to cause misclassification at the boundary between snow and snow-free pixels. The larger cloud fraction usually includes smaller proportion of boundary pixels, and so a better precision can be achieved. Meanwhile, the gap-filling accuracy of STCPI is similar for different CF, which indicates that the STCPI method can remove large gaps in binary snow cover effectively. Figure 6 shows an example: STCPI removes CFs of 25.06%, 57.77%, and 77.38% and has an OA value of 97.72%, 98.41%, and 98.43%, respectively. Compared with the true images, the snow cover information is perfectly recovered visually.

4.1.2. Validation Based on Landsat–8 OLI

Validation based on the cloud assumption can effectively evaluate the accuracy of the gap-filling method, but the accuracy of the generated daily cloud-free MODIS SCE products is influenced by the gap-filling method and the original MODIS snow classification method, so we further evaluate the daily cloud-free MODIS SCE products using Landsat–8 OLI. OE, UE, OA, and Fs were calculated to validate the accuracy of the daily cloud-free MODIS SCE products acquired in this study.

We selected a total of 20 Landsat–8 OLI images (CF < 5%) in northern Xinjiang, China, from 2018 for validation.

Table 2. Landsat–8 OLI validation of the performance of cloud-free MODIS SCE products.

Landsat–8 Path/Row	Acquisition Day	CF of MOYD (%)	OE (%)	UE (%)	OA (%)	Fs (%)
142/30	7 January 2018	26.08	13.21	1.18	85.61	88.79
143/26	30 January 2018	29.74	9.75	2.70	87.55	92.56
143/30	30 January 2018	33.79	13.45	5.69	80.86	87.60
147/30	10 January 2018	22.33	10.14	1.80	88.06	93.28
141/28	17 February 2018	65.99	5.46	4.73	89.81	93.74
141/30	1 February 2018	18.45	7.47	0.74	91.79	91.55
145/26	13 February 2018	8.54	4.10	1.88	94.02	96.41
145/27	13 February 2018	22.33	5.24	4.56	90.20	92.10
143/26	3 March 2018	31.82	2.25	3.56	94.19	95.65
142/27	13 April 2018	5.87	1.49	3.57	94.94	94.82
145/26	27 October 2018	4.41	2.66	1.53	95.81	97.01
145/30	27 October 2018	8.65	7.32	6.56	86.12	88.66
141/29	16 November 2018	46.85	4.65	0.89	94.46	94.63
143/30	14 November 2018	31.03	9.78	3.43	86.79	90.56
147/29	26 November 2018	45.48	6.31	3.56	90.13	91.71
147/30	26 November 2018	33.36	5.83	1.11	93.07	95.14
143/27	16 December 2018	28.56	6.48	1.20	92.32	95.12
144/28	7 December 2018	45.53	4.66	1.42	93.91	92.59
145/27	14 December 2018	46.73	6.64	2.33	91.03	93.85
145/30	14 December 2018	19.01	12.26	1.64	86.09	91.63
Mean value		28.73	6.96	2.70	90.34	92.87

shows the accuracy of the four assessment indicators and the CF of MOYD. We found that the average OE, 6.96%, of the daily cloud-free MODIS SCE products is higher than the UE, 2.70%, based on Landsat–8 OLI validation, which may be affected by the accuracy of the binary snow cover generated by the NDSI threshold. However, on the whole, the generated SCE has a high agreement with the Landsat–8 OLI image, with OA and Fs values of 90.34% and 92.87%, respectively.

4.2. Snow Cover Variability

Based on the generated daily cloud-free MODIS SCE products, we analyzed the variability of snow cover over northern Xinjiang, China, from 2001 to 2018. The main characteristics include the intra-annual variability, spatial distribution, and interannual variability.

4.2.1. Intra-Annual Variability of Snow Cover

The SCA is widely used to analyze intra-annual snow cover changes. Figure 7 shows the time series of the daily SCA (%) for the whole region and six altitudinal zones from 2001 to 2018. The SCA in the study area has a strong seasonal variation. For all of northern Xinjiang, China, the SCA from December to February is greater than 60% as a stable snow period; in the March–April period, the SCA progressively decreases as snow melts; in the May–August period, the SCA is less than 5%; and in the September–November period, the SCA progressively increases as snow accumulation. From the perspective of the change in the standard deviation, the value in the May–August period is low, so the interannual fluctuation is slow, and the change is relatively large the rest of the time (Figure 7a).

At different altitudinal zones, the difference in SCA in terms of seasonal fluctuation shows obviously climatic characteristics (Figure 7b). On the whole, with the increase in elevation, the SCA increases gradually, the snowmelt period moves back, and the snow accumulation period moves forward. In the subregion with elevation > 3000 m, the lowest SCA is 23.88%, even in summer. Furthermore, 69.70% of northern Xinjiang, China, has an elevation < 1500 m, and the SCA variation in these areas has an obvious seasonal characteristic: the snow fraction is more than 80% in winter and close to 0 in summer.

4.2.2. Spatial Distribution of Snow Cover

The mean SCD from 2001 to 2018 reflects the spatial distribution of the snow cover in northern Xinjiang, China (Figure 8). The average SCD and elevation have a positive correlation (r = 0.76), which shows that the SCD gradually increases as the elevation increases. SCD > 60 days, which is regarded as the stable SCA, occupies 96.64% of northern Xinjiang, China. Compared with the stable SCA, the unstable SCA, with SCD ≤ 60 days, occupies only 3.36%. The area with SCD > 200 days, occupying 10.85%, is mainly distributed in the Altai Mountains and the Tianshan Mountains.

4.2.3. Interannual Variability of Snow Cover

We used the Mann–Kendall and Sen’s median methods to analyze the trend and significance of SCD for northern Xinjiang, China, from 2001 to 2018 at the pixel scale. As shown in Figure 9, the SCD in 34.39% of the area decreased (Z < 0) during the 18 years and was mainly distributed in the Junggar Basin. The area of decrease was significant in 0.86% of the total area (Z < −1.96). In addition, the area with increasing SCD accounted for 59.96% of the total area, and 3.80% of the area experienced a significant increase; this region was mainly distributed in the Altai Mountains and the Tianshan Mountains. On the whole, SCD has a slight increasing trend in northern Xinjiang, China, from 2001 to 2018.

The slope of the SCD variation in northern Xinjiang, China, is shown in Figure 10. The SCD variation within 2 days/year has a proportion of 97.75%. In particular, the SCD shows an increased trend (β > 0) over 57.21% of the area, and the SCD shows a decreasing trend (β < 0) over 31.88% of the area. This also shows that the SCD has a slightly increased trend in northern Xinjiang, China, from 2001 to 2018, which is similar to the analysis result of the Mann–Kendall method. The increasing tendency of the SCD is consistent with previous studies that show snowfall has significantly increased in northern Xinjiang, China in the last few decades [12,64,65].

5. Discussion

MODIS version 6 snow products are used to generate daily cloud-free SCE products, and the accuracy of version 6 products is greatly improved compared with version 5 products due to better handling of atmospheric correction, cloud and snow misclassification, recovering Aqua band 6 [48]. Meanwhile, the original NDSI threshold of 0.4 applied in version 5 products has been proved by many studies to cause serious underestimation of snow cover in forest areas [55,66,67], which will reduce the credibility of the research on the spatial and temporal variation of snow cover. Based on MODIS version 6 products and NDSI threshold of 0.1, the accuracy of snow products is greatly improved. Therefore, a new cloud removal algorithm of STCPI is used to generate higher accuracy daily cloud-free MODIS SCE products.

STCPI can effectively fill large cloud gaps to generate SCE scenes, which may be attributed primarily to the following two advantages: (1) this method makes full use of adjacent space-time information; the spatial filter and temporal filter method uses spatiotemporal information step by step, which loses a large amount of information in the space-time cube, so it does not completely fill the cloud gaps in the image [13,68]. Conversely, STCPI makes full use of adjacent space-time information that can completely fill cloud gaps in one step and maintain high accuracy even with large cloud gaps in the image. (2) STCPI first calculates the conditional probability of the central pixel and every neighboring pixel in the space-time windows has the same snow condition. Generally, the closer the conditional probability is to 1, the more consistent the snow conditions of the neighboring pixels and central pixels [40]. The conditional probability established for fixed-position pixels is equivalent to considering the land cover, elevation, local environment information, and so on in the spatial neighborhood. Additionally, only a conditional probability > 0.8 was selected to calculate the snow probability of the central pixel because many neighboring pixels have different snow conditions from the central pixel, especially in mountain areas.

In addition, the STCPI shows that Fs in summer is lower than in winter in the validation of the cloud assumption. A difficult problem in image restoration is when there is little snow cover in summer and is covered by clouds [69], because there is no effective neighborhood information that can be used to recover the snow conditions under cloud coverage. This problem needs further study in the future. Moreover, the accuracy of the daily cloud-free MODIS SCE products is commonly influenced by the gap-filling method and the original MODIS snow classification method. Therefore, a higher-precision snow classification method can effectively improve the quality of the final daily cloud-free MODIS SCE products.

Under the background of global warming, snow cover change has a significant impact on social economy and regional climate. Therefore, it is very important to study the temporal and spatial changes of snow cover in northern Xinjiang, China. The snow cover has obvious seasonal variation trend in this region. The amount of snow cover is large in winter but little distribution in high altitude mountainous areas in summer. Chen et al. [43] and Li et al. [68] analyzed the relationship between snow cover and altitude, and showed that the distribution of snow cover had a clear positive correlation with altitude, we also obtained an correlation coefficient 0.76 in northern Xinjiang, China. Based on multivariate data fusion method among MODIS version 5 snow products, the advanced microwave scanning radiometer for earth observing system daily snow water equivalent products and IMS products, the analyses of the snow cover changes in northern hemisphere show that SCD of Altay Mountain and Tianshan Mountain is significantly reduced [32]. Tang et al. [44] generated cloud-free MODIS fractional snow cover products based on cubic spline function interpolation method and analyzed the SCD variation of Tianshan Mountain, they found that 26.39% of that region shows a downward trend and 34.26% shows a upward trend from 2001 to 2015. It is rare to use MODIS snow products to research the spatial and temporal variation of snow cover in northern Xinjiang, China. Yang et al. [42] and Wang et al. [12] analyzed the spatial-temporal changes of snow cover using the site data, and obtained the conclusion that the snowfall increased in northern Xinjiang, China, in the past decades, which is consistent with the slight increase trend of SCD in northern Xinjiang, China in this study.

6. Conclusions

This study proposed a new gap-filling method that produced high-precision daily cloud-free MODIS SCE products for northern Xinjiang, China, from 2001 to 2018 through TAC and STCPI. Through the validation of the cloud assumption, the OA of STCPI reached 97.44%. The OA of the final daily cloud-free MODIS SCE products is 90.34%, as validated by Landsat–8 OLI.

Based on the generated daily cloud-free MODIS SCE products, the snow cover variability was analyzed. The spatiotemporal variation characteristics of the snow cover from 2001 to 2018 provide further evidence for climate change in northern Xinjiang, China. It was found that the interannual change of SCA gradually decreases as the elevation increases and that the SCD has a positive correlation with elevation. Furthermore, based on Mann–Kendall and Sen’s median methods, we analyzed the changing tendency of the interannual snow cover. The results show that the area of SCD increase is higher than that of decrease during the 18 years. In the future, the relationship between snow cover variation and temperatures needs to be analyzed in detail under the conditions of a warming climate.