Articles | Volume 26, issue 8
Research article
19 Apr 2022
Research article |  | 19 Apr 2022

Development and validation of a new MODIS snow-cover-extent product over China

Xiaohua Hao, Guanghui Huang, Zhaojun Zheng, Xingliang Sun, Wenzheng Ji, Hongyu Zhao, Jian Wang, Hongyi Li, and Xiaoyan Wang

Based on MOD09GA/MYD09GA surface reflectance data, a new MODIS snow-cover-extent (SCE) product from 2000 to 2020 over China has been produced by the Northwest Institute of Eco-Environment and Resources (NIEER), Chinese Academy of Sciences. The NIEER MODIS SCE product contains two preliminary clear-sky SCE datasets – Terra-MODIS and Aqua-MODIS SCE datasets and a final daily cloud-gap-filled (CGF) SCE dataset. The first two datasets are generated mainly through optimizing snow-cover discriminating rules over land-cover types, and the latter dataset is produced after a series of gap-filling processes such as aggregating the two preliminary datasets, reducing cloud gaps with adjacent information in space and time, and eliminating all gaps with auxiliary data. The validation against 362 China Meteorological Administration (CMA) stations shows that during snow seasons the overall accuracy (OA) values of the three datasets are larger than 93 %, all of the omission error (OE) values are constrained within 9 %, and all of the commission error (CE) values are constrained within 10 %. Bias values of 0.98, 1.02, and 1.03 demonstrate on a whole that there is no significant overestimation nor a significant underestimation. Based on the same ground reference data, we found that the new product accuracies are obviously higher than standard MODIS snow products, especially for Aqua-MODIS and CGF SCE. For example, compared with the CE of 23.78 % that the MYD10A1 product shows, the CE of the new Aqua-MODIS SCE dataset is 6.78 %; the OA of the new CGF SCE dataset is up to 93.15 % versus 89.54 % of MOD10A1F product and 84.36 % of MYD10A1F product. Besides, as expected, snow discrimination in forest areas is also improved significantly. An isolated validation at four forest CMA stations demonstrates that the OA has increased by 3–10 percentage points, the OE has dropped by 1–8 percentage points, and the CE has dropped by 4–21 percentage points. Therefore, our product has virtually provided more reliable snow knowledge over China; thereby, it can better serve for hydrological, climatic, environmental, and other related studies there.

1 Introduction

Snow is one of the most active elements on the land surface. Because of its unique properties of high shortwave reflectivity, low longwave emissivity, and high phase-change latent heat, its covering may significantly alter the surface radiation budget, exchanges of energy and moisture between the atmosphere and the surface, and thereby our climate and weather systems (Henderson et al., 2018; Huang et al., 2019; Warren, 1982). In addition, seasonal snow is an important supply source, providing precious freshwater for many arid and semiarid regions (Li et al., 2019). Snow-cover extent (SCE), therefore, is an indispensable parameter in climatic, weather, hydrological, environmental, and other related studies.

With the development and progress of space technology, satellite remote sensing has become a primary option to monitor snow-cover conditions on the Earth (Frei et al., 2012). In particular, as one of the most successful satellite missions, since it launched in 1999, the Moderate Resolution Imaging Spectroradiometer (MODIS) has been widely used to acquire SCE information from regional to global scales (Gafurov, et al., 2016; Hall and Riggs, 2007; Hall et al., 2002; Muhammad and Thapa, 2020). The reasons for this arise from the following: (1) MODIS is a multispectral sensor with a specific snow-detecting band around 1.64 µm; (2) features of nearly global coverage, double-satellite observations (Terra and Aqua), and suitable spatiotemporal resolutions also make it reasonable to acquire regional and global SCE.

The US National Snow and Ice Data Center (NSIDC) routinely produces and continually updates the standard MODIS snow products. Before the C6 version, there were only two sets of standard snow products – MOD10A1 and MYD10A1, which provide conventional SCE information only under clear skies. Since there are abundant cloud-induced gaps in the products, they in fact can not give the complete SCE knowledge. This is an obvious flaw for which reason the previous standard products were criticized most (Liang et al., 2008). As such, among the latest C6.1 version that was released very recently, another two sets of new cloud-gap-filled (CGF) products, MOD10A1F and MYD10A1F, are introduced and generated (Hall et al., 2019). But as of now, this update has not been completed yet, and these products are only available in some years.

In the standard MOD10A1 and MYD10A1 products, the core is datasets of the Normalized Difference Snow Index (NDSI) that is defined as

(1) NDSI = ( b4 - b6 ) / ( b4 + b6 ) ,

where b4 and b6 represent the reflectance of MODIS band 4 (0.55 µm) and band 6 (1.6 µm), respectively. The versions before C5 adopted NDSI ≥0.4 as the threshold to discriminate snow-cover pixels and produced binary snow-cover datasets (Riggs et al., 2006), whereas the versions after C6 only provide NDSI datasets but indicate 0.1 as the new threshold (Riggs et al., 2016). For clear skies, such a simple strategy is generally effective to distinguish snow cover, because the snow has a distinct NDSI characteristic relative to other common land covers. But there are also several appreciable shortcomings. First, several studies (Hall and Riggs, 2007; Hao et al., 2019; Zhang et al., 2019) have shown that the optimal NDSI threshold may vary with land-cover types, and using a fixed value for global products, regardless of 0.4 or 0.1, may lead to considerable uncertainty. Second, standard MOD10A1 and MYD10A1 products often exhibit large errors in forest areas (Maurer et al., 2003; Parajka et al., 2012; Poon and Valeo, 2006). For example, an evaluation conducted by Chen et al. (2014) indicates that in many forest regions of Northeast China the omission error or commission error may be up to 40 %. Last but not the least, for snow discrimination, theoretically using surface reflectance is also more reasonable than using top-of-atmosphere (TOA) reflectance (which is what the standard products do), as the atmospheric contribution may virtually reduce the NDSI contrasts between snows and other land covers (Hall et al., 1995).

Although the new MOD10A1F and MYD10A1F products are already able to give the complete daily snow information and thus would promote standard product applications dramatically, there are still several deficiencies, as far as our experience. First, because of MODIS tradition, they are produced for Terra-MODIS and Aqua-MODIS separately, but many studies (Huang et al., 2014; Gao et al., 2010; Gafurov and Bárdossy, 2009) have revealed that their synergy may be a better way to mitigate cloud influence. Second, their CGF mapping algorithm virtually only replaces cloud gaps in the current day with the previous most-recent clear-sky observation (Hall et al., 2019). Implications indicated by the subsequent clear-sky observation or spatially adjacent observations are all neglected. Under many scenarios, this seems unreasonable. For example, the snowfalls often happen on cloudy days, and if only the previous clear-sky observation is used, these days would be very likely mistaken for snow-free days. More advanced algorithms appearing in these years often infer snow-cover conditions beneath clouds comprehensively considering all available adjacent implications both in time and space (Huang et al., 2018).

Therefore, we decided to produce a new MODIS SCE product over China based on the Google Earth Engine (GEE) platform (Gorelick et al., 2017). Compared with the standard snow products, the following improvements are reached: (1) varying NDSI thresholds with land-cover types are obtained by a volume of training data; (2) the approach that combines the Normalized Difference Vegetation Index (NDVI) and the Normalized Difference Forest Snow Index (NDFSI, Wang et al., 2015, 2020) is adopted to improve snow discrimination in forest areas; (3) NDSI is computed using surface reflectance instead of TOA reflectance; (4) a gap-filling technique based on a hidden Markov random field (HMRF) (Huang et al., 2018) is used to reduce cloud-induced gaps, which can assimilate temporally and spatially adjacent information simultaneously.

2 Data

2.1 MODIS products

The MODIS products we use as the input data to generate new SCE data include the following: MOD09GA, MYD09GA, and MCD12Q1. MOD09GA and MYD09GA are the standard land surface reflectance products that are derived from Terra-MODIS and Aqua-MODIS, respectively, after the so-called atmospheric correction. They provide the 500 m land surface reflectance from MODIS bands 1 to 7, as well as the mask information (e.g., cloud and water masks), and they are our main inputting data. MCD12Q1 is the Terra\Aqua composite land-cover-type product, providing the annual land-cover information that is generated according to the International Geosphere–Biosphere Programme (IGBP) land-cover classification system. In this study, it is another important input which is used to indicate the detailed land-cover types. For all three products, the newest C6.1 version is adopted.

2.2 Landsat-8 OLI snow maps and MOD09GA\MYD09GA training samples

To refine the decision rules (which will be elaborated on in Sect. 3), we must obtain the quality MOD09GA and MYD09GA training samples in advance. Toward this, two groups of Landsat-8 Operational Land Imager (OLI) snow maps across China during the 2013–2018 snow seasons (beginning on 1 November through 31 March of the next year) are produced here. The first group is derived from 1509 scenes of OLI images, which will be regarded as “true” values to acquire the Terra-MODIS training samples; the second group comes from 1648 scenes of OLI images, which will be used to acquire the Aqua-MODIS training samples. These snow maps are generated by the improved “SNOMAP” algorithm developed by Chen et al. (2020), and they have a spatial resolution of 30 m. In every OLI snow map, there are only three classes – snow, snow free, and cloud.

Then 30 m snow maps will be aggregated within the 500 m spatial window to indicate whether the corresponding MOD09GA and MYD09GA pixel is covered by snow or not. Within a spatially and temporally (in the same day) colocated MOD09GA/MYD09GA pixel, if no less than 50 % of OLI pixels are snow covered, then it is a snow sample. If over 50 % of OLI pixels are snow free, it is a snow-free sample. If most of the OLI pixels are the cloud class, then it will be deemed as an invalid sample and is subsequently discarded. Of course, all of these must be done under the condition that the colocated MOD09GA/MYD09GA pixel is a clear pixel.

Finally, for MOD09GA, totally 21.20 million snow samples and 17.66 million snow-free samples are obtained; for MYD09GA, 12.05 million snow samples and 12.65 million snow-free samples are obtained in total. Note that it is necessary to obtain the Terra-MODIS and Aqua-MODIS training samples separately, as MYD09GA band 6 is not the directly observing data (many sensor detectors of this band have broken since the Aqua launch) but the restored data using the algorithm by Wang et al. (2006).

2.3 Auxiliary data

The auxiliary data we used include a snow-depth product generated from passive microwave satellite observations, a reanalysis land surface temperature (LST) product, and a digital elevation model (DEM) product.

The snow-depth product is available at (last access: 14 April 2022), providing daily and 0.25 snow-depth data over China from 1979 to 2020, which is generated by Che et al. (2008) and Dai et al. (2015) through a combination of multiple satellite passive-microwave observations. The reanalysis LST product is available on the GEE platform, providing daily and 0.25 LST data over China derived from ERA-5 global reanalysis (Munoz, 2019). The DEM product is one of the Shuttle Radar Topography Mission (SRTM) DEM products with a resolution of 90 m and directly accessible on the GEE platform. To match with MODIS data, all three products have been resampled or aggregated into 500 m.

2.4 Ground snow-depth measurements

Daily ground snow-depth observations up to 362 stations from the China Meteorological Administration (CMA) since 2000 will be used to validate or assess all associated MODIS SCE products in this study. To ensure the quality of the measurements, most CMA stations obey the following observing specifications: (1) snow depth is measured manually in an open spot near the station using a professional ruler, (2) the measurements were conducted only when the fractional snow cover in the field of view is larger than 50 %, (3) all observations were done at 08:00 LT (local time in Beijing) every morning.

At each station, the surface true condition, snow cover or snow free, is determined by the criterion proposed by Klein and Barnett (2003). That is, if measured snow depth is ≥1 cm, it is covered by snow; otherwise, it will be snow free. Because ground measurements are not limited by weather conditions, they are a better choice to independently validate all satellite-based SCE products.

3 Algorithm

3.1 Snow discrimination algorithm under clear skies

Guided by the algorithm of MODIS standard snow products (Hall et al., 2002; Riggs et al., 2006, 2016) and our motivations that are mentioned in Sect. 1, we developed a new snow discrimination algorithm for clear skies which is shown in Fig. 1. Approximately, it is divided into three steps to determine snow-cover conditions under clear skies from MODIS surface reflectance data. The first step is the preliminary screening with the purpose of precluding the cases completely impossibly covered by snow. The second step provisionally determines snow-cover conditions over non-forest land-cover types using the optimized NDSI thresholds and those over forest land-cover types through importing a new decision rule. Step three is the postprocessing based on surface temperature and DEM, which is designed to reverse those false snow pixels determined by the previous step into snow-free pixels.

Figure 1Flowchart of the new snow discrimination algorithm under clear skies.


3.1.1 Preliminary screening

As mentioned already, the purpose of the preliminary screening is to preclude the pixels that are impossibly covered by snow completely. Snow has a distinct spectral characteristic relative to other common land-cover types. Generally, its reflectance is high in the visible spectrum but rapidly drops in the infrared spectrum. As done by standard MODIS snow products (Riggs et al., 2006), we can use the combination of MODIS bands 2 and 4 within the visible spectrum and band 6 within the infrared spectrum to preliminary screen out the pixels that must be snow free but keep all possible snow pixels (even with a very low possibility) for a further discrimination.

For that purpose, we investigate all available snow samples and find for Terra-MODIS more than 99 % of the snow samples are constrained in the condition of band 2 ≥0.15, band 4 ≥0.05, and band 6 ≤0.45. Therefore, the preliminary screening rule of Terra-MODIS is adjusted for all possible snow pixels that must meet the conditions of band 2 ≥0.15, band 4 ≥0.05, and band 6 ≤0.45, and pixels that do not meet these conditions will be identified as snow free immediately. Similarly, for Aqua-MODIS 99 % of the snow samples are constrained in the conditions of band 2 ≥0.12, band 4 ≥0.07, and band 6 ≤0.40. The preliminary screening rule of Aqua-MODIS is set for all possible snow pixels that should meet these conditions, and pixels that do not meet these conditions will be deemed as snow free immediately.

3.1.2 Optimized NDSI thresholds over non-forest land-cover types

NDSI threshold usually plays a key role in snow discrimination under clear skies. Pixels whose NDSI value is greater than this threshold will be identified as snow cover; otherwise, snow-free pixels are assigned. Therefore, here the optimal NDSI threshold is crucial.

To obtain optimal NDSI thresholds for snow discrimination over different land-cover types, all of the training samples are first grouped according to their land-cover types that are indicated by the MCD12Q1 products. With the type of “barren or sparsely vegetated” as an example, the upper left of Fig. 2 presents the NDSI frequency distribution from Terra-MODIS training samples over this land cover, and the upper right presents the fluctuation of overall accuracy (OA; see Sect. 5 for details) with different NDSI thresholds. From the figure, one can see (around the crossing point of snows and barren lands) that the OA is highest (up to 98.34 %) when the NDSI threshold =0.08. Thus, 0.08 will be regarded as the optimal NDSI threshold for this land-cover type. Similarly, the bottom left and the bottom right of Fig. 2 present the NDSI frequency distribution and the fluctuation of OA derived from Aqua-MODIS training samples, respectively. It is apparent that for Aqua-MODIS the optimal NDSI threshold has changed to 0.06, and the highest OA is 95.20 %.

Figure 2Optimal NDSI thresholds over the land-cover type of “barren or sparsely vegetated”.


Following the same philosophy, the optimal NDSI thresholds for the other seven land-cover types, such as “grasslands”, “croplands”, “urban and built-up lands”, etc., are also obtained and listed in Table 1, together with their corresponding OAs. One can see that over these land-cover types the OAs are all larger than 95 %, which demonstrates the effectiveness of the new thresholds.

Table 1Optimal NDSI thresholds over eight non-forest land-cover types.

Download Print Version | Download XLSX

However, as we expected, only using the NDSI criterion seems not accurate enough to discriminate snow cover over those forest land-cover types, except the “evergreen broadleaf forest” (due to its sparse distributions in China). For example, Fig. 3 gives the Terra-MODIS NDSI frequency distribution and the fluctuation of OAs over the “evergreen needleleaf forest”. Obviously, here the commission error (CE; see Sect. 5 for details) is very large (>35 %), and the OAs drop severely (<76 %). Therefore, for the remaining seven forest land-cover types, we will adopt a new decision rule as elaborated next.

Figure 3Optimal NDSI thresholds over the land-cover type of “evergreen needleleaf forest”.


3.1.3 Optimized NDVI-NDFSI decision rules for forests

The study by Wang et al. (2015) had showed, compared with NDSI, that their so-called Normalized Difference Forest Snow Index (NDFSI, namely, using band 2 to substitute band 4 in Eq. 1) may be better for snow discrimination in forest areas. Our preliminary investigation reinforces this, and the improvement is, in particular, evident when using NDFSI in conjunction with NDVI (Wang et al., 2020).

Following their studies, NDVI in our study is divided into seven segments, and the optimal NDFSI at each NDVI segment will be computed with the highest OA as the object (similar to the acquirements of optimal NDSI thresholds). For a given pixel, if computed NDFSI is greater than the threshold in its corresponding NDVI segment, it will be regarded as a snow pixel; otherwise, it will be a snow-free pixel. With “deciduous broadleaf forest” as an example, Fig. 4 shows the NDVI-NDFSI number-density scatterplot that is derived from Terra-MODIS training samples over this land cover, as well as computed optimal NDFSI thresholds at every NDVI segment. Here the OA of the new NDVI-NDFSI decision rule is 99.91 % versus 88.77 % of the optimal NDSI threshold.

Figure 4NDVI-NDFSI number-density scatterplot over “deciduous broadleaf forest” from Terra-MODIS training samples (colors represent the number density within a 0.01×0.01 bin).


Using this approach, optimized NDVI-NDFSI decision rules for the six left forest land-cover types are all determined. Tables 2 and 3 give the detailed NDVI-depending NDFSI thresholds for Terra-MODIS and Aqua-MODIS, respectively. Besides these thresholds, the tables also give the OAs of the optimized NDVI-NDFSI decision rules in the second last column and the OAs of the optimal NDSI thresholds in the last column. We can see that, compared with optimal NDSI thresholds, optimized NDVI-NDFSI decision rules increase OAs greatly. Using the optimal NDSI thresholds, most OAs are ≤90 %, whereas when using the optimized NDVI-NDFSI decision rules, most OAs are ≥96 %, in line with the accuracies of non-forest land-cover types in Table 1.

Table 2Optimized Terra-MODIS NDVI-NDFSI decision rules for forest land-cover types (ENF: evergreen needleleaf forests; DNF: deciduous needleleaf forests; DBF: deciduous broadleaf forests; PWS: permanent wetland savannas).

a refers to the OA of optimized NDVI-NDFSI decision rules. b refers to the OA of optimal NDSI thresholds.

Download Print Version | Download XLSX

Table 3Optimized Aqua-MODIS NDVI-NDFSI decision rules for forest land-cover types.

Download Print Version | Download XLSX

3.1.4 Postprocessing based on surface temperature and DEM

Ice clouds have similar optical properties to snows. From time to time, pixels contaminated by broken ice clouds are mistaken for snow pixels, which will lead to some false snow existing even in the tropics. Like standard MODIS products (Riggs et al., 2006), we also use postprocessing based on surface temperature and DEM to reverse these snow pixels into snow-free pixels, but here ERA5 reanalysis LST (Munoz, 2019) is used rather than the brightness temperature from band 31.

According to our previous investigation (Hao et al., 2021), for lowlands of DEM <1300 m, snow at the surface basically impossibly exists when its temperature is ≥275 K (2 C), but for highlands of DEM ≥1300 m, a higher temperature threshold, 281 K (8 C), seems more appropriate due to the possible existence of warm snow on highlands (Riggs et al., 2016). Thereby, if snow is detected in a pixel when its height <1300 m and LST ≥275 K, this snow decision would be reversed to snow free. If snow is detected in a pixel when its height ≥1300 m and LST ≥281 K, this snow decision would also be reversed to snow free. Such a postprocessing is very effective to alleviate the problem of false snow detections, which are often seen in summers or over south China (when and where snow is impossible).

3.2 Cloud-gap-removing algorithm

Figure 5 describes the flow of the new cloud-gap-removing algorithm we developed in this study. It can be divided into three steps. The first step is preliminarily excluding some cloud gaps by the synergy of Terra-MODIS and Aqua-MODIS; the second step is further filling gaps according to the implication of the nearby clear-sky pixels using the hidden Markov random field (HMRF) technique; and in step three all left gaps will be filled using the auxiliary passive microwave product.

Figure 5Flowchart of the new cloud-gap-removing algorithm.


3.2.1 Aggregation of Terra-MODIS and Aqua-MODIS

There may be different weather conditions when Terra and Aqua pass in a day, and their synergy can reduce cloud gaps approximately by 10 %–30 % over China (Huang et al., 2016). Thus, in this study, we will first aggregate the Terra-MODIS SCE and Aqua-MODIS SCE under clear skies to preliminarily exclude some cloud gaps. The aggregating rule is as follows. In a day, if Terra-MODIS indicates a clear sky, its SCE information will be adopted; Otherwise, if Aqua-MODIS indicates a clear sky, Aqua-MODIS SCE information will be adopted. If both indicate a cloudy sky, the cloud gap is kept in that day.

3.2.2 Gap-filling based on the HMRF technique

For the aggregated SCE data, we will adopt the method based on the spatiotemporal hidden Markov random field (HMRF) proposed by Huang et al. (2018) to further fill gaps. The core of this method is computing the total spatiotemporal “energy” (probability) of belonging to snow pixels or snow-free pixels at one specific gap, using


where Ust is the spatiotemporal neighborhood energy function. Ns and Nt denote the spatial neighborhood and temporal neighborhood cantered with the gap, respectively.

Figure 6 illustrates our gap-filling process based on the HMRF technique. For a given gap, we will first calculate the U(Snow) and U(Snow-free) based on the 3 (row) ×3 (column) × 3 (day) spatiotemporal neighborhood. If the U(Snow) is greater than the U(Snow-free), the gap will be classified as snow pixels; otherwise, it will be classified as snow-free pixels. If there are not sufficient valid pixels for calculating the U(Snow) or U(Snow-free), we will extend the spatiotemporal neighborhood to 3 (row) × 3 (column) × 5 (day). If there are still not sufficient valid pixels, the spatiotemporal neighborhood will expand into 5 (row) × 5 (column) × 5 (day). If the above tries all fail, the gap will be kept still.

Figure 6Diagram of the HMRF gap-filling process used in this study (dark grey represents high weight, and light grey represents low weight).


The gap-filling technique based on the HMRF provides a rigorous analytical framework for integrating spatial and temporal contexts in an optimal manner, which is accurate and highly efficient for MODIS cloud-gap filling (Huang et al., 2018). Our investigation shows that this method can remove approximately 80 % of the gaps in the aggregated data.

3.2.3 Eliminating residual gaps with auxiliary data

Although most cloud gaps have been filled after the two above processes, there are still a few exceptions. In this phase, we will use the auxiliary passive microwave snow-depth product provided by Che et al. (2008) to eliminate all remaining gaps. If colocated snow depth is ≥2 cm, then the gap will be assigned a snow pixel; otherwise, it will be a snow-free pixel (Hao et al., 2019).

4 Product

Using the new snow discrimination algorithm developed for clear skies in Sect. 3.1, we produce a new Terra-MODIS SCE dataset and a new Aqua-MODIS SCE dataset over China with MOD09GA and MYD09GA products as the main inputs. There are four classes of pixels in the new datasets, namely, snow, snow free, water body, and gap (see Fig. 7). Because optical remote sensing is affected by clouds severely, there are often abundant gaps (up to 70 %) in such clear-sky products, which limit their application to a very great extent.

Figure 7NIEER MODIS SCE product on 4 January 2020: (a) Terra-MODIS SCE dataset; (b) Aqua-MODIS SCE dataset; (c) CGF-MODIS SCE dataset.

Next, based on the two clear-sky datasets above, we produce a new CGF SCE dataset using the new cloud-gap-removing algorithm developed in Sect. 3.2. In the new dataset, all gaps have been removed, and only three classes (snow, snow free, and water body) are included (see Fig. 7). Comparing with the first two datasets, this dataset undoubtedly has a wider application prospect.

The new Terra-MODIS, Aqua-MODIS, and CGF-MODIS SCE datasets together constitute the NIEER MODIS SCE product (Hao et al., 2021). It has a 500 m spatial resolution and a daily temporal resolution, providing the snow-cover condition in the past two decades over China. An instance of snow coverage revealed by the product is presented in Fig. 7.

5 Validation

A confusion matrix similar to Table 4 is used to validate or assess all associated MODIS SCE data in this study. The overall accuracy (OA) represents the percentage that snow and snow-free pixels are correctly classified, the producer's accuracy (PA, corresponding to the omission error, OE=1-PA) represents the percentage that real snow pixels are classified as such, and the user's accuracy (UA, corresponding to the commission error, CE=1-UA) represents the percentage that pixels classified as snow are also snow in reality. Because the kappa coefficient has caused too much controversy recently (Pontius and Millones, 2011; Stehman and Foody, 2019), we no longer use it as an accuracy metric here. Meanwhile, to more directly show whether SCE is overestimated or underestimated and to what extent, a metric termed as “bias” is introduced. The bias, defined as the ratio of pixels classified as snow to real snow pixels (Hüsler et al., 2012; Wu et al., 2021), is a relative measure to detect overestimation (values >1) or underestimation of snow (values <1).

Table 4Description of the confusion matrix used to validate or assess the associated SCE products (SS, SN, NS, and NN are the numbers in each case).

Download Print Version | Download XLSX

Note that nominally 20-year and 362-station snow-depth measurements will be used as reference data to validate our product, but in fact validations are conducted only during the snow seasons at one station when the snow-cover days are ≥20; otherwise, some accuracy metrics would become artificially high due to too many non-snow pixels in a mid-latitude region like China (Hao et al., 2008).

5.1 Overview of the new product's accuracies

In Table 5, an overview is given of validation results for the new MODIS SCE product. We can see on a whole that the two new clear-sky datasets are very accurate. For the Terra-MODIS SCE dataset, the OA, OE, and CE are 95.51 %, 5.56 %, and 8.26 %, respectively; for the Aqua-MODIS SCE dataset, the OA, OE, and CE are 94.73 %, 8.47 %, and 6.78 %, respectively. The former has a lower omission error, while the latter has a lower commission error. From the bias point of view, the Terra-MODIS dataset slightly overestimates snow cover, and the Aqua-MODIS dataset slightly underestimates snow cover. But from the specific values (1.03 and 0.98), the overestimation and underestimation are marginal. Comparatively speaking, the Aqua-MODIS dataset is inferior to the Terra-MODIS dataset. This is not surprising considering the problem of the Aqua-MODIS band 6 and the fact that CMA snow-depth measurement is only conducted in the mornings.

Although the CGF-MODIS SCE dataset is slightly worse than the two clear-sky datasets, there is no significant difference no matter what accuracy metric is used. Here the OA, OE, and CE are 93.15 %, 8.25 %, and 9.83 %, respectively. Compared with the metrics shown above, the accuracy drops are basically within 1–2 percentage points. Moreover, the bias of 1.02 also demonstrates an inconspicuous overestimation. Therefore, our cloud-gap-removing processes are effective, and the accuracy may be sustained at a consistent level as the clear-sky datasets.

Table 5Validation results of the new NIEER MODIS SCE product with reference to CMA ground measurements.

Download Print Version | Download XLSX

5.2 The stability in time

To assess the stability of the NIEER MODIS SCE product in time, the accuracy metrics in every snow season since 2000 are calculated separately. Figure 8 presents the interannual fluctuation of these accuracy metrics. From the figure, we can see that all three datasets have a similar fluctuation characteristic. They always perform a consistent high accuracy in one snow season, and a coincident low accuracy in another snow season. Such an accuracy fluctuation may be attributed to the different validation numbers available in different years caused by ground measurements.

Figure 8Annual accuracy metrics' fluctuation of the NIEER MODIS SCE product: (a) OA, (b) PA, (c) UA, and (d) Bias.


The highest accuracy happens in the snow season of 2012, where the OEs of the three datasets are 2.84 %, 5.76 %, and 4.60 %, and the CEs are 5.02 %, 3.72 %, and 6.01 %. Meanwhile, the worst accuracy happens in the snow season of 2018, where the OEs of the three datasets are 9.49 %, 12.91 %, and 14.95 %, and the CEs are 15.42 %, 14.96 %, and 17.15 %. But except this snow season, the OEs and CEs in other seasons are all less than 15 %. Relative to OE and CE (namely, PA and UA), OA fluctuations of the three datasets are smaller. The Terra-MODIS dataset fluctuates within 93 %–97 %, the Aqua-MODIS dataset fluctuates within 93 %–96 %, and the CGF-MODIS dataset fluctuates within 92 %–95 %.

In addition, from the bias point of view, the values ranging from 0.94 to 1.07 demonstrate that there are neither substantial overestimations nor substantial underestimations in all datasets and snow seasons. Therefore, on average the new product is stable in time and promising to better serve the climatic, hydrological, and other related studies in China (Huang et al., 2020).

5.3 The reliability in space

To clarify the new product's reliability in space, we calculate accuracies at each CMA station. With the Terra-MODIS SCE dataset as an example, Fig. 9 presents the detailed accuracies. From the figure, we find, in spite of a very high accuracy seen at most stations, the errors may be relatively larger in the southwest of Northeast China and the northeastern edges of Qinghai–Tibet Plateau. This can be attributed to their patchy snow-cover features and rugged terrains (Xiao et al., 2020), where Loess Plateau and Qinghai–Tibet Plateau rapidly transit toward Sichuan Basin and Northeast Plain, respectively. The Aqua-MODIS SCE dataset performs a very similar accuracy pattern like Fig. 9 but slightly worse.

Figure 9Accuracies of the Terra-MODIS SCE dataset at each CMA station: (a) OA, (b) PA, (c) UA, and (d) Bias.

Figure 10 further details the accuracies of the CGF-MODIS SCE dataset at all stations. From the figure, we see that in North Xinjiang and the north of Northeast China, where the snowpack is stable, the accuracy is high, whereas in the northeast of Inner Mongolia Autonomous Region, the northwest of North China, and Qinghai–Tibet Plateau, where snow may melt rapidly even in winter, the accuracy is relatively lower. Besides the propagated errors from the two original clear-sky SCE datasets, it is supposed that an unstable snowpack in warmer areas would result in another large uncertainty.

Figure 10Accuracies of the CGF-MODIS SCE dataset at each CMA station: (a) OA, (b) PA, (c) UA, and (d) Bias.

6 Discussion

6.1 Comparisons with standard MODIS products

Using the same ground reference data, we also evaluate the standard SCE data derived from the newest C6.1 MOD10A1, MYD10A1, MOD10A1F, and MYD10A1F products (see Table 6). Comparing the accuracy metrics listed here and those in Table 5, we can see that these products' accuracies are significantly lower than our product's accuracy. If ground measurements are seen as true values, the MOD10A1 product's OE and CE values are 7.22 % and 12.97 %, respectively, versus 5.56 % and 8.26 % of our Terra-MODIS dataset; the MYD10A1 product's OE and CE values are 13.78 % and 23.78 %, respectively, versus 8.47 % and 6.78 % of our Aqua-MODIS dataset. Our improvement is particularly significant for Aqua-MODIS SCE, where the CE is improved by over 17 percentage points. Besides, from the bias point of view, there are appreciable overestimations to snow cover in standard MOD10A1 and MYD10A1 products.

Note that here factual improvement may be much more obvious than shown by the accuracy metric of OA, because there are a large number of non-snow pixels (far beyond snow pixels) in a mid-latitude region like China, which would pull it up dramatically. Although the new product's improvement in OA is within 1 % relative to MOD10A1 products, in fact the OE has decreased by 37 %, and the CE has decreased by 22 %.

Table 6Accuracies of SCE products derived from MOD10A1, MYD10A1, MOD10A1F, and MYD10A1F using the same ground measurements.

Download Print Version | Download XLSX

Due to possible wider application, the CGF product comparison is emphasized here. Relative to the better standard product (MOD10A1F), our product's OA increases by nearly 4 percentage points, the OE drops from 11.04 % to 8.43 %, and the CE drops from 13.26 % to 9.83 %; relative to the worse standard product (MYD10A1F), the OA even increases by nearly 9 percentage points, the OE drops from 21.07 % to 8.25 %, and the CE drops from 14.47 % to 9.83 %. It is clear for the standard CGF products that all of the OEs and CEs are larger than 10 %, while for our product the OE and CE are both within 10 %. Meanwhile, the two standard CGF products' biases are also further away from 1. Therefore, our improvement to CGF SCE is obvious, too. Besides, comparing the two standard clear-sky products with the two standard CGF products, we find their cloud-gap-removing strategy indeed lowers the accuracies dramatically.

6.2 Improvement in forest areas

To figure out the improvement in forest areas, the validations and comparisons at four forest CMA stations are isolated and listed in Table 7. From the table, one can see the OA increases by 3–10 percentage points, the OE drops by 1–8 percentage points, and the CE drops by 4–21 percentage points. As well, more significant improvements are seen for the Aqua-MODIS and the CGF SCE datasets. For example, the CE of the new Aqua-MODIS SCE dataset has already dropped to 9.62 %, comparing with 30.13 % of MYD10A1 product; the OA of the new CGF-MODIS dataset increases to 91.23 %, compared with 81.79 % of MOD10A1F product. Thus, the improvement of the new product in forest areas is also considerably significant, even more significant than in non-forest areas.

Table 7Validation and comparison at four forest CMA stations.

Download Print Version | Download XLSX

Two more intuitional examples are shown in Fig. 11. One is a typical forest area in Northeast China on 1 January 2018; another comes from southeastern forest regions of Qinghai–Tibet Plateau on 1 February 2018. It is clear that the new CGF-MODIS SCE dataset agrees much better with the higher-resolution snow maps than the standard CGF product does. The standard MOD10A1F perceivably overestimates snow coverage in these forest areas.

Figure 11Intuitional improvements of CGF SCE in two representative forest areas: (a) an example in a Northeast China forest region; (b) an example in the Qinghai–Tibet Plateau forest region.

Of course, the new product is far from perfect. In particular, the product's performance is limited severely by the cloud mask provided by standard MOD09GA/MYD09GA products, which is virtually caused by the cloud/snow confusion problem. On the one hand, some broken ice clouds are susceptible to being mistaken for snow pixels because of their similar optical properties, which will result in some artificial snow pixels in South China (although it has been mitigated largely compared with the standard MODIS snow products). On the other hand, snow pixels are also possibly mistaken for clouds, which will result in some omitted snow covers. During our validations or comparisons, we found this phenomenon is somewhat common in the edges of snow-cover areas and the forest areas of Northeast China.

7 Summary and outlook

In this study, we produce a new MODIS SCE product over China, which is committed to addressing currently known problems of the standard snow products. Therefore, the optimal NDSI thresholds varying with land-cover types are extracted, and the NDVI-NDFSI decision rules proposed by Wang et al. (2015) specific for snow discrimination in forest areas are optimized by a number of training samples indicated by the higher-resolution snow maps from Landsat-8 OLI images. Meanwhile, to obtain a daily continuous cloud-free SCE dataset, the clear-sky snow SCE data from Terra-MODIS and Aqua-MODIS are aggregated first; an HMRF gap-filling technique, which can simultaneously assimilate temporally and spatially neighboring information, is imported second; finally, the residual gaps are all filled according to the implication given by an auxiliary passive microwave snow-depth dataset.

The new product is named the NIEER MODIS SCE product, which contains three individual datasets – the Terra-MODIS SCE dataset, the Aqua-MODIS SCE dataset, and the MODIS CGF SCE dataset. The first two provide the SCE knowledge under clear skies directly derived from the MODIS surface reflectance products, while the last provides daily and totally continuous SCE knowledge.

The comprehensive validations against 362 CMA stations across China have revealed that our product is not only very accurate but also considerably stable in time and reliable in space. On a whole, there are no significant overestimations nor any significant underestimations in all datasets, and the cloud-gap-removing processes do not lower the SCE accuracy severely. Compared with SCE indicated by the standard MODIS snow products, our product's accuracies are obviously higher, especially for the Aqua-MODIS and CGF SCE datasets. Relative to the MYD10A1 product, the CE drops 17 percentage points, from 23.78 % to 6.78 %. Relative to the better MOD10A1F CGF product, the OA increases by nearly 4 percentage points, from 89.54 % to 93.15 %. Meanwhile, as we expected, the improvement in forest areas is also evident. If validations at forest stations are isolated, OA increases by 3–15 percentage points. Therefore, the new product has provided more reliable knowledge on snow coverage over China and thereby will be a better choice for climatic, hydrological, and other related studies there.

The problem of cloud/snow confusion may contribute to the largest uncertainty in the new product. There are many cases when snow pixels are mistaken for cloudy pixels or the opposite case due to the inaccurate cloud mask provided by MOD09GA/MYD09GA products. In the future, we will consider further improving our product from this aspect.

Code and data availability

The new MODIS SCE product is available from the National Tibetan Plateau Data Center: (Hao, 2021). The code used to generate the product can be obtained from the authors without conditions.

Author contributions

XH, GH, and ZZ conceived and designed the study. XS, WJ, and HZ developed the code and generated the product on the GEE platform. GH and XH prepared the initial draft of the article. All authors contributed to the data analysis, results interpretation and validation, and revision of the article.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We would like to thank the China Meteorological Administration (CMA) for their reliable ground measurements and the Google Earth Engine (GEE) for its high-quality services.

Financial support

This research has been supported by the National Key Research and Development Program of China (grant no. 2019YFC1510503), the National Natural Science Foundation of China (grant nos. 41971325, 42171322, and 42171391), and the National Key Research and Development Program of China (grant no. 2017FY100502).

Review statement

This paper was edited by Hongkai Gao and reviewed by four anonymous referees.


Che, T., Li, X., Jin, R., Armstrong, R., and Zhang, T. J.: Snow depth derived from passive microwave remote-sensing data in China, Ann. Glaciol., 49, 145–154,, 2008. 

Chen, S., Wang, X., Guo, H., Xie, P., Wang, J., and Hao, X.: A Conditional Probability Interpolation Method Based on a Space-Time Cube for MODIS Snow Cover Products Gap Filling, Remote Sens., 12, 3577,, 2020. 

Chen, S. B., Yang, Q., Xie, H. J., Zhou, C., and Lu, P.: Time series of snow cover data of Northeast China (2004–2013), Acta Geographica Sinica, 69, 178–184, 2014. 

Dai, L. Y., Che, T., and Ding, Y. J.: Inter-calibrating SMMR, SSM/I and SSMI/S data to improve the consistency of snow-depth products in China, Remote Sens., 7, 7212–7230,, 2015. 

Frei, A., Tedesco, M., Lee, S., Foster, J., Hall, D. K., Kelly, R., and Robinson, D. A.: A review of global satellite-derived snow products, Adv. Space Res., 50, 1007–1029,, 2012. 

Gafurov, A. and Bárdossy, A.: Cloud removal methodology from MODIS snow cover product, Hydrol. Earth Syst. Sci., 13, 1361–1373,, 2009. 

Gafurov, A., Lüdtke, S., Unger-Shayesteh, K., Vorogushyn, S., Schöne, T., Schmidt, S., Kalashnikova, O., and Merz, B.: MODSNOW-Tool: an operational tool for daily snow cover monitoring using MODIS data, Environ. Earth Sci., 75, 1078,, 2016. 

Gao, Y., Xie, H. J., Yao, T. D., and Xue, C. S.: Integrated assessment on multi-temporal and multi-sensor combinations for reducing cloud obscuration of MODIS snow cover products of the Pacific Northwest USA, Remote Sens. Environ., 114, 1662–1675,, 2010. 

Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., and Moore, R.: Google Earth Engine: Planetary-scale geospatial analysis for everyone, Remote Sens. Environ., 202, 18–27,, 2017. 

Hall, D. K. and Riggs, G. A.: Accuracy assessment of the MODIS snow products, Hydrol. Process., 21, 1534–1547,, 2007. 

Hall, D. K., Riggs, G. A., and Salomonson, V. V.: Development of Methods for Mapping Global Snow Cover Using Moderate Resolution Imaging Spectroradiometer Data, Remote Sens. Environ., 54, 127–140,, 1995. 

Hall, D. K., Riggs, G. A., Salomonson, V. V., DiGirolamo, N. E., and Bayr, K. J.: MODIS snow-cover products, Remote Sens. Environ., 83, 181–194,, 2002. 

Hall, D. K., Riggs, G. A., DiGirolamo, N. E., and Román, M. O.: Evaluation of MODIS and VIIRS cloud-gap-filled snow-cover products for production of an Earth science data record, Hydrol. Earth Syst. Sci., 23, 5227–5241,, 2019. 

Hao, X.: A new MODIS snow cover extent product over China (2000–2020), National Tibetan Plateau Data Center (NTPDC) [data set],, 2021. 

Hao, X., Huang, G., Che, T., Ji, W., Sun, X., Zhao, Q., Zhao, H., Wang, J., Li, H., and Yang, Q.: The NIEER AVHRR snow cover extent product over China – a long-term daily snow record for regional climate research, Earth Syst. Sci. Data, 13, 4711–4726,, 2021. 

Hao, X. H., Wang, J., and Li, H. Y.: Evaluation of the NDSI threshold value in mapping snow cover of MODIS – a case study of snow in the Middle Qilian Mountains, Journal of Glaciology and Geocryology, 30, 132–138, 2008. 

Hao, X. H., Luo, S. Q., Che, T., Wang, J., Li, H. Y., Dai, L. Y., Huang, X. D., and Feng, Q. S.: Accuracy assessment of four cloud-free snow cover products over the Qinghai-Tibetan Plateau, Int. J. Digit. Earth, 12, 375–393,, 2019. 

Henderson, G. R., Peings, Y., Furtado, J. C., and Kushner, P. J.: Snow-atmosphere coupling in the Northern Hemisphere, Nat. Clim. Change, 8, 954–963,, 2018. 

Huang, G. H., Li, Z. Q., Li, X., Liang, S. L., Yang, K., Wang, D. D., and Zhang, Y.: Estimating surface solar irradiance from satellites: Past, present, and future perspectives, Remote Sens. Environ., 233, 111371,, 2019. 

Huang, G. H., Li, X., Lu, N., Wang, X. F., and He, T.: A General Parameterization Scheme for the Estimation of Incident Photosynthetically Active Radiation Under Cloudy Skies, IEEE T. Geosci. Remote, 58, 6255–6265,, 2020. 

Huang, X., Deng, J., Ma, X., Wang, Y., Feng, Q., Hao, X., and Liang, T.: Spatiotemporal dynamics of snow cover based on multi-source remote sensing data in China, The Cryosphere, 10, 2453–2463,, 2016. 

Huang, X. D., Hao, X. H., Feng, Q. S., Wang, W., and Liang, T. G.: A new MODIS daily cloud free snow cover mapping algorithm on the Tibetan Plateau, Sciences in Cold and Arid Regions, 6, 116–123, 2014. 

Huang, Y., Liu, H. X., Yu, B. L., We, J. P., Kang, E. L., Xu, M., Wang, S. J., Klein, A., and Chen, Y. N.: Improving MODIS snow products with a HMRF-based spatio-temporal modelling technique in the Upper Rio Grande Basin, Remote Sens. Environ., 204, 568–582,, 2018. 

Hüsler, F., Jonas, T., Wunderle, S., and Albrecht, S.: Validation of a modified snow cover retrieval algorithm from historical 1-km AVHRR data over the European Alps, Remote Sens. Environ., 121, 497–515,, 2012. 

Klein, A. G. and Barnett, A. C.: Validation of daily MODIS snow cover maps of the Upper Rio Grande River Basin for the 2000-2001 snow year, Remote Sens. Environ., 86, 162–176,, 2003. 

Li, H. Y., Li, X., Yang, D. W., Wang, J., Gao, B., Pan, X. D., Zhang, Y. L., and Hao, X. H.: Tracing snowmelt paths in an integrated hydrological model for understanding seasonal snowmelt contribution at basin scale, J. Geophys. Res.-Atmos., 124, 8874–8895,, 2019. 

Liang, T. G., Huang, X. D., Wu, C. X., Liu, X. Y., Li, W. L., Guo, Z. G., and Ren, J. Z.: An application of MODIS data to snow cover monitoring in a pastoral area: A case study in Northern Xinjiang, China, Remote Sens. Environ., 112, 1514–1526,, 2008. 

Maurer, E. P., Rhoads, J. D., Dubayah, R. O., and Lettenmaier, D. P.: Evaluation of the snow-covered area data product from MODIS, Hydrol. Process., 17, 59–71,, 2003. 

Muhammad, S. and Thapa, A.: An improved Terra–Aqua MODIS snow cover and Randolph Glacier Inventory 6.0 combined product (MOYDGL06*) for high-mountain Asia between 2002 and 2018, Earth Syst. Sci. Data, 12, 345–356,, 2020. 

Munoz, S. J.: ERA5-Land hourly data from 1981 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set],, 2019. 

Parajka, J., Holko, L., Kostka, Z., and Blöschl, G.: MODIS snow cover mapping accuracy in a small mountain catchment – comparison between open and forest sites, Hydrol. Earth Syst. Sci., 16, 2365–2377,, 2012.  

Pontius Jr., R. G. and Millones, M.: Death to Kappa: birth of quantity disagreement and allocation disagreement for accuracy assessment, Int. J. Remote Sens., 32, 4407–4429,, 2011. 

Poon, S. K. M. and Valeo, C.: Investigation of the MODIS snow mapping algorithm during snowmelt in the northern boreal forest of Canada, Can. J. Remote Sens., 32, 254–267,, 2006. 

Riggs, G. A., Hall, D. K., and Salomonson, V. V.: MODIS Snow products user guide to collection 5, NSIDC – US National Snow and Ice Data Center, Boulder, CO, (last access: 13 April 2022), 2006. 

Riggs, G. A., Hall, D. K., and Roman, M. O.: MODIS snow products user guide for collection 6, NSIDC – US National Snow and Ice Data Center, Boulder, CO, (last access: 13 April 2022), 2016. 

Stehman, S. V. and Foody, G. M.: Key issues in rigorous accuracy assessment of land cover products, Remote Sens. Environ., 231,111199,, 2019. 

Wang, L. L., Qu, J. J., Xiong, X. X., Hao, X. J., Xie, Y., and Che, N. Z.: A new method for retrieving band 6 of Aqua MODIS, IEEE T. Geosci. Remote S., 3, 267–270,, 2006. 

Wang, X. Y., Wang, J., Jiang, Z. Y., Li, H. Y., and Hao, X. H.: An effective method for snow-cover mapping of dense coniferous forests in the upper Heihe River Basin using Landsat Operational Land Imager data, Remote Sens., 7, 17246–17257,, 2015. 

Wang, X. Y., Chen, S. Y., and Wang, J.: An adaptive snow identification algorithm in the forests of Northeast China, IEEE J-Stars, 13, 5211–5222,, 2020. 

Warren, S. G.: Optical-Properties of Snow, Rev. Geophys., 20, 67–89,, 1982. 

Wu, X., Naegeli, K., Premier, V., Marin, C., Ma, D., Wang, J., and Wunderle, S.: Evaluation of snow extent time series derived from Advanced Very High Resolution Radiometer global area coverage data (1982–2018) in the Hindu Kush Himalayas, The Cryosphere, 15, 4261–4279,, 2021. 

Xiao, P. F., Li, C. X., Zhu, L. J., Zhang, X. L., Ma, T. Y., and Feng, X. Z.: Multitemporal ensemble learning for snow cover extraction from high-spatial-resolution images in mountain areas, Int. J. Remote Sens., 41, 1668–1691,, 2020. 

Zhang, H. B., Zhang, F., Zhang, G. Q., Che, T., Yan, W., Ye, M., and Ma, N.: Ground-based evaluation of MODIS snow cover product V6 across China: Implications for the selection of NDSI threshold, Sci. Total Environ., 651, 2712–2726,, 2019. 

Short summary
We develop and validate a new 20-year MODIS snow-cover-extent product over China, which is dedicated to addressing known problems of the standard snow products. As expected, the new product significantly outperforms the state-of-the-art MODIS C6.1 products; improvements are particularly clear in forests and for the daily cloud-free product. Our product has provided more reliable snow knowledge over China and can be accessible freely