Gap-Filling of MODIS Fractional Snow Cover Products via Non-Local Spatio-Temporal Filtering Based on Machine Learning Techniques

Abstract

Cloud obscuration leaves significant gaps in MODIS snow cover products. In this study, an innovative gap-filling method based on the concept of non-local spatio-temporal filtering (NSTF) is proposed to reconstruct the cloud gaps in MODIS fractional snow cover (SCF) products. The ground information of a gap pixel was estimated by using the appropriate similar pixels in the remaining known part of an image via an automatic machine learning technique. We take the MODIS SCF product cloud gap filling data from 2001 to 2016 in Northern Xinjiang, China as an example. The results demonstrate that the methodology can generate almost continuous spatio-temporal, daily MODIS SCF images, and it leaves only 0.52% of cloud gaps long-term, on average. The validation results based on “cloud assumption” exhibit high accuracy, with a higher $R^2$ exceeding 0.8, a lower RMSE of 0.1, an overestimated error of 1.13%, an underestimated error of 1.4%, and a spatial efficiency (SPAEF) of 0.78. The validation based on 50 in situ snow depth observations demonstrates the superiority of the methodology in terms of accuracy and consistency. The overall accuracy is 93.72%. The average omission and commission error have increased approximately 1.16 and 0.53% compared with the original MODIS SCF products under a clear sky term.

1. Introduction

Snow cover is one of the fastest changing objects on the Earth’s surface. It strongly affects the climate system, radiation, and energy balance between the atmosphere and the Earth, hydrological circulation, biogeochemical cycling, and even human activities [1,2]. In the Northern Hemisphere, the snow coverage ranges from approximately 45.2 to 1.9 million square kilometers, which occurs in January and August, respectively. The drastic variation of snow cover spatial extent is a sensitive indicator of climate change [3]. The melting of seasonal snow cover supplies the dominant water source for river runoff and ground water discharge, particularly for arid and semi-arid regions [4]. Approximately 1/6 of the world population and more than 1/5 of the Chinese population depend on water resources supplied by snow cover or glaciers [5]. Snow cover is also a critical parameter in hydrological and ecological process models [6]. There is an urgent requirement to accurately map the spatial distribution of snow cover. This is not only because of a rising demand for fresh water management but also due to the concerns about the effects of climate change.

Several remotely sensed snow cover products provide the long series of global/regional snow cover information [7]. The suite of MODIS snow cover products, available since 2000, has been applied widely in mapping snow cover area because of its moderate spatial resolution and high temporal resolution [8,9]. Extensive studies have demonstrated that the standard MODIS snow products have high snow classification accuracy (>90%) under a clear sky term. The evaluation results of MODIS snow cover product report that it has an average classification accuracy of 93% for MOD10A1 in the Upper Rio Grande Basin [10]; 95% accuracy for eight-days composited MOD10A2 in Northern Xinjiang, Northwest China [11], an accuracy exceeding 90% for MOD10A1 in the Qinghai-Tibet Plateau [12,13], an accuracy of 93% for MOD10A1 in Canada [14], an accuracy of 95% for MOD10A1 in Europe [15] and an accuracy exceeding 92% for MOD10A1 and MYD10A1 in Central Asia [16].

Nevertheless, cloud contamination in MODIS snow cover products has inevitably caused great data gaps which have seriously hindered its application in total environment research. Lin et al. (2013) suggested that the global annual mean cloud cover was approximately 66% from the project of International Satellite Cloud Climatology [17]. Parajka et al. (2006) demonstrated that the clouds have obscured, on average, 63% of MODIS snow cover maps in Austria from 2000 to 2005 [15]. Similar average cloud coverages of approximately 70% in the Quesnel River Basin and 50–60% in Alaska are also reported [18,19]. The mean cloud coverage of MODIS snow cover products was 46.86% (Terra) and 48.47% (Aqua) in Europe from 2000 to 2011 [20], 39% in the Upper Salt River Basin from 1 October 2004 to 31 May 2005 [21]. The mean overall classification accuracy of MODIS snow cover products was approximately 34–50% in all weather conditions [22,23,24]. Wang et al. (2009) verified that the mean annual cloud coverage of MODIS snow cover products was 44% (Terra) and 47% (Aqua) in Northern Xinjiang, China. The overall classification accuracy was below 40% in all-sky conditions [25].

A series of efforts have been made to fill the cloud gaps in the MODIS snow cover maps. The combination of MODIS Terra and MODIS Aqua snow cover products can reduce approximately 10–20% of cloud obscuration and achieve an accuracy of 90% [21,25]. Temporal filtering directly replaces cloud pixels by using ground information of a pixel in the previous few days or/and subsequent days [26]. Previous research has demonstrated that seven days of temporal filtering can eliminate more than 95% of clouds while maintaining an overall accuracy exceeding 92% over Austria [27]. Spatial filtering is designed to replace cloud-contaminated pixels by using information on surrounding adjacent pixels that are not obscured by clouds on the same date, which can fill gaps by approximately 6–13% [23,27,28,29,30]. SNOWL is also a spatial filtering method based on snow transition elevation, which is capable of reducing cloud coverage in standard MODIS snow cover products from 60 to 20% [30,31]. In practice, the combined use of temporal and spatial filtering is the most common way for MODIS cloud mitigation, which can produce cloudless images as well as maintain the accuracy of the source products under clear-sky conditions. The combination of MODIS and passive microwave (e.g., AMSR-E) snow cover products take advantage of MODIS’s high spatial resolution and microwave penetration into the cloud, which is capable of eliminating the cloud-contaminated pixels but has a reduced accuracy due to the low spatial resolution of the passive microwave [23,24,32].

The additional limitation of existing methods often leaves out a percentage of cloud pixels or has low effectiveness when cloud obstruction fills the same area for several continuous days or is widely distributed. Furthermore, these studies were only applicable to MODIS binary snow products, and incapable of delineating the fractional snow cover (SCF) of a pixel, which might be more valuable for hydrological and ecological process models [33]. Dozier et al. (2008) first proposed the notion of treating the time-varying snow cover as a space-time cube. They filled the gaps caused by clouds via cubic spline temporal interpolation at each grid [34]. Tang et al. (2013) also used this cubic spline interpolation algorithm to eliminate the cloud-contaminated pixels in MODIS SCF products of the Tibetan Plateau [12]. This one-dimensional temporal interpolation generates a large number of outliers in the case of continuous cloud coverage whereas the spatial distribution characteristics of snow cover was not considered.

Cheng et al. (2014) developed a new approach based on the idea of similar pixel replacement [35]. The spatio-temporal Markov random fields (STMRF) function was constructed to search for the most appropriate similar non-cloud pixels to replace a cloud-contaminated pixel in MODIS land surface temperature products. Inspired by this work, we presented an innovative gap-filling method, based on the concept of the non-local spatio-temporal filter (NSTF), to eliminate the cloud contamination in daily MODIS SCF products. The ground information of a cloud pixel is estimated through selected similar pixels in the remaining known region. We anticipate that twofold goals can be achieved: (1) reconstructing the ground information of MODIS SCF products obscured by cloud cover and (2) achieving the overall accuracy of MODIS SCF source products under clear a sky term after filling the gaps.

2. Study Area and Data

2.1. Study Area

North Xinjiang, located between 42–50°N and 79–92°E and covering a total area of 0.39 million km² (Figure 1), is one of the three major seasonal snow cover areas in China. The area is characterized by distinct topographic relief, with the Altai Mountains in the north, the Tianshan Mountains in the south, and the Junggar basin in between. The elevation of the study area ranges from 187 to 5950 m. The land cover, which can severely influence the accuracy of snow cover mapping [8], mainly consists of grasslands and shrubs (42.6%), as well as bare ground (40.7%). Cropland (7.7%), water bodies (2.5%), urban and built up land (0.6%), and forests (5.9%) cover only 16.7% of the total area [34]. Xinjiang has a temperate continental climate and is a typical arid and semi-arid region. The annual mean precipitation is approximately 145 mm, and the evaporation capacity is about 200 mm [36]. The study area has a long, cold winter and a short, hot summer, with an extreme minimum temperature of −40 °C. Over half of the precipitation in the study area falls as snow in the winter. The seasonal snow cover period is approximately 120 days, with a mean snow depth of 0.6 m [25]. The considerably heterogeneous nature of the topography, diverse land cover types, and changeable climatic conditions have created complex spatial and temporal snow cover patterns.

2.2. Data

2.2.1. MODIS SCF Products

MODIS Terra and MODIS Aqua daily SCF products (MOD10A1 and MYD10A1, Collection 005) from 1 January 2001 to 31 December 2016 and 4 July 2002 to 31 December 2016 were collected from NASA’s EOSDIS “Earthdata” (https://earthdata.nasa.gov/), respectively. Three tiles (i.e., h23v03, h23v04, and h24v04) were required to cover the whole study area. The MODIS Reprojection Tool (MRT) was used to mosaic each of the three tiles via nearest neighbor resample method and reprojected the primitive sinusoidal projection into the Universal Transverse Mercator (UTM zone 45) projection with the WGS84 datum and 500 m resolution.

2.2.2. Meteorological Observation, Land Cover and Digital Elevation Model data

The daily snow depth (SD) recorded by 50 meteorological stations were provided by the Chinese Meteorological Data Sharing Service System (CMDS) (http://cdc.cma.gov.cn/). Most of these meteorological stations are located in the inhabited and low elevation valleys. Only three stations are located at elevations above 2000 m (Figure 1). The SD observations were recorded in centimeters, and the fractional part is rounded up.

The land cover data was extracted from land cover products of China, which was provided by the Cold and Arid Regions Science Data Center at Lanzhou (http://westdc.westgis.ac.cn). The overall classification accuracy of this land cover map was 71%, which was higher than the accuracies of other land cover maps (i.e., global land cover dataset of IGBP, MODIS land cover map, a global land cover map produced by the University of Maryland, and the land cover map produced by the global land cover for the year 2000 (GLC 2000) project. The overall classification accuracy of these four land cover maps is between 54.17–59.28%) [37,38]. The land cover map uses a unified 1000 m resolution and ALBERS projection. We reprojected it into UTM projection with the WGS84 datum and resampled it into 500 m resolution so that it has consistent projection and resolution with MODIS SCF data.

The NASA Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) data (SRTM3, 90 m) was provided by the Consortium for Spatial Information (CGIAR-CSI) (http://srtm.csi.cgiar.org). The eighteen DEM tiles are mosaicked, reprojected, and resampled to achieve consistency with the MODIS SCF images. The elevation, aspect, and slope maps were extracted from processed DEM data.

3. Methodology

3.1. Investigation on Cloud Gaps in MODIS SCF Products

How much cloud gaps are in the MODIS SCF products? How long is the cloud cover duration? To answer these two questions, we first reclassified the MODIS SCF products into two categories: fractional snow (0 to 100) and cloud (250). After reclassification, annual cloud coverage is calculated to analyze the influence of cloud gaps on MODIS SCF products. Furthermore, to explore the relationship between elevation and cloud cover frequency, we also take elevation into account in our calculations. For example, the daily cloud coverage for elevations above the specified threshold $Z_T$ of day $d$ can be calculated using the following formula.

$${\mathrm{C}}_{\mathrm{d}}\left(Z>Z_T\right)=\frac{N_d\left(Z>Z_T\right)}{N_{Tol}\left(Z>Z_T\right)}$$

(1)

where $N_d\left(Z>Z_T\right)$ and $N_{Tol}\left(Z>Z_T\right)$ are the numbers of cloud pixels and total number of pixels in the elevation above $Z_T$ for date $d$ in the study area. If $Z_T$ is taken as the minimum elevation ($Z_{min}$) of study area, the study area is considered in its entirety.

3.2. Gap-Filling Procedure

3.2.1. Theoretical Basis

Existing studies have found that there are always numerous similar pixels in a remote sensing image. Similar pixels do not exist only in an adjacent location but also appear in non-local neighboring regions [35,39]. Similar pixels have the following characteristics. (1) They have approximately equal spectral characteristics and, (2) despite images obtained at different times having different spectral characteristics, similar pixels have approximately the same reflectance changes over time and relative coincident positions in multi-temporal images [40]. For example, Figure 2 shows the snow distribution of MOD10A1 in a small region (the size is $31\times 31$ pixels, with the upper left corner is located at (45.05°N, 86.00°E)) of study area on (a) 30 January 2003 and (b) 18 February 2003. The pixels displayed in the red box in Figure 2a are the similar pixels of the target pixel (black box), which are also similar pixels of target pixels in Figure 2b. The land cover types of these pixels are cropland, and these pixels have close reflectance within an image and consistent changes with time. Taking this basic concept, the ground information of a cloud pixel within a MODIS SCF image can be estimated by using the appropriate similar pixels in the remaining known region with the help of multitemporal reference images.

However, some other studies point out that a pixel would also show considerable reflectance change where land cover type changes, therefore causing the phenomenon of temporal difference measurement [41,42]. For example, regarding the pixel we picked out with a green square in Figure 2a, its land cover type is bare land. It is also the similar pixel of the target pixel before snow melts in Figure 2a, but it is no longer the similar pixel of target pixels in Figure 2b after partial snow melting. This may be caused by different snow melting rates under changed underlying surface conditions. To tackle this problem, the pixel which had the same land cover type with the target pixel was regarded as a similar pixel [41,42,43].

3.2.2. Gap-Filling Method

An innovative gap-filling method via NSTF is proposed to reconstruct the ground information of the cloud gaps in the daily MODIS SCF products. As there are significant gaps in large cloud regions, it is unreliable and unrealistic to rely on the remaining known regions of the image to match similar pixels. Therefore, a combination of Terra and Aqua, and temporal filtering are executed first to fill some gaps before NSTF. However, the NSTF method can be used independently. The detailed schematic of the MODIS SCF products gap-filling procedure is presented in Figure 3.

1. Daily MOD10A1 and MYD10A1 Combination

Daily MOD10A1 and MYD10A1 products were combined first by using the composite rules where priority was assigned to MOD10A1, as the validation studies demonstrated that it had more accurate retrievals than MYD10A1 [13]. Hereafter we refer to MOD10A1, MYD10A1, and the combination products to MOD, MYD, and MOYD, respectively. The specific combination strategy is the following. When a pixel is cloud-free in both MOD and MYD, the value in MOD is taken. When a pixel is cloud-free in only one of the products, the cloud free observation is taken. In addition, because the MYD product was not available until July 2002, we only used the MOD product prior to that date.

2. Adjacent Temporal Filtering

Adjacent temporal filtering (ATF) rests on the fundamental assumption that snow will persist on the surface for a certain period of time. Previous studies of one-day to eight-day temporal filters have confirmed that the overall accuracy of multi-days combined snow cover products improves with the increase of composition days [26,27,32,44,45]. However, Gao et al., (2010) found that the overall accuracy of multi-day combined snow cover products decreases with increasing composition days [19]. To determine the optimum composition days, we calculated a cloud persistence index (CPI) for each pixel via finding the average of the continuously cloud covered days [25]. The result indicates that three days is the optimum number of composition days, as the mean CPI for the study area from 2001 to 2016 is 2.72 days. In addition, to guarantee the real-time application ability of the proposed method, we only use images from the previous date. Thus, we determined the three-days backward temporal filtering method used in this study. For a cloudy pixel, we searched images from the previous three days and, if a cloud-free observation was obtained from two or more days, the last valid SCF at this location was assigned to this cloud pixel.

3. NSTF

To match the similar pixels in the remaining non-cloud covered regions to fill the gap pixel, the following three criteria need to be considered. (1) The similar pixels have similar SCF values and have an approximately similar SCF variation in the time domain; (2) Considering the continuity of snow distribution in the spatial domain, the similar pixels have an approximately similar SCF proximity with their respective surrounding neighbors; (3) The similar pixels have approximately similar topographic features. According to the above three criteria, three objective functions are defined. The similar pixels are the solution of this three-objective optimization problem which requires the simultaneous minimization of three objective functions under equality constraints. However, the three objective functions are almost impossible to minimize simultaneously because there is often conflict between multiple objectives. To provide a means to assess trade-off between three conflicting objectives, an automatic machine learning technique of the fast elitist and non-dominated sorting multi-objective genetic algorithm (NSGA-II) is adopted [46,47].

First, for a cloud pixel $P$ on target image, it is assumed that $P^{\prime }$ is the pixel at the same position as $P$ on an auxiliary reference image. The historical image closest to the target image with $B\left(P^{\prime }\right)$ is a $3\times 3$ non-cloud block centered at $P^{\prime }$ (Figure 4). This was selected as the most appropriate auxiliary reference image. We then extracted the $M$ non-cloud common pixels of the target image and reference image within a 51 × 51 pixels search window centered at location $P$, as is too time-consuming to search for similar pixels in the entire study area [41]. To guarantee the precision of the NSTF, we let $M$ account for at least 30% of the total number of pixels in searching window. If this condition was not satisfied, the search window was expanded to 101 × 101 pixels then 151 × 151 pixels, etc., until the entire study area is encompassed.

Second, the three-objective MODIS SCF products gap-filling problem can be expressed as:

$$Minimize:\begin{array}{c}\left[f_1\left(P_S\right),f_2\left(P_S\right),f_3\left(P_S\right)\right]\end{array}$$

(2)

$$Subject\begin{array}{c}to:\end{array}LC(P_S)=LC(P)$$

where $P_s$ is a candidate similar pixel of $P$. $LC$ is the land cover type.

A: Minimize Variation in Time Domain

The objective considers that similar pixels having similar change trends over time, and relative coincident positions in multi-temporal images can be mathematically defined as

$$f_1\left(P_S\right)={\Vert S_T\left(P_S\right)-S_T\Vert }^2$$

(3)

where the subscript T represents the target image, $S_T\left(P_S\right)$ is the SCF value of $P_S$, $S_T$ is the average SCF value of pixels in a corresponding location ${P^{\prime }}_S$ in the target image, and ${P^{\prime }}_S$ is a candidate similar pixel of $P^{\prime }$ which was calculated in advance (Figure 4). The calculation procedure of ${P^{\prime }}_S$ is

$${\Vert S_R\left(B\left(P^{\prime }\right)\right)-S_R\left(B\left({P^{\prime }}_S\right)\right)\Vert }^2\le d\times 2^{S_R\left(P^{\prime }\right)}$$

(4)

$$Subject to: LC({P^{\prime }}_S)=LC(P^{\prime })$$

where $B\left({P^{\prime }}_S\right)$ is a $3\times 3$ non-cloud block centered at location ${P^{\prime }}_S$; subscript $R$ represents the reference image; $S_R\left(P^{\prime }\right)$ is the SCF value of $P^{\prime }$; and $S_R\left(B\left(P^{\prime }\right)\right)$ and $S_R\left(B\left({P^{\prime }}_S\right)\right)$ represent the average SCF value of blocks $B\left(P^{\prime }\right)$ and $B\left({P^{\prime }}_S\right)$, respectively. $d$ is a free parameter related to the sensor [41]; it is set as 0.01 for the MODIS case. Equation (4) is the similarity measurement to ensure that the difference between the two blocks is extremely slight.

B: Minimize Discontinuity in Spatial Domain

The objective function f₂ considers the spatial continuity of snow cover distribution. The continuity between a cloud pixel and surrounding spatial neighbors can be presented as

$$f_2\left(P_S\right)={\Vert S_T\left(P_S\right)-\overline{S_T\left(P^{″}\right)}\Vert }^2+{\Vert S_R\left({P^{\prime }}_S\right)-\overline{S_R\left({P^{″}}_{}\right)}\Vert }^2$$

(5)

where $P^{″}$ is the spatial non-cloud pixel of first order adjacency to $P$ (Figure 4), and $\overline{S_T\left(P^{″}\right)}$ and $\overline{S_R\left(P^{″}\right)}$ represent the average SCF values of $P^{″}$ in the target and reference image, respectively.

C: Minimize the Topographic Differences

Similar to [21], we also assume that there is equivalent exposure and the same amount of snow cover under the same terrain conditions. Thus, the objective function f₃ is used to enforce similar pixels having similar topographic features.

$$f_3\left(P_S\right)={\Vert E\left(P_S\right)-E\left({P^{\prime }}_S\right)\Vert }^2+{\Vert S\left(P_S\right)-S\left({P^{\prime }}_S\right)\Vert }^2+{\Vert A\left(P_S\right)-A\left({P^{\prime }}_S\right)\Vert }^2$$

(6)

where $E$, $S$, and $A$ represent normalized elevation, slope and aspect, respectively. Among them, $E$ and $S$ adopt the Min-Max normalization method, and the normalization method of $A$ was: $A=\left|A_0-{180}^{{}^{\circ }}\right|/{180}^{{}^{\circ }}$ ($A_0$ was the original aspect value).

Third, the NSGA-II algorithm was executed to search for $N$ (the value is 10 in the present study) similar pixels among the $M$ non-cloud common pixels. We encoded the chromosome ($CM$, a solution to this three-objective minimal optimization problem) as the combination of the location of the $N$ similar pixels.

$$CM=[x_1{y}_1{x}_2{y}_2{x}_3{y}_3\dots x_N{y}_N]$$

(7)

where $\left(x_i,y_i\right)$ is the location of i-th similar pixels, i = 1, 2, …, N.

The detailed execution process of NSGA-II can be referred to in [46].

Finally, the arithmetic mean SCF value of $N$ similar pixels (that is, the average SCF value of 10 pixels in corresponding locations of $\left(x_i,y_i\right)$) is arranged for cloud gap pixel $P$.

3.3. Validation Methodology

It is noticeable that the most appropriate validation technique is to use in situ SD observations as the true ground information. However, as most of the 50 meteorological stations in the area are located in the inhabited and low elevation valleys (Figure 1), the in situ SD observations can hardly provide the overall picture of the entire region because of its low spatial density and because there are no records from some inaccessible, high-mountainous areas. Therefore, another methodology based on “cloud assumption” was adopted. The cloud gaps from some more clouded MODIS SCF images were borrowed to mask the “cloud-free” MODIS SCF image, which is regarded as the true ground information [27,30,45].

3.3.1. Validation Based on “Cloud Assumption”

Considering that percentages of cloud obstruction and the duration of the clouds will affect the accuracy assessment results directly, two strategies are adopted to carry out validation [25].

1. One-Day Cloud Mask

The MODIS SCF image with the lowest cloud coverage per month from 2001 to 2016 is considered the “cloud-free” image. Then, to ensure that the artificial fictitious cloud gaps would occur as close to real life as possible, we calculated the cloud coverage level of 25, 50, and 75% probability of occurrence (PO) via the empirical cumulative distribution function (ECDF) to describe the cloud pattern of low, moderate, and high occurrence in each month. The closest images in time to 25% PO, 50% PO, and 75% PO (abbreviated to P25, P50, and P75 here after) in each month were identified, and the corresponding cloud gaps were extracted to mask the selected “cloud-free” image of this month. Consequently, each selected “cloud-free” SCF image was related to three mask images, and it co-produced a total of 576 validation images during the study period of 16 years.

2. Multi-Day Cloud Mask

To imitate the real natural conditions of consecutive overcast days, the spatio-temporal additional cloud masked images were generated. Because the mean CPI of the study area during the period from 2001 to 2016 is 2.72 days, we chose two groups of three consecutive days characterized by minimum cloud cover every snow season (i.e., November to March next year as a snow season) to be the “cloud-free” images. Additionally, the cloud gaps in the corresponding images of the previous snow cover season records were extracted to overlap onto the selected “cloud-free” images. Moreover, we chose additional groups of two, four, five, and seven consecutive days every snow season to launch more extensive validation. This co-produced a total of 90 validation groups.

Four performance evaluation criteria (i.e., root mean square error ($RMSE$), coefficient of determination ($R^2$), overestimated error ($OE$), and underestimated error ($UE$)) are defined as follows.

$$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{(x_i^s-x_i^T)}^2}}$$

(8)

$$R^2={\left(\frac{1}{N-1}{\displaystyle \sum _{i=1}^N(\frac{x_i^s-\overline{x_i^s}}{{\sigma }_s})(\frac{x_i^T-\overline{x_i^T}}{{\sigma }_T})}\right)}^2$$

(9)

$$OE=\left({\displaystyle \sum _{i=1}^N{\alpha }_i}/N\right)\begin{array}{cc}\times 100\%& \begin{array}{c}{\alpha }_i=\{\begin{array}{c}1\begin{array}{cc},& x_i^s-x_i^T\ge 5\%\end{array}\\ 0\begin{array}{cc},& x_i^s-x_i^T<5\%\end{array}\end{array}\end{array}\end{array}$$

(10)

$$UE=\left({\displaystyle \sum _{i=1}^N{\beta }_i}/N\right)\begin{array}{cc}\times 100\%& \begin{array}{c}{\beta }_i=\{\begin{array}{c}1\begin{array}{cc},& x_i^s-x_i^T\le -5\%\end{array}\\ 0\begin{array}{cc},& x_i^s-x_i^T>-5\%\end{array}\end{array}\end{array}\end{array}$$

(11)

where $x_i^s$ is the estimated SCF value of the i-th fictitious cloud pixel, $x_i^T$ is the corresponding ground truth MODIS SCF value, and $N$ is the total number of fictitious cloud masked pixels.

In addition to the four above traditional evaluation criteria, a new bias-insensitive metric called spatial efficiency (SPAEF) was used for snow cover spatial pattern validation, the SPAEF is defined as follows [48]:

$$SPAEF=1-\sqrt{{\left(A-1\right)}^2+{\left(B-1\right)}^2+{\left(C-1\right)}^2}$$

(12)

where $A$ is the correlation coefficient between the gap-filled and ground truth MODIS SCF image, $B$ is the fraction of coefficient of variation, and $C$ is the percentage of histogram intersection. The optimal value of $SPAEF$ is 1, which means that the gap-filled fictitious MODIS SCF image is exactly the same as the ground truth MODIS SCF image.

3.3.2. Validation Based on In Situ SD Observations

The original and gap-filled MODIS SCF images were compared with the in situ SD observations corresponding to the 50 meteorological stations. We applied the validation strategy based on the confusion matrix (

Table 1. Confusion matrix comparing MODIS SCF value with snow depth (SD) observation.

		Observed SD
		$\mathbf{Snow}(\ge {\mathit{\epsilon }}_1\mathbf{cm})$	$\mathbf{No}\mathbf{Snow}(<{\mathit{\epsilon }}_1\mathbf{cm})$
MODIS SCF	Snow ($\ge {\epsilon }_2$)	a	b
	No snow ($<{\epsilon }_2$)	c	d
	Cloud	e	f

), which is recommended by extensive studies such as [49].

$a$, $b$, $c$, and $d$ are the number of pixels correctly hit, commission alarmed, omission alarmed, and correctly rejected, respectively. $e$ and $f$ are the cloud pixels in the MODIS SCF image when the meteorological station reports snow and no snow, respectively. The ${\epsilon }_1$ and ${\epsilon }_2$ are the defined threshold values of observed SD and MODIS SCF to determine whether a station or a MODIS pixel is covered by snow, respectively.

The four validation indices are defined as follows.

$$\mathrm{OA}=\left(\left(a+d\right)/\left(a+b+c+d+e+f\right)\right)\times 100\%$$

(13)

$$\mathrm{MU}=\left(c/\left(a+b+c+d\right)\right)\times 100\%$$

(14)

$$\mathrm{MO}=\left(b/\left(a+b+c+d\right)\right)\times 100\%$$

(15)

$$\mathrm{Bias}=\left(\left(a+b\right)/\left(a+c+e\right)\right)\times 100\%$$

(16)

Overall accuracy ($\mathrm{OA}$) is the probability a MODIS pixel that correctly classified; $\mathrm{MU}$ and $\mathrm{MO}$ represent the underestimated and overestimated snow event, respectively; and $\mathrm{Bias}$ refers to the probability of snow event being correctly identified by MODIS. The perfect estimation of MODIS is $\mathrm{OA}=1$, $\mathrm{Bias}=1$ (we also set $\mathrm{Bias}$ to 1 when $a=b=c=e=0$), and $\mathrm{MU}={\rm M}O=0$. The imperfect estimation of MODIS is $\mathrm{OA}=0$, $\mathrm{MU}=\mathrm{MO}=1$, $\mathrm{Bias}=0$ or $\infty$ (we set $\mathrm{Bias}$ to 0 when $a=c=e=0$ and $b\ne 0$, which means no snow is observed whereas MODIS indicates snow cover).

4. Experimental Results

4.1. Cloud Gaps in MODIS SCF Products Over Northern Xinjiang

After investigating the cloud gaps in MODIS SCF products of the study area, we found that the extent of cloud obstruction is even worse than we had desired. Annual average cloud coverage of the study area is 43.53–49.85% for Terra and 42.43–52.31% for Aqua from 2001 to 2016, with a mean cloud coverage of 46.48 and 47.78%, respectively (

Table 2. Mean cloud coverage in the MODIS SCF products.

Year	Mean Cloud Coverage (%)
	MOD10A1	MYD10A1
2001	47.20	-
2002	46.95	42.43 *
2003	48.39	45.08
2004	46.52	48.37
2005	44.70	45.97
2006	44.81	46.81
2007	43.60	47.31
2008	43.53	46.06
2009	47.46	50.11
2010	48.69	49.88
2011	46.15	48.39
2012	44.21	46.40
2013	45.84	48.49
2014	46.48	48.60
2015	49.85	52.31
2016	49.33	50.54
Average	46.48	47.78

* represents only the data from 4 July 2002 to 31 December 2002 was used to calculate the annual average cloud coverage.

). The results are effectively consistent with previous studies of [24]. The cloud cover duration maps demonstrate that the mean number of cloud observations is 94–275 days for MODIS Terra SCF images and 96–316 days for MODIS Aqua SCF images, with mean cloud days of 169 and 175, respectively (Figure 5).

We took MODIS Terra SCF products as an example (though there were similar patterns for MODIS Aqua). The calculated annual cloud coverage in different elevation zones are shown in Figure 6. From Figure 6, the annual average cloud coverage for each elevation zone demonstrates that the cloud obstruction probability increases with elevation. The annual average cloud percentages are 42–48.94, 43.05–49.32, 45.07–52.12, 50.26–57.31, and 43.53–49.91 for elevations of <1000, 1000–2000, 2000–3000, >3000 m, and the entire study area, with average values of 45.31, 46.30, 48.61, 53.09, and 46.52%, respectively. The systematic trends of cloud obstruction for the entire study area are the most consistent with those of the elevation zone of 2000–3000 m. The influence of topography on cloud cover is also confirmed in Figure 5, in which the number of cloud cover days have shown a positive pattern in relation to elevation. This result is partly attributable to the water vapor uplift caused by orography, making its condensation into clouds more likely. Meanwhile, atmospheric circulation may be contributing to clearing the clouds at low-altitude regions [37].

4.2. The Effectiveness of Gap-Filling Methodology

The proposed gap-filling methodology generated almost continuous spatio-temporal daily MODIS SCF products, leaving only 0.52% of cloud gaps on long-term average. The maximum gaps remaining is 5.92%. The mean cloud coverage remaining in different elevation zones after the execution of each gap-filling step is presented in Figure 7. We found an extremely consistent change trend between the single year of 2004 and the average after 16 years. Similar conditions were found in other years, which will not be explored here. MOYD, ATF, and NSTF reduced cloud cover by an average of 8.19% (compared to the original Terra data) or 9.72% (compared to the original Aqua data), 21.4%, and 16.43%, respectively.

Figure 8 shows an example of the procedure applied to gap filling over four days for 1 January 2001, 1 January 2005, 1 January 2009 and 1 January 2013. For these four days, Terra and Aqua maps had a cloud coverage of 85.6, 86.92, 70, and 47.58% and 84.72, 91.87, 76.1, and 61.43%, respectively. MOYD left out a cloud coverage of 84.94, 83.09, 62.73, and 41.85%. The intervention of ATF affects 25.39, 70.87, 45.25, and 29.78% of the cloud coverage on those days, respectively. NSTF removes the remaining 59.55, 12.22, 17.48, and 12.07% of the clouds on those days, respectively, and closes the gap-filling procedure. By consecutively performing three gap-filling steps, even if more than 90% of the ground information is obscured by clouds, extensive cloud gaps can be reconstructed and more reasonable ground information is obtained. This agrees with the general knowledge of snow spatial distribution characteristics in Northern XinJiang.

4.3. Accuracy Validation Results

4.3.1. Validation Based on “Cloud Assumption”

1. Accuracy: An Example of One-Day Cloud Mask

We take the special case of January 2004 as an example. The clearest MODIS Terra SCF image appeared on 21 January with 18.18% cloud gaps. This was regarded as the “cloud-free” MODIS SCF image (

Table 3. The cloud coverage of selected “cloud-free” MODIS SCF images, three extracted mask images, three artificially cloud-filled images, and after the execution of each gap-filling step for the one-day cloud mask.

	Cloud Coverage (%)
	P25	P50	P75
Original terra	18.18
Mask image	47.55	55.08	69.00
Cloud-filled image	53.59	59.78	71.87
MOYD	38.84	39.93	54.56
ATF	19.89	14.76	22.40
NSTF	0	0	0

Original Terra represents the MODIS Terra SCF image with the lowest cloud cover in January 2004; Mask image represents three extracted images closest to P25, P50, and P75 in January, respectively; Cloud-filled image represents mask images were artificially assigned to the original terra.

). The ECDF of cloud coverage in each month for the MODIS Terra SCF product in 2004 is shown in Figure 9. The closest images to P25, P50, and P75 in January were the images from 1 January 2004, 20 January 2004, and 8 January 2004, with cloud coverages of 47.55, 55.08 and 69%, respectively. After these three extracted images were assigned to the clearest SCF image on 21 January, the percentage of cloud cover in the three artificially cloud-filled images were 53.59, 59.78 and 71.87%, respectively. Note that we also performed this cloud-filled procedure for the corresponding MODIS Aqua SCF image.

After the gap-filling procedure was finished, the ground information was perfectly reconstructed visually. Despite the clouds obscuring more than half of the study area, the SCF values of the reconstructed image were extremely close to the true value of the “cloud-free” MODIS SCF image in all three cloud patterns (Figure 10). MOYD reduced 14.75, 19.85, and 17.31% of the cloud cover and ATF removed 33.7, 45.02, and 49.47% of the cloud cover for P25, P50, and P75 cloud-masked conditions, respectively. NSTF eliminates the remaining clouds and completely cloudless images were produced. As seen from

Table 4. Quantitative performance evaluation results for the one-day cloud mask.

	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$			${\mathit{R}}^2$			$\mathit{O}\mathit{E}(\%)$			$\mathit{U}\mathit{E}(\%)$			$\mathit{S}\mathit{P}\mathit{A}\mathit{E}\mathit{F}$
	P25	P50	P75	P25	P50	P75	P25	P50	P75	P25	P50	P75	P25	P50	P75
MOYD	0.07	0.11	0.13	0.95	0.89	0.84	0.37	0.89	1.46	0.36	0.89	1.48	0.89	0.86	0.82
ATF	0.07	0.13	0.14	0.96	0.86	0.85	0.38	1.99	2.57	0.27	0.67	0.77	0.66	0.66	0.72
NSTF	0.08	0.11	0.10	0.95	0.90	0.91	1.2	2.36	1.66	0.28	0.45	0.51	0.70	0.76	0.83
Final images	0.07	0.12	0.13	0.96	0.88	0.86	0.59	1.65	2.07	0.3	0.72	0.93	0.91	0.75	0.80

, the proposed gap-filling approach can achieve the preferred accuracy. The $RMSE$ was 0.07–0.13, $OE$ was 0.59%–2.07%, $UE$ was below 1%, and $R^2$ was 0.86–0.96. The spatial pattern validation metric of $SPAEF$ for the MODY, ATF, and NSTF steps, and final images exceeded 0.8, 0.65, 0.7, and 0.75, respectively. This illustrates that the gap-filled fictitious MODIS SCF image is very consistent with assumed ground truth MODIS SCF image. As expected, there was better gap-filling performance under lower cloud coverage level conditions. In addition, the performance of MOYD was superior to ATF, while ATF was superior to NSTF in most cases. This was because the implementation of the second step was based on the results of the previous step. As such, the imprecision of the previous step was passed through the input information of the second step and error was accumulated.

2. Accuracy: An Example of Multi-Day Cloud Mask

In view of the fact that the reconstruction of SCF information under consecutive cloud cover in the snow season has more significant importance for practical applications, we take the special case of the 2015 to 2016 snow season as an example. Three consecutive days characterized by lowest average cloud obstruction for MODIS Terra SCF images on February 26, February 27, and February 28 of 2016 are selected as the “cloud-free” true SCF images. The cloud coverages during these three days were 17.50, 13.12, and 64.68%, respectively (

Table 5. Cloud coverage of selected “cloud-free” MODIS Terra SCF images, three extracted mask images, three artificially filled images, and after the implementation of each gap-filling step for a consecutive three-days cloud mask. D1–D3 refers to the selected first to third days.

	Cloud Coverage (%)
	D1	D2	D3
Original Terra	17.50	13.12	64.68
Mask image	66.03	44.06	45.50
Filled image	70.43	50.44	79.59
MOYD	61.63	33.52	62.96
ATF	26.92	21.94	27.37
NSTF	0.9	1	0.95

). The corresponding images on the same dates in 2015 had cloud coverages of 66.03, 44.06, and 45.50% were used to mask the selected “cloud-free” images. The three artificially cloud-filled images were produced with cloud coverages of 70.43, 50.44, and 79.59%. The same was also done for the corresponding Aqua images.

From

Table 6. Quantitative performance evaluation results for a consecutive three-day cloud mask.

	$\mathit{R}\mathit{M}\mathit{S}\mathit{E}$			${\mathit{R}}^2$			$\mathit{O}\mathit{E}(\%)$			$\mathit{U}\mathit{E}(\%)$			$\mathit{S}\mathit{P}\mathit{A}\mathit{E}\mathit{F}$
	D1	D2	D3	D1	D2	D3	D1	D2	D3	D1	D2	D3	D1	D2	D3
MOYD	0	0.06	0.19	1	0.97	0.93	0	0.15	0.35	0	0.34	7.82	0.94	0.90	0.95
ATF	0.07	0.13	0.08	0.97	0.93	0.95	0.74	0.85	0.47	0.15	1.80	1.43	0.90	0.97	0.84
NSTF	0.15	0.13	0.12	0.90	0.88	0.93	4.39	2.76	2.92	0.23	1.12	1.57	0.74	0.73	0.70
Final images	0.09	0.10	0.11	0.96	0.93	0.94	1.58	1.15	1.27	0.15	0.93	2.41	0.87	0.88	0.80

, we can see that the performance of this consecutive cloud cover case is slightly worse, compared to the one-day cloud mask case. MOYD reduced the cloud cover by 8.8, 16.92, and 16.73% on the three days, respectively. ATF reduced the cloud cover by 34.71, 11.58, and 35.59% on the three days, respectively. Although the NSTF step eliminated a lot of cloud cover (26.02, 20.94, and 26.42% on the three days, respectively), it does not completely eliminate the cloud cover. Tiny amounts (i.e., 0.9, 1, and 0.95% on the three days, respectively) of cloud gaps were still retained. Higher accuracy was also achieved with lower overall $RMSE$ values of 0.09–0.11, $OE$ values of 1.15–1.58%, $UE$ values of 0.15–2.41%, higher $R^2$ values of 0.93–0.96, and higher $SPAEF$ values of 0.8–0.88. Despite the significant differences in terms of both cloud percentage and cloud patterns, there was little performance difference among these three days.

3. Validation Results Based on “Cloud Assumption”

The generalization validation results of 576 one-day and 90 groups of multi-day (2–7 days) cloud masked images are shown in Figure 11. The comparison of the result of each gap-filing step with the selected “cloud-free” true MODIS SCF image confirms the efficacy of the proposed gap-filling method. For the one-day case, most cloud-free images on each month can be found with an average cloud coverage less than 1%. The average cloud coverage of the artificially cloud-filled images was 30, 47, and 66%, respectively, when P25, P50 and P75 additional cloud cover was introduced. The average cloud coverage of the selected 2–7 consecutive “cloud-free” days was approximately 11–27% and 50–64% after artificially spatiotemporal additional cloud cover was incorporated, as shown in Figure 11a. The gap-filling efficiency is further confirmed by Figure 11b. More than 99% of gaps were reconstructed. The average amount of cloud gaps reduced by MODY, ATF, and NSTF are 9.21, 25.19, and 17.28%, which are basically similar with the results averaged from 16 years given in Section 4.2. Based on Figure 11c–f and Figure 12, although there are significant differences in both cloud percentage and cloud patterns among different validation strategies, all of them are characterized by lower values of $RMSE$, $OE$, and $UE$, as well as higher values of $R^2$ and $SPAEF$. The average $RMSE$ was about 0.1, $R^2$ was greater than 0.8, $OE$ was 1.13%, $UE$ was 1.4%, and $SPAEF$ was 0.78 for the final images. The performance of MOYD and ATF was very prominent, with the $SPAEF$ below 0.1, the $OE$ and $UE$ below 1%, and $SPAEF$ exceeding 0.8 in all validation strategies. However, the performance of NSTF was slightly worse, particularly for the validation strategy of four to seven consecutive cloud covered days. The $RMSE$ reached up to 0.15, some disagreement errors exceeded 2%, and the average $SPAEF$ was approximately 0.65. This is because, after three days of backward temporal filtering in ATF, it will inevitably leave out a larger percentage of gaps. There will also be less available non-cloud pixels for NSTF. Together with the accumulative effect of error, these factors lead to a slightly inferior performance of NSTF.

4.3.2. Validation Based on In-Situ SD Observations

Four validation indices were calculated by using an SD threshold value of ${\epsilon }_1=3$ and MODIS SCF threshold value of ${\epsilon }_2=50\%$. We identified a pixel as snow covered when the SD or SCF exceeds the threshold value, otherwise it was classified as a land pixel. As shown in Figure 13, we found that, for each elevation zones, the performance of the MODIS Terra SCF product was superior to MODIS Aqua’s. The MODIS Terra SCF product had a slightly higher $\mathrm{OA}$, lower $\mathrm{MU}$ and $\mathrm{MO}$, and a $\mathrm{Bias}$ closer to 1. This may be attributed to the MODIS Aqua SCF mapping algorithm that uses band 7 instead of band 6 and the changed underlying surface conditions [50]. The performance of each gap-filling step is varied. For $\mathrm{OA}$, the average value increased from 50.86% (Terra)/48.55% (Aqua) to 58.15% after MOYD, to 78% after ATF, and 93.72% after NSTF. The improvement of $\mathrm{OA}$ is more significant in elevation zones of <2000 m. The $\mathrm{OA}$ remains almost unchanged in the zones of elevation above 2000 m. This may be because only three in situ observation stations were located in this area, and it may not be representative enough. For $\mathrm{Bias}$, three consecutive gap-filling steps bring it closer to 1. This proves that the probability of MODIS observed snow events is increasing. The $\mathrm{MU}$ and $\mathrm{MO}$ errors for the study area for 16 years were 2.16 and 2.43%, respectively, for the MODIS SCF product under clear-sky conditions, and they were 3.32 and 2.96%, respectively, for the final, gap-filled SCF product.

5. Discussion

5.1. The Sensitivity of Different SCF and SD Threshold Combinations

When comparing the MODIS SCF value with in situ SD observations, the most suitable SD threshold value (${\epsilon }_1$) and SCF threshold value (${\epsilon }_2$) must be chosen. ${\epsilon }_1=1$ [26,51,52], ${\epsilon }_1=3$ [21,28], ${\epsilon }_1=4$ [22,24,30,31], ${\epsilon }_2=15\%$ [53,54], and ${\epsilon }_2=50\%$ [12,55] have been used in the literature. We use these thresholds to analyze the sensitivity of different thresholds to validate results, as shown in Figure 14.

The results indicated that a higher $\mathrm{OA}$ (above 90%) and lower $\mathrm{MU}$ and $\mathrm{MO}$ (below 2%) are achieved for all the different threshold combinations in the non-snow season, but $\mathrm{Bias}$ is closest to 1 when the threshold value combination of ${\epsilon }_1=3$ and ${\epsilon }_2=50\%$ is used. This illustrates that the inconsistency between the final gap-filled MODIS SCF product and in situ SD observations is minimal. In snow seasons, $\mathrm{OA}$ is over 80% and $\mathrm{Bias}$ is also extremely close to 1 for all different threshold combinations, which further confirms that there is favorable consistency between the final gap-filled MODIS SCF product and in situ SD observations. The MODIS SCF threshold value of ${\epsilon }_2=50\%$ produced a relatively a higher $\mathrm{MU}$ and lower $\mathrm{MO}$. When the ${\epsilon }_2$ is fixed, a large value of ${\epsilon }_1$ produces relatively lower $\mathrm{MU}$ and higher $\mathrm{MO}$ values. For example, $\mathrm{MU}$ is 12% and $\mathrm{MO}$ is 0 for the threshold combination of ${\epsilon }_1=1$ and ${\epsilon }_2=50\%$, and $\mathrm{MU}$ is 8.2% and $\mathrm{MO}$ is 9.08% for the threshold combination of ${\epsilon }_1=4$ and ${\epsilon }_2=50\%$ in March. The above results demonstrate that the most suitable ${\epsilon }_1$ should be between 1 and 4cm and ${\epsilon }_2$ should be between 15 and 50%. Therefore, we drew the inference that the threshold combination of ${\epsilon }_1=3$ and ${\epsilon }_2=50\%$ is the best choice for the MODIS accuracy assessment in the study area.

5.2. Comparison with Another Gap-Filling Method

We compared the proposed method with the existing cubic spline temporal interpolation approach. Figure 15 is the monthly/annual variation of the four validation indices. Although the cubic spline interpolation method uses information from both preceding and following days and the proposed method only uses information from preceding days, it was found that the proposed method performs more excellently, particularly during the snow season. $\mathrm{Bias}$ was closer to 1, and the average $\mathrm{OA}$ was 87%, $\mathrm{MU}$ was 6.85%, and $\mathrm{MO}$ was 6.3%. These values are 2% higher, 1.92% lower, and 0.04% lower, respectively than the corresponding values for the cubic spline interpolation method. It must be noted that the $\mathrm{Bias}$ for the period from June to September is far more than 1 for both gap-filling methods. The value even exceeds 3 for the spline interpolation approach. This means that both methods have poor consistency with in situ SD observations during this period. On the other hand, the $\mathrm{OA}$ values of both gap-filling methods are very high. Their average values even reach 99%. This seems to be contradictor. The cause of this phenomenon is that there are only a few permanent snow covers in the summer and some slight overestimations occurred.

5.3. Dependence of Gap-Filling Performance on Land Cover Types

The spatio-temporal distribution and “accumulation-melting” pattern of snow cover is closely related to land cover types [8,10,45]. Thus, we examined the performance of each gap-filling step among different land cover types, as shown in

Table 7. The effect of the land cover types on the four validation indices of $\mathrm{OA}$, $\mathrm{Bias}$, $\mathrm{MU}$, and $\mathrm{MO}$.

Land Cover Type	SCF Images	L-C	L-L	L-S	S-C	S-L	S-S	Cloud (%)	$\mathbf{OA}(\%)$	$\mathbf{Bias}$	$\mathbf{MU}(\%)$	$\mathbf{MO}(\%)$
Cropland (22)	MOD	33,854	53,865	940	20,831	249	10,183	45.60	53.41	0.36	0.38	1.44
	MYD	35,813	52,180	666	20,843	227	10,193	47.24	52.01	0.35	0.36	1.05
	MOYD	27,589	59,971	1099	17,494	328	13,441	37.59	61.22	0.47	0.44	1.47
	ATF	9022	77,985	1652	10,715	638	19,910	16.46	81.63	0.69	0.64	1.65
	NSTF	0	85,526	3133	0	1378	29,885	0	96.24	1.05	1.15	2.61
Grasslands and shrubs (11)	MOD	18,499	23,806	1421	10,666	1274	4295	48.64	46.87	0.35	4.14	4.61
	MYD	19,929	22,974	823	12,034	1424	2777	53.31	42.95	0.22	5.09	2.95
	MOYD	15,328	26,790	1608	9564	1715	4956	41.51	52.94	0.4	4.89	4.59
	ATF	5677	35,727	2322	6186	2664	7385	19.78	71.9	0.6	5.54	4.83
	NSTF	0	40,216	3510	0	4236	11,999	0	87.08	0.96	6.96	5.85
Urban and built up (9)	MOD	14,317	21,174	75	10,004	142	3347	49.58	49.98	0.25	0.57	0.3
	MYD	15,102	20,434	30	11,038	248	2207	53.28	46.15	0.16	1.08	0.13
	MOYD	11,708	23,765	93	9134	276	4083	42.48	56.76	0.31	0.98	0.33
	ATF	4058	31,351	157	6688	506	6299	21.90	76.74	0.48	1.32	0.41
	NSTF	0	34,914	652	0	1396	12,097	0	95.83	0.95	2.85	1.33
Forests (3)	MOD	3861	7171	25	3713	134	1449	46.32	52.71	0.28	1.53	0.28
	MYD	4118	6923	16	3831	168	1297	48.61	50.27	0.25	2	0.19
	MOYD	3070	7954	33	3232	193	1871	38.54	60.08	0.36	1.92	0.33
	ATF	738	10,267	52	2093	335	2868	17.31	80.32	0.55	2.48	0.38
	NSTF	0	10,929	128	0	524	4772	0	96.01	0.93	3.2	0.78
Unused (4)	MOD	6183	9470	119	4069	515	1448	47.02	50.07	0.26	4.46	1.03
	MYD	6511	9181	80	3938	517	1577	47.92	49.34	0.27	4.55	0.7
	MOYD	4882	10,751	139	3443	632	1957	38.18	58.26	0.35	4.69	1.03
	ATF	1654	13,923	195	2315	919	2798	18.20	76.69	0.5	5.15	1.09
	NSTF	0	15,268	504	0	1421	4611	0	91.17	0.85	6.52	2.31

Note: The figures in parentheses represent the number of in situ SD observation stations. S, L, and C represent snow, land, and cloud pixels respectively. the first letter indicates in situ SD observation and the second letter represents MODIS observation.

. We found that the $\mathrm{OA}$ of original MODIS SCF products is approximately 50%. The highest accuracy appears in cropland, and relatively low accuracies occur in grasslands and shrubs. The disagreement between $\mathrm{MU}$ and $\mathrm{MO}$ errors occurred in the various land cover types, particularly in grassland and unused land. The performance of $\mathrm{OA}$ and $\mathrm{Bias}$ improved a lot, whereas there was no significant increase in $\mathrm{MU}$ and $\mathrm{MO}$ error. The significant improvements were found in cropland, forests, and urban and built up land, with $\mathrm{OA}$ values exceeding 95%. The improvements of grasslands and shrubs, and unused land are also obvious, with $\mathrm{OA}$ values of 87.08 and 91.17%, respectively. Note that, only 49 observation stations were used in our analysis because another observation station (Daxigou) was located at the eastern edge of the Glacier No. 1 in Tianshan, may be due to the errors of geolocation or the land cover type data, the corresponding land cover type of this observation station was permanent glacier, but we found that the cumulative number of snow covered days in 16 years at this station was 1517 days, this means that this observation station was not sufficiently representative as it is generally believed that the annual average number of snow days in permanent snow areas is more than 320 days, we therefore give it up. In addition, the final result is related to the number and location of the observation sites, it may not be enough to explain the essence of the problem. The various land cover type may lead to different spatio-temporal distribution and “accumulation-melting” patterns of snow cover. Varying intensities of illumination of the heterogeneous terrain in complex mountainous areas, wind-drifting snow, avalanche, and snow cover left in garden and streets caused by man-made operations, all these factors make the problem more complicated, we will explore the dependence of gap-filling performance on land cover types by introduce higher resolution data from Landsat OLI 8, Sentinel-2, geostationary satellites in the future.

6. Conclusions

Due to atmospheric disturbances, much ground information usually cannot be acquired. This is what mainly limit MODIS SCF products from being widely applied. Addressing this challenge, an innovative gap-filling method for MODIS SCF products based on the concept of NSTF is proposed. Ground information of a gap pixel is estimated through selected similar pixels in the remaining known region via an automatic machine learning technique of NAGA-II. As there are significant gaps in large cloud regions, it is unreliable and unrealistic to rely on the remaining known regions of the image to match similar pixels. The widely confirmed method of MOYD and ATF are first executed sequentially to reconstruct some gap pixels. However, the NSTF method can be used independently.

In this study, we first answered the two questions on how much cloud gaps are in the MODIS SCF products and how long is the cloud cover duration. The influence of topography on cloud cover is also evident. Almost cloud-free daily MODIS SCF images can be produced after consecutive implementation of three cloud removal steps. The long-term mean cloud coverage is 0.52% on average. The generalization validation results of 576 one-day cloud masked images and 90 groups of multi-day (2–7 days) spatiotemporal additional cloud masked images prove the reliability and feasibility of the proposed methodology. The average $RMSE$ is approximately 0.1, $R^2$ exceeds 0.8, $OE$ error is 1.13%, $UE$ error is 1.4%, and $SPAEF$ is 0.78, even if more than 80% of the surface covering information is missing. The validation based on 50 in situ SD observations also demonstrates the superiority of the method in terms of accuracy and consistency. The average $\mathrm{MU}$ and $\mathrm{MO}$ error have increased about 1.16 and 0.53% for the study area in 16 years. In addition, the validation test in terms of land cover types also proves the efficiency and accuracy of the proposed gap-filling method. In general, existing MODIS cloud gap-filling methods have satisfactory accuracy in the stable snow cover phase, but they cannot achieve a satisfying accuracy in the transition phase of snow cover accumulation and ablation [10,14]. Different from the existing methods, the proposed method has a no significant performance difference in seasonal or different elevation zones. NSTF can consider both the spatio-temporal correlation of multitemporal MODIS SCF products and auxiliary information, such as topography and land cover, in an integrated manner, instead of analyzing each individually. Thus, it can more fully exploit the space-time correlation within a pixel. Nevertheless, there are still some limitations that should be improved. For example, we simply used a fixed value of 10 as the number of similar pixels. This is obviously defective. We will use adaptive variable values according to the cloud coverage level in subsequent studies.