Validation of Cloud-Gap-Filled Snow Cover of MODIS Daily Cloud-Free Snow Cover Products on the Qinghai–Tibetan Plateau

Abstract

Accurate daily snow cover extent is a significant input for hydrological applications in the Qinghai–Tibetan Plateau (QTP). Although several Moderate Resolution Imaging Spectroradiometer (MODIS) daily cloud-free snow cover products over the QTP are openly accessible, the cloud-gap-filled snow cover from these products has not yet been validated. This study assessed the accuracy of cloud-gap-filled snow cover from three open accessible MODIS daily products based on snow maps retrieved from Landsat TM images. The F1-score (FS) from daily cloud-free MODIS snow cover for the combined MOD10A1F and MYD10A1F (SC1) was 64.4%, which was 7.4% points and 5.3% points higher than the other two commonly used products (SC2 and SC3), respectively. The superior accuracies from SC1 were more evident in regions with altitudes lower than 5000 m, with a weighted average FS by the area percentage of the altitude regions of 58.3%, which was 6.9% points and 9.1% points higher than SC2 and SC3. The improved SC1 accuracies also indicated regional clustering characteristics with higher FS values compared to SC2 and SC3. The lower accuracies of cloud-gap-filled snow cover from SC2 and SC3 were mainly due to the limitation in determining snow cover based on the method of the inferred snow line and the overestimation of the Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) snow water equivalent (SWE). These results indicate that the temporal filter approach used in SC1 is a good solution to produce daily cloud-gap-filled snow cover data for the QTP because of its higher accuracy and simple computation. The findings can be helpful for the selection of cloud-removal algorithms for determining snow cover dynamics and phenological parameters on the QTP.

1. Introduction

The Qinghai–Tibetan Plateau (QTP) is known as the earth’s ‘‘third pole’’ and the Asian Water Tower [1,2,3,4]. The QTP is the source region of major rivers, such as the Yangtze River, Yellow River, and Lancang River [5,6], and numerous lakes are distributed in this region. Water from these river basins fulfills the needs of drinking water, irrigation, energy, industry, and sanitation for approximately 210 million people within the region and approximately 1.3 billion people downstream of the region [7,8,9,10,11]. Snow cover dynamics on the QTP influence water availability downstream in the major river basins of Asia [12], which plays a significant role in the hydrological cycle [13,14,15,16]. Snow cover phenological parameters have become important indices that indicate the response of the hydrological cycle to climatic change [17]. Accurate daily snow cover extent on the QTP is a significant input for hydrological applications [18].

Since 2001, daily global snow cover products from MODIS (Moderate Resolution Imaging Spectroradiometer) data with a spatial resolution of 500 m are available [19,20] and have become one of the most popular data resources for snow-related studies [12,21,22,23,24,25,26,27]. However, MODIS snow cover products are retrieved from optical bands, and a significant portion of the daily snow cover data is missing because clouds obscure the area in MODIS images. To remove the data gaps caused by cloud cover, a series of methods has been developed to rebuild cloud-free MODIS daily snow products [28,29,30,31,32,33,34,35,36,37]. Presently, several long-term MODIS daily cloud-free snow cover products over the QTP based on these methods are openly accessible, including MOD10A1F (global Level-3 dataset of daily cloud-free snow cover derived from the MODIS/Terra Snow Cover Daily L3 Global 500m SIN Grid dataset (MOD10A1), https://nsidc.org/data/mod10a1f/versions/61, accessed on 2 November 2022) and MYD10A1F (global Level-3 dataset of daily cloud-free snow cover derived from the MODIS/Aqua Snow Cover Daily L3 Global 500m SIN Grid dataset (MYD10A1), https://nsidc.org/data/myd10a1f/versions/61, accessed on 8 November 2022) [36] from the National Snow and Ice Data Center (https://nsidc.org/home, accessed on 8 November 2022), and daily cloud-free snow cover products conducted by [28,29]. However, the validation of these commonly used products for QTP is limited.

Except for the validation conducted by Huang et al. [28] and Qiu et al. [29], only Hao et al. [38] compared the accuracies of different daily cloud-free snow products over the QTP, including the products developed by Huang et al. [28] and Qiu et al. [29]. Hao et al. [38] assessed the accuracies of TAI (Terra–Aqua–IMS), IMS (Interactive Multi-sensor Snow products), MOD-SSM/I (combination of the MODIS and Special Sensor Microwave Image snow products, similar to the product conducted by Huang et al. [28], except that the AMSR-E data were replaced by SSM/I snow products), and MOD-B (blending method based on MODIS snow cover products, similar to the product conducted by Qiu et al. [29]). However, the validation results differed from those reported by Huang et al. [28] and Qiu et al. [29]. The assessment of MOD10A1F on the QTP has not been conducted to date. Although Hao et al. [39] assessed MOD10A1F over China during the snow season from 2000 to 2020, the study did not provide a detailed validation of the QTP.

All these previous validations were carried out under all-weather conditions, implying that the validation results were obtained under clear sky and cloud conditions. Thus, the accuracies of the cloud-gap-filled snow cover from these products cannot be evaluated in this manner as it is biased by data from clear sky conditions. As a result, it is difficult for users to determine the efficiency of cloud removal algorithms. In addition, existing validations were conducted based on in situ observation data [28,29,38,39,40]. Owing to the harsh natural environment, meteorological and hydrological observation stations on the QTP are scarce, and their distribution is unbalanced. The limited stations are only located in the valleys with lower elevations (Figure 1). The sparse distribution of meteorological stations can directly induce significant uncertainties in the validation of snow cover. In addition, the snow cover in the low-altitude regions is less than in the higher-altitude regions. The spatial bias of meteorological stations cannot clarify the regional differences in the accuracies of snow cover. Thus, large uncertainties may exist in the validation of snow cover products [5,18,41].

Therefore, the objective of this study was to assess the accuracy of the cloud-gap-filled snow cover of three open accessible daily cloud-free snow cover products derived from MODIS data. To decrease the impact of the spatial bias of meteorological stations, snow maps retrieved from Landsat TM images were used as validation data in this study because of their high spatial resolution [20,27,42]. This study addressed three questions: (1) Which product achieves the highest precision for cloud-gap-filled snow cover? (2) Are there spatial differences in precision between these products? (3) What are the major limitations of the cloud removal algorithms used in these products?

2. Study Area and Data

2.1. Study Area

The QTP is located in the central Eurasian continent (Figure 1) and covers an area of approximately 2.6 × 10⁶ km² with an average elevation of over 4000 m, the largest and highest elevation plateau in the world [43,44,45]. The QTP ranges from approximately 26°00′12″N to 39°46′50″N from north to south and from 73°18′52″E to 104°46′59″E from west to east [46]. Over the entire QTP, the annual mean temperature is below 0 ℃ and indicates spatial heterogeneity, which generally decreases from east to west [44]. The total annual precipitation varies from approximately 1500 mm in the southeast to <100 mm in the northwestern parts of the QTP [45]. Approximately 90% of the total annual precipitation occurs in the southeastern QTP during the wet season (March to August) [47]. The main vegetation type in the QTP is grassland. Alpine meadows, steppes, and desert steppes cover 63% of the QTP [48] and are distributed from the southeast to the northwest.

2.2. Data and Data Preprocessing

2.2.1. MODIS Daily Cloud-Free Snow Cover Products

MOD10A1F and MYD10A1F products for eight hydrological years from 1 September 2002 to 31 August 2010 were downloaded from the National Snow and Ice Data Center (https://nsidc.org/data/, accessed on 8 November 2022). Snow cover from MOD10A1F and MYD10A1F products was determined by a normalized difference snow index (NDSI) of 0.4 as a threshold based on the method proposed by Hall et al. [36]. In addition, the combined product MOYD10A1F (hereafter labeled as SC1) was conducted based on both MOD10A1F and MYD10A1F products (if clouds cover MOD10A1F and the corresponding grid in MYD10A1F is under clear-sky conditions, then the grid code in MOD10A1F is replaced by MYD10A1F).

Daily cloud-free snow cover products developed by Huang et al. [28] (hereafter labeled as SC2) and Qiu et al. [29] (hereafter labeled as SC3) during the same period were downloaded from the Qinghai–Tibetan Plateau Data Center (http://data.tpdc.ac.cn/zh-hans/data, accessed on 8 November 2022). SC2 and SC3 used MOD10A1 and MYD10A1 to establish cloud-free snow-cover products. The NDSI threshold values adopted in SC2 and SC3 are 0.4 and 0.36, respectively. MODIS snow cover data used in SC2 were Version 5, and in SC1 and SC3, they were Version 6.

The cloud condition for a pixel was determined based on the “cloud-persistence count” (CPC) map in MOYD10A1F for all three daily cloud-free snow cover products [36]. The CPC represented the number of consecutive days of cloud obscuration since the latest cloud-free day of the pixel. Pixels with a CPC value of zero were cloud-free ones, otherwise not. Thus, pixels with values greater than zero in the CPC map were considered under the cloud.

2.2.2. Landsat TM

A total of 4868 Landsat TM surface reflectance images between 2002 and 2010 with a cloud cover of less than 10% were obtained from the United States Geological Survey (USGS) website (https://espa.cr.usgs.gov/index/, accessed on 8 November 2022). Data gaps (no-data) in these images were removed. As the cloud cover area was low in the used images and we only compared the accuracies of the cloud-gap-filled snow cover of three open accessible MODIS daily products, the influence of the cloud pixels on the validation results should be limited. Thus, the cloud pixels on the Landsat TM images were not further processed. These images cover almost the entire QTP, except for a small area in the South (Figure 2). For the same place (the same path and row of TM), about 32 TM images were used on average in this study.

TM snow cover maps were derived from these images using the snow-mapping algorithm (SNOMAP) [49]. Landsat TM bands 2 and 5 were used to calculate the NDSI, and band 4 was used to remove the effect of water bodies. The pixel with an NDSI greater than 0.4 and band 4 reflectance greater than 0.11 was classified as snow cover [27,42]; alternatively, it was classified as non-snow cover. The Landsat TM snow cover maps were then resampled to 500 m based on the abundance of snow cover or land area within the window to match the resolution of the MODIS daily snow cover products [27,42].

2.2.3. Other Auxiliary Data

Digital elevation model (DEM) data with a spatial resolution of 90 m were downloaded from the Shuttle Radar Topography Mission (SRTM) database (http://srtm.csi.cgiar.org/, accessed on 8 November 2022). Lake and river data downloaded from the National Catalogue Service for Geographic Information (https://www.webmap.cn/, accessed on 8 November 2022) were used to mask the water bodies over the QTP. Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) snow water equivalent (SWE) products were obtained from the National Snow and Ice Data Center (https://nsidc.org/data/amsre/, accessed on 8 November 2022) to clarify the cloud removal method used in SC2, where AMSR-E was used to determine whether cloud pixels in the patch snow area were snow cover or not.

3. Methodology

Not all clouds stay in the same place for a whole day and quite a lot of them move quickly [50]. The receiving times of Landsat and MODIS images (Aqua and Terra) are different on the same day. Thus, sometimes, there is no cloud on Landsat TM images while there is on MODIS ones at the same place due to their time lag. In this situation, the cloud-free TM images can be used as the validation data for cloud-gap-filled snow cover when we assume that the snow cover change is small within the short time gap that is less than one day.

Commonly used accuracy assessment indices were used for validation, including the balanced accuracy (BA), F1-score (FS), precision (PC), and recall (RC) [38,40,51,52]. BA is the balanced accuracy that accounts for the unbalance class distributions. It can overcome the drawbacks of the commonly used overall accuracy, which may lead to false conclusions for an imbalanced dataset [52]. PC shows omission errors with a higher PC, indicating lower omission errors. RC shows commission ones with a higher RC, indicating lower commission ones. FS shows overall accuracy considering the balance of precision and recall. For each Landsat TM snow map, the PC, RC, FS, and BA were calculated based on the confusion matrix as follows:

$${\mathrm{PC}}_{\mathrm{single}}^{\mathrm{sc}}=\frac{TP}{TP+FP}$$

(1)

$${\mathrm{RC}}_{\mathrm{single}}^{\mathrm{sc}}=\frac{TP}{TP+FN}$$

(2)

$${\mathrm{FS}}_{\mathrm{single}}^{\mathrm{sc}}=\frac{2\times {\mathrm{PC}}_{\mathrm{single}}^{\mathrm{sc}}\times {\mathrm{RC}}_{\mathrm{single}}^{\mathrm{sc}}}{{\mathrm{PC}}_{\mathrm{single}}^{\mathrm{sc}}+{\mathrm{RC}}_{\mathrm{single}}^{\mathrm{sc}}}$$

(3)

$${\mathrm{PC}}_{\mathrm{single}}^{\mathrm{nsc}}=\frac{TN}{TN+FN}$$

(4)

$${\mathrm{RC}}_{\mathrm{single}}^{\mathrm{nsc}}=\frac{TN}{TN+FP}$$

(5)

$${\mathrm{FS}}_{\mathrm{single}}^{\mathrm{nsc}}=\frac{2\times {\mathrm{PC}}_{\mathrm{single}}^{\mathrm{nsc}}\times {\mathrm{RC}}_{\mathrm{single}}^{\mathrm{nsc}}}{{\mathrm{PC}}_{\mathrm{single}}^{\mathrm{nsc}}+{\mathrm{RC}}_{\mathrm{single}}^{\mathrm{nsc}}}$$

(6)

$${\mathrm{BA}}_{\mathrm{single}}=\frac{1}{2}(\frac{TP}{TP+FP}+\frac{TN}{TN+FN})$$

(7)

${\mathrm{PC}}_{\mathrm{single}}^{\mathrm{sc}}$, ${\mathrm{RC}}_{\mathrm{single}}^{\mathrm{sc}}$, ${\mathrm{FS}}_{\mathrm{single}}^{\mathrm{sc}}$, ${\mathrm{PC}}_{\mathrm{single}}^{\mathrm{nsc}}$, ${\mathrm{RC}}_{\mathrm{single}}^{\mathrm{nsc}}$, ${\mathrm{FS}}_{\mathrm{single}}^{\mathrm{nsc}}$, and ${\mathrm{BA}}_{\mathrm{single}}$ are indices for one Landsat TM snow map. The superscripts ‘sc’ and ‘nsc’ indicate the index for snow cover or non-snow cover, respectively. TP, FP, FN, and TN are true positive, false positive, false negative, and true negative numbers, respectively (

Table 1. The confusion matrix for MODIS snow cover product validation.

		Reference
		Snow	Non-Snow
Map	Snow	TP	FP
	Non-snow	FN	TN

).

The distribution of snow cover varies significantly over the study area [31]. Therefore, if snow cover in the whole study area was directly considered for validation, the Landsat TM snow maps with extremely high or low percentages of snow cover pixels may lead to a biased overall accuracy over the entire study area. Therefore, the PC, RC, FS, and BA for the entire QTP or sub-region are the average values. For example, the average of PCs can be calculated as follows:

$${\mathrm{PC}}_{\mathrm{region}}^{\mathrm{sc}}=\frac{{{\displaystyle \sum }}_1^n{\mathrm{PC}}_{\mathrm{single}}^{\mathrm{sc}}}{n}$$

(8)

where n is the number of Landsat TM snow maps covered over the entire QTP or sub-region of the QTP. The pixel numbers of TP, FP, FN, and TN were calculated by overlay analysis using the resampled TM snow cover maps and the corresponding MODIS daily snow cover product.

Except for the assessment of the accuracies over the entire study area, the QTP was also divided into different regions based on altitude (

Table 2. The altitude regions over the whole Qinghai–Tibetan Plateau.

Region ID	Altitude Range (m)	Area Percent (%)
1	≤3000	9.9
2	3000–3500	8.1
3	3500–4000	10.4
4	4000–4500	16.7
5	4500–5000	31.2
6	5000–5500	19.8
7	˃5500	4.0

) and a 10 × 10 km grid window to clarify the spatial differences of validation indices on the QTP. First, the areas within each altitude zone and 10 × 10 km grid were trimmed for every TM snow cover map. Then, the PC, RC, FS, and BA (Equations (1)–(7)) were recalculated for every trimmed TM snow cover map. Finally, the average of these indices was computed for each altitude zone and 10 × 10 km grid.

4. Results

4.1. The Accuracies of Daily Cloud-Free Snow Cover over the Whole QTP

A total of 41,379,105 pixels under cloudy conditions and 459,464,960 pixels under clear conditions were used for validation, respectively. For the entire snow cover product under clear and cloudy skies, the accuracies of SC1 were better than those of SC2 and SC3. The BA of SC1 was 71.9%, which was 1.2% points and 0.5% points higher than SC2 and SC3. For snow cover, although the PC of SC1 was 47.4%, which was 0.6% points and 0.4% points lower than SC2 and SC3, the RC and FS of SC1 were 64.0% and 54.7%, which was 13.9% points and 6.6% points higher than the RC of SC2 and SC3, and 6.7% points and 3.3% points higher than the FS of SC2 and SC3. For the non-snow cover, the RC of these products were close to each other, but the PC and FS of SC1 were 96.4% and 95.2%, which were 3.0% points and 1.4% points higher than the PC of SC2 and SC3, and 1.5% points and 0.8% points higher than the FS of SC2 and SC3 (

Table 3. Validation results for the whole Qinghai–Tibetan Plateau. BA, FS, PC, and RC were balanced accuracy, F1-score, precision, and recall, respectively. SC1, SC2, and SC3 were the snow cover products of the combined MOD10A1F and MYD10A1F (SC1) developed by Huang et al. [28] and Qiu et al. [29], respectively.

Sky Condition for MODIS Data	Product Name	Snow Cover (%)			Non-Snow Cover (%)			BA (%)
		PC	RC	FS	PC	RC	FS
Clean and cloudy	SC1	47.4	64.0	54.7	96.4	94.4	95.2	71.9
	SC2	48.0	50.1	48.1	93.4	94.5	93.7	70.7
	SC3	47.8	57.3	51.4	95.0	94.4	94.4	71.4
Cloudy	SC1	64.8	62.1	64.4	77.9	82.6	78.7	71.4
	SC2	62.5	51.1	57.0	71.7	80.8	74.1	67.1
	SC3	60.3	56.4	59.1	78.3	80.6	78.0	69.3

).

For cloud-gap-filled snow cover, the accuracies of SC1 were also better than those of SC2 and SC3, especially for snow cover. The BA of SC1 was 71.4%, which was 4.3% points and 2.1% points higher than SC2 and SC3. The FSs of snow and non-snow cover for SC1 were 64.4% and 78.7%, which were 7.4% points and 4.6% points higher than those for SC2 and 5.3% points and 0.7% points higher than SC3 (Table 3).

For cloud-gap-filled snow cover, the PC and RC of SC1 for snow cover were higher than those from SC2 and SC3. The PC of SC1 was 64.8%, which was 2.3% points and 4.5% points higher than SC2 and SC3. The RC of SC1 was 62.1%, which was 10.0% points and 5.7% points higher than SC2 and SC3. SC1 also indicated a better balance between commission and omission errors, particularly for snow cover. The difference between PC and RC for snow cover for SC1 was 2.7%, which was 8.6% points and 1.1% points lower than those for SC2 and SC3 (Table 3).

4.2. The Accuracies of Cloud-Gap-Filled Snow Cover within Different Altitude Zones

The PC, RC, and FS for snow cover for all products increased with increasing altitude (Figure 3a–c). The RCs of SC1 changed more smoothly than those of SC3 and SC2. The RC of SC1 increased by 38.4% points, from 41.0% to 79.4%, which was less than SC2 by 54.4% points, from 22.5% to 76.9 %, and less than SC3 by 45.9% points, from 30.7% to 76.6%. The PC of SC1 and SC2 changed more smoothly than that of SC3. The PC of SC1 and SC2 increased by 43.4% points, from 41.8% to 85.2%, and by 41.2% points, from 45.1% to 86.3%, respectively, less than SC3 by 68.2% points, from 22.5% to 90.7%. The FS of SC1 also changed more smoothly than those of SC2 and SC3. FS of SC1 increased by 29.3% point, from 53.2% to 82.5%, which was less than SC2 by 39.9% points, from 42.1% to 82.0 %, and SC3 by 46.7% points, from 36.5% to 83.2%.

The better accuracy for the cloud-gap-filled snow cover of SC1 was more prominent in regions with altitudes lower than 5000 m compared with SC2 and SC3. The average FS (weighted by the area percentage of the altitude regions) of SC1 was 58.3%, which was 6.9% points and 9.1% points higher than SC2 and SC3 (Figure 3c). In regions with altitudes higher than 5000 m, the average FS of SC1 was 72.7%, which was 3.4% points and 0.6% points higher than SC2 and SC1 (Figure 3c).

The better accuracy for the cloud-gap-filled snow cover of SC1 was related to the higher PC compared with SC3 and better RC compared with SC2 and SC3. The average PC of SC1 was 58.0%, which was 7.6% points higher than SC3 (Figure 3a). Notably, the higher PC of SC1 was mainly in regions with altitudes lower than 5000m. On the contrary, in regions with altitudes higher than 5000 m, the PC of SC3 was higher than that of SC1 (Figure 3a). The average RC of SC1 was 55.7%, which was 10.4% points and 6.0% points higher than SC2 and SC3, especially at lower altitudes (Figure 3b).

SC1 showed the best balance between commission and omission errors for cloud-gap-filled snow cover in each altitude region. The difference between the PC and RC for SC1 was 0.5–6.9%, lower than SC2 of 7.8–22.6% (decreased by 38.1–96.3% points) and SC3 of 1.0–14.2% (decreased by 7.6–90.0% points). Especially in lower altitude regions (<5000 m), the average differences between PC and RC of SC1, SC2, and SC3 were 1.6%, 12.3%, and 3.5%, respectively, decreased by 85.0% and 54.3% points compared with SC2 and SC3.

4.3. Accuracy of Cloud-Gap-Filled Snow Cover Based on the 10 × 10 km Grid

The higher accuracy of SC1, compared with SC2 and SC3, indicates regional clustering characteristics. There were two regions in SC2 (red frame 1, 2 in Figure 4a and Figure S1a) where the FS was evidently lower than SC1 (Figure S1c). The first region is located between the Qaidam Basin and Tanggula Mountains (red frame 1), and the other region is located in the Tanggula and Gangdese mountains (red frame 2). For SC3, there were also two regions in which the FS was lower than that of SC1 ((red frame 3, 4 in Figure 4b and Figure S1b)). One region was located in the northwest of the Qaidam Basin (red frame 3), and the other was in the northeast of the QTP around Qinghai Lake (red frame 4).

For PC, there were no evident spatial differences between SC1 and SC2 (Figure 5a and Figure S2a). However, there were two regions (blue frames 1 and 2 in Figure 5b and Figure S2b) in SC3, where the PC was evidently lower than SC1 and SC2. One region is in the Qaidam Basin (blue frame 1), and the other is northeast of the QTP, around Qinghai Lake (blue frame 2). Both regions were located at low altitudes.

For RC, the number of grids with high values in SC1 was evidently higher than in SC2 and SC3 (Figure 6 and Figure S3). There were four regions in SC2 (black frames 1, 2, 3, and 4 in Figure 6a and Figure S3a) and three regions in SC3 (black frames 5, 6, and 7 in Figure 6b and Figure S3b), where RC was lower than SC1. For SC2, the first region was in the Qaidam Basin (black frame 1), the second region was between the Qaidam Basin and Tanggula Mountains (black frame 2), the third region was southwest of the QTP, including the mountains of Tanggula and Gangdese (black frame 3), and the fourth region was in the Hengduan Mountains on the eastern edge of the QTP (black frame 4). For SC3, the regions were the same as for SC2, except for the first one (black frame 1).

5. Discussion

5.1. Comparison with Previous Studies

Previous studies on the validation of daily cloud-free snow cover products were mostly conducted under all-weather conditions based on in situ snow depth data; that is, the validation results were a mixture of clear sky and cloud conditions, and the assessment accuracies for these products varied prominently from different studies. For SC2, the PC and RC reported by Huang et al. [28] were 87.1% and 79.6%, which were 27.0% points and 33.5% points higher than those reported by Hao et al. [38]. The PC and RC for SC3 reported by Qiu et al. [29] were 85.4% and 78.4%, which were 49.5% points and 11.5% points higher than those by Hao et al. [38]. The assessment of MOD10A1 or MOD10A1F focusing on the QTP has not yet been conducted. Hao et al. [39] assessed the accuracies of these two snow cover products under all-weather conditions with a PC and RC of over 85.0% for China, but the validation results for the QTP alone were not listed. The accuracies were much higher than those of SC1 in this study for PC (47.4%) and RC (64.0%). This may be related to the differences in the validation data and regions. The accuracy of MODIS snow cover products on the QTP was lower than in other areas in China [40].

The accuracies under all weather conditions in this study were close to those reported by Hao et al. [38]. For example, the PC and RC of SC2 under all-weather conditions were 48.0% and 50.1% in this study, which were 39.1% points and 29.5% points lower than the results in Huang et al. [28], but 12.1% points lower and 4.0% points higher than the PC and RC from Hao et al. [38]. In contrast, for SC3, the PC and RC under all-weather conditions were 47.8% and 57.3% in this study, which was 37.6% points and 21.1% points lower than the results reported by Qiu et al. [29] but 11.9% points higher and 9.6% points lower than the PC and RC from Hao et al. [38].

This study assessed the accuracy of cloud-gap-filled snow cover from commonly used daily snow cover products for the first time. Our study found that SC1 performed better than SC2 and SC3 (Table 3), particularly in regions with altitudes lower than 5000 m (Figure 3, Table 3), although the algorithm used in SC1 was the simplest among these products. SC1 indicated a better balance between commission and omission errors, particularly for snow cover (Table 3). The better accuracy for snow cover of SC1 was related to a higher PC than SC3 and better RC than SC2 and SC3 (Figure 3, Table 3). The improvement in the accuracies of SC1, compared with SC2 and SC3, indicated regional clustering characteristics (Figure 4, Figure 5 and Figure 6).

This study indicated that the error of all three products in the lower altitude areas was higher than that in higher altitude areas, based on Landsat TM snow cover maps (Figure 3). The accuracy assessments from previous studies were all based on in situ snow depth data, and most of the meteorological stations (approximately 85%) were located in the region with an altitude lower than 4000 m. Based on in situ snow depth data, these validation results represent the local area with low altitude on the QTP. Therefore, accuracies based on in situ snow depth data may underestimate the accuracy of snow cover products.

5.2. Limitations in the Determination of Cloud-Gap-Filled Snow Cover Based on Inferred Snow Lines

SC2 and SC3 used inferred snow lines to determine cloud pixel type (snow or non-snow). SC3 applied a multiple linear regression method to determine the expected snow line based on snowline samples. If the actual elevation of the cloud is higher than the expected snow lines, the cloud pixel is reclassified as snow; otherwise, it is reclassified as land. However, snow commonly occurs in regions under the snow line on the QTP [53]. Thus, the snow cover of SC3 may have been underestimated, and omission errors may have occurred at lower elevations (Figure 3a, Table 3).

Thus, compared to SC1 and SC2, the PC distribution for cloud-gap-filled snow cover showed that the lower PC of SC3 was located in the Qaidam Basin and northeast of the QTP around Qinghai Lake (Figure 4), which are the lowest areas on the QTP. For example, in the northern region with an altitude lower than 5000 m, the TM image from 6 January 2003 (Figure 7a) shows snow cover distribution (Figure 7b). SC1 (Figure 7f) and SC2 (Figure 7d) successfully reclassified the cloud pixels (Figure 7c) as snow cover, but SC3 (Figure 7e) misclassified almost all the cloud pixels as land.

The snow line (SNOWL) method [54] taken by SC2 classified the cloud as snow in the snow area and as land in the land area, where the cloud in the patch snow area was determined as snow or land by whether the snow actually existed or not based on the AMSR-E SWE [46]. This method decreases the omission errors of SC2 at lower elevations. Thus, the PC of SC3 was lower than that of SC2 (Table 3), although they both used snowline methods.

The high spatial heterogeneity of snow line elevation on the QTP induced more clustered errors from SC2 and SC3, compared with SC1. This is because the climate on the QTP transitioned from a tropical and subtropical humid climate in the southeast and south to a temperate arid climate in the northwest; thus, the altitude of the snow line varied greatly [55,56]. More cloud pixels are directly classified as snow cover when the actual snow line is higher than the calculated one, leading to more snow cover commission errors. This explains why the RC of SC2 and SC3 was lower than that of SC1 (Table 3) and why the errors mainly occurred in the mountains of Tanggula and Gangdese (black frames 3 and 6 in Figure 6a,b). This is because these areas are located on the northern slope of the Himalayas and the interior of the QTP, and the altitude of the snow line is higher than that of other regions because of less precipitation due to huge shielding from the mountains [57,58].

5.3. Impact of AMSR-E SWE on Cloud-Gap-Filled Snow-Cover Accuracy

The AMSR-E SWE product was used to reclassify the cloud pixels in the patch snow zone in SC2. However, the pixel size of the AMSR-E SWE is approximately 25 × 25 km, which is too coarse compared to the MODIS pixel size. As a result, the precision of the AMSR-E SWE product is not high, and the snow cover is easily overestimated over the QTP compared with other snow cover data [59,60]. This is because the dry frozen ground and cold deserts cannot be distinguished from snow by the snow retrieval algorithm for passive microwave remote sensing data [61,62,63]. This may explain why the PC for snow cover of SC2 was high (Figure 3a) but RC was the lowest among the three products (Figure 3b). Therefore, the cloud pixels have more chances of being classified as snow cover based on the AMSR-E SWE product.

In the Qaidam Basin, the RC of SC2 was lower than that of SC3 and SC1 (Figure 5). The snow cover may have been overestimated in the Qaidam Basin because of the impact of the AMSR-E SWE owing to the cold desert in winter in this area [61,62,63]. For example, the Landsat TM image in Figure 8a evidently shows that snow cover appears only in the mountains south of the Qaidam Basin. Except for SC2 (Figure 8d), there was no snow cover south of the mountains in the other products (Figure 8e,f). Therefore, the snow cover in the area highlighted by the red frame in Figure 8d is incorrect. The shape of the incorrect snow cover was close to that of clouds (the red frame in Figure 8c). Consequently, it can be inferred that the areas under the clouds were identified as patch snow by the SNOWL approach and then classified as snow cover based on the corresponding AMSR-E SWE data, which showed that the AMSR-E SWE was greater than zero (Figure 8g).

5.4. MOYD Cloud-Removal Algorithm Deficiency and Future Prospect

The temporal filter approach used by SC1 performed better than the snow line-based method, although the algorithm was simple (Table 3, Figure 3, Figure 4, Figure 5 and Figure 6). The cloud-removal algorithm adopted in SC1 uses a temporal filter approach in which the cloud pixels are replaced by the previous clean-sky pixels of the same location. This method can decrease the omission errors compared to the other products (Figure 3a, Table 3). However, this cloud-removal method may be ineffective during snow melting or accumulation periods, especially when the snowpack disappears or accumulates within only one or two days [64]. In other words, this method may not capture temporal snow under clouds. The accuracy of these products depends on the age of the observation, that is, the number of days since the last cloud-free observation [65].

Statistics-based methods for producing cloud-free snow cover products have recently been developed [32,33,35,39]. The spatial (neighborhood for a given pixel), temporal (time series), and environmental (topography and temperature) information were integrated into a statistical model (conditional probability and hidden Markov random field) to calculate the possibility of snow cover or non-snow cover for each cloud pixel. These methods achieved high accuracy for cloud removal based on the validation results of these studies, but these methods were not validated in regions other than the study area. Once the study area or period changed, the training data needed to be reselected, and the optimized parameters were recalibrated accordingly.

Thus, statistics-based methods are absent in generality, complicated, and computationally intensive. This study indicated higher accuracies of the temporal filter approach used in SC1 than the statistics-based method proposed by Huang et al. [28] and Qiu et al. [29]. Therefore, the temporal filter approach used in SC1 is the best solution to obtain daily cloud-gap-filled snow cover because of its higher accuracy and simple computation. Considering that in high-altitude regions (>5000 m), SC1 is not outstanding compared with other products, in the future, the spatial filter and environmental information should be assimilated into the cloud removal algorithm adopted in SC1 to improve the accuracy. Both the last and subsequent cloud-free observation can also be considered. In addition, integration with other satellite data of high spatial resolution such as Sentinel-1 active synthetic aperture radar (SAR) data may also be a promising solution for snow detection in snow melting and accumulation periods [66,67].

6. Conclusions

This study assessed the accuracy of cloud-gap-filled snow cover for commonly used daily snow-cover products. The results showed that SC1 was better than SC2 and SC3 and indicated a better balance between commission and omission errors, although the algorithm used in SC1 was the simplest among these products. The better accuracy for snow cover of SC1 was more prominent in regions with altitudes lower than 5000 m and indicated regional clustering characteristics compared with SC2 and SC3. The improvement in the accuracy of SC1, compared with SC2 and SC3, was related to lower omission errors compared with SC3 and commission errors compared with SC2 and SC3. The lower accuracies from SC2 and SC3 were mainly due to the limitation in the determination of cloud-gap-filled snow cover based on inferred snow lines and the overestimation of the AMSR-E SWE. These results indicate that the temporal filter approach used in SC1 should be the best solution to obtain daily cloud-gap-filled snow cover, owing to its higher accuracy and simple computation. However, this cloud-removal method may be ineffective during snow melting or accumulation periods, which needs to be further researched in the future. These findings can be helpful in the selection of cloud-removal algorithms to determine snow cover dynamics and phenological parameters.