The recent developments in cloud removal approaches of MODIS snow cover product

The snow cover products of optical remote sensing systems play an important role in research into global climate change, the hydrological cycle, and the energy balance. Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover products are the most popular datasets used in the community. However, for MODIS, cloud cover results in spatial and temporal discontinuity for long-term snow monitoring. In the last few decades, a large number of cloud removal methods for MODIS snow cover products have been proposed. In this paper, our goal is to make a comprehensive summarization of the existing algorithms for generating cloud-free MODIS snow cover products and to expose the development trends. The methods of generating cloud-free MODIS snow cover products are classified into spatial methods, temporal methods, spatio-temporal methods, and multi-source fusion methods. The spatial methods and temporal methods remove the cloud cover of the snow product based on the spatial patterns and temporal changing correlation of the snowpack, respectively. The spatio-temporal methods utilize the spatial and temporal features of snow jointly. The multi-source fusion methods utilize the complementary information among different sources among optical observations, microwave observations, and station observations.

1 Introduction

Because of the high albedo, high thermal emissivity, low thermal conductivity, and water storage ability (Tait et al., 2000; Tekeli and Tekeli, 2012), snow has a significant influence on the energy balance (Robinson et al., 1993; Crawford et al., 2013), the hydrological cycle (Şorman et al., 2007; Kostadinov and Lookingbill, 2015), and climate change (Cohen and Entekhabi, 1999; Brown, 2000). In recent years, increasing attention has been focused on monitoring the spatial and temporal change of snow cover (Gao et al., 2012). The traditional in situ snow monitoring approach is conducted only sparsely and is limited due to the large gaps in both space and time (Brown and Braaten, 1998). In contrast, remote sensing data have the advantage of a wide coverage area and relatively short revisit period (Lopez et al., 2008; Zeng et al., 2013) and have been an effective and alternative supplement for in situ data since April 1960 through the Television Infrared Observation Satellite (TIROS; Singer and Popham, 1963). According to the data source, snow cover products based on remote sensing mainly include microwave-based products, optically based products, and combined products (Frei et al., 2012), as shown in Fig. 1.

Figure 1 Snow cover products based on remote sensing.

Microwave-based products are derived according to the relationship between microwave energy and snow depth (SD) or the snow water equivalent (SWE) when the snowpack is dry (Tait et al., 2000; Wulder et al., 2007). Microwave-based products are free from cloud cover contamination and can capture the snow information with an all-time and all-weather ability. Since microwaves penetrate most of the snow cover, it is possible to detect the SD and the SWE. The typical products include the Scanning Multichannel Microwave Radiometer (SMMR) SD product (Chang et al., 1987), the Special Sensor Microwave/Imager (SSM/I) SD product (Grody and Basist, 1996), the Microwave Radiation Imager (MWRI) SD product (Che et al., 2016), and the Advanced Microwave Scanning Radiometer – Earth Observing System (AMSR-E) – SWE product (Gao et al., 2010a; Anthony et al., 2008). Microwave-based products, which include both active and passive forms, have a high temporal resolution and can quickly cover the Earth's surface in 3 to 5 d. Compared to the passive microwave platform, active microwave remote sensing has a higher spatial resolution (<1 km) and can provide the details of snow cover. However, active microwave remote sensing is rarely used in research into the SD because it features a highly transparent frequency to snow (Baghdadi et al., 1997) and can only recognize wet snow reliably (Dietz et al., 2012a). In contrast, passive microwave remote sensing is effective in retrieving both the SD and SWE. However, the spatial resolution of passive microwave remote sensing systems is usually so low (>1 km) that they are not able to acquire the detailed information of snow (Foster et al., 1984). Generally speaking, passive microwave-based products are not suitable for mesoscale or micro-scale research and may result in confusion between wet snow, thin snow, and forest (Tait et al., 2000). In addition, because of the limitations of the imaging orbit gap, microwave-based products are subject to spatial gaps.

Compared to the microwave-based products, optically based products have no spatial gaps resulting from the imaging orbit gap. Optically based products are derived according to the differences in visible and infrared spectra between snow and cloud, bare land, and vegetation. The representative products are mainly derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard the Aqua and Terra satellites (Hall et al., 2002), AVHRR data (Simpson et al., 1998), VEGETATION data (Xiao et al., 2004), Thematic Mapper (TM) data (Rosenthal and Dozier, 1996), Enhanced Thematic Mapper Plus (ETM+) data (McFadden et al., 2011), Operational Land Imager (OLI) data (Crawford, 2015), and Visible and Infrared Radiometer (VIRR) data (Zhang et al., 2017). Moreover, the weekly snow cover product (Brown et al., 2010) from NOAA is also derived from optical remote sensing imagery. The optically based products have formed long observation time series (from 1960) and have a high spatial resolution (≤1 km). The unique spectral signature (e.g., high albedo or high thermal emissivity) of snow derived from optical observations is commonly used in the snow cover mapping algorithms (Simic et al., 2004). Optical observations can also be used to extract useful snow information. For example, the snow cover area (SCA) or snow cover extent (SCE) is usually extracted by binary classification or the use of the subpixel snow-fraction algorithm with the optical observations (Rittger et al., 2013; Salomonson and Appel, 2004). The optically based products can not only be applied for monitoring the spatio-temporal variation of snow cover but can also be set as the input of atmospheric forecast and hydrological models. As a result, the optically based products are widely used in practical applications. However, they are inevitably subject to cloud cover contamination. Especially in the snow period, the cloud fraction of snow product is usually more than 10 % (Liang et al., 2008b).

In order to derive cloud-free snow cover products, combined products are generated, which are combinations of satellite data (optical and microwave), climate station observations, and models (Frei et al., 2012; Tait et al., 2000). For example, the Canadian Meteorological Centre (CMC) snow product (Drusch et al., 2004) is a combination of station observations and models, and GlobSnow (Metsämäki et al., 2015), from the European Space Agency (ESA), is a multiple-dataset snow cover product generated from satellite data, station observations, and models. In addition, the Interactive Multisensor Snow and Ice Mapping System (IMS) produced by the US National Ice Center (NIC; Ramsay, 1998) is the fusion of many kinds of optical and microwave data, including Advanced Very High Resolution Radiometer (AVHRR) data, SSM/I data, AMSR-E data, Geostationary Operational Environmental Satellite (GOES) data, Polar Operational Environmental Satellite (POES) data, european geostationary meteorological satellite (METEOSAT) data, Japanese Geostationary Meteorological Satellite (GMS) data, the National Centers for Environmental Prediction (NCEP) model data, US Air Force (USAF) snow and ice analysis data, and so on. The IMS is also jointly supported by the US National Oceanic and Atmospheric Administration (NOAA), the US Navy, and the US Coast Guard. Regardless of cloud cover, the IMS produces near-real-time products with spatial resolutions of ∼1, ∼4, and ∼24 km, which provide the input for atmospheric forecast models (Brown et al., 2014). In addition, the Air Force Weather Agency (AFWA) and National Aeronautics and Space Administration (NASA) snow algorithm (ANSA) blends AMSR-E, MODIS, and Quick Scatterometer (QuikSCAT) data products (Foster et al., 2011). Since ANSA integrates optical, passive, and active microwave data, it can map the snow-covered area (SCA), the fractional snow cover (FSC), the SWE, and the snowmelt area. The combined products draw together the respective advantages of each of the component products to improve the accuracy, quality, and spatio-temporal continuity of snow cover. Several problems encountered when the component products are used alone, including cloud cover and low accuracy, have been solved. The combined products provide more information on the state of snow cover than each of the component products. However, the primary disadvantage of the combined products is the poor spatial resolution (>1 km; Tait et al., 2000).

Although the combined products have the advantage of spatio-temporal continuity, they need many different data sources, and their spatial resolution is restricted to the lowest spatial resolution among the data sources. As a result, many attempts have been made to derive spatio-temporally continuous snow cover products from optical remote sensing data. Among the existing optically based snow cover products (e.g., AVHRR, VEGETATION, and VIRR), the MODIS products have become one of the main data sources for ice and snow research. The MODIS products have advantages in the spatial and temporal resolutions, global coverage, long time series, and free access, which together allow real-time, accurate, and large-area snow cover variation monitoring. The MODIS products are available through the National Snow and Ice Data Center Distributed Active Archive Center (NSIDC DAAC). The snow cover products of MODIS are derived from the SNOMAP algorithm (Hall et al., 1995). SNOMAP automatically uses the normalized difference snow index (NDSI) and decision strategies to identify snow cover (Hall et al., 2002), which makes this type of snow product more consistent with station observation data and other higher-spatial-resolution remote sensing data (Parajka and Blöschl, 2006; Wang et al., 2008; Bitner et al., 2002; Klein and Barnett, 2003; Chelamallu et al., 2014; Huang et al., 2011). According to accuracy assessments, under clear-sky conditions, the overall accuracy (OA) of MODIS snow cover products ranges from 85 % to 99 % (Parajka et al., 2012), and the overall absolute accuracy is ∼93 % (Hall and Riggs, 2007). Version 5 of MODIS snow cover products provides a binary product and a fractional snow product (Rittger et al., 2013; Li et al., 2018). The binary product is the SCA and the fractional snow product is the FSC, which are both derived from the NDSI (Singh et al., 2014; Liang et al., 2017). The SCA can be used to estimate the regional snow resource when combined with the SWE, and the FSC is a very important parameter for estimating the SWE (Molotch and Margulis, 2008). The MODIS products have temporal resolutions of 1 to 8 d and 1 month and spatial resolutions from 500 m to 0.05^∘ (Déry et al., 2005). Like other optically based products, the MODIS products are prone to contamination by a large percentage of cloud cover (Wang et al., 2009; Maskey et al., 2011), which limits the continuity of snow cover monitoring in space and time. Cloud cover also causes confusion in snow discrimination (Ault et al., 2006). As a result, it is of importance to remove the cloud coverage from the MODIS snow cover products.

In the past decades, there have been a large number of algorithms developed to improve the spatio-temporal continuity of the MODIS SCA products. The SCA products are usually just classification data, so the cloud removal method is notably different from the methods for traditional remote sensing images (Li et al., 2016; Shen et al., 2015; X. Li et al., 2014). In this review, we comment on the recent developments in producing cloud-free MODIS snow cover products (without any special instructions, the product is MODIS SCA). The algorithms for generating cloud-free MODIS snow cover products are mainly categorized into temporal methods, spatial methods, spatio-temporal methods, and multi-source fusion methods.

The remainder parts of this paper are arranged as follows. Section 2 introduces the spatial methods for generating spatio-temporally continuous snow cover products. The temporal methods, spatio-temporal methods, and multi-source fusion methods are then introduced in Sects. 3–5, respectively. Section 6 summarizes the current validations and evaluations of the cloud removal methods for MODIS snow cover products. Section 7 exposes the future direction. Finally, the conclusions are provided in Sect. 8.

2 Spatial methods

The spatial methods remove the cloud cover of the snow product based on the spatial patterns of the snowpack. The main spatial methods are the spatial filter (SF), the snow line mapping (SNOWL) approach, and the locally weighted logistic regression (LWLR).

2.1 The spatial filter (SF)

The most common spatial method is the SF (Gafurov and Bárdossy, 2009; Parajka and Blöschl, 2008; Paudel and Andersen, 2011), where the idea is to replace a cloudy pixel using the four or eight neighboring non-cloud pixels. There have been many rules put forward to replace the cloudy pixel, as follows: (1) the cloud pixel should be reassigned as a snow pixel on the condition that three of its four direct “side-bordering” neighboring pixels are snow (Paudel and Andersen, 2011; Lindsay et al., 2015). (2) The cloudy pixel is replaced by the main classification (land or snow); i.e., the class of the majority of the non-cloud pixels in a neighborhood is used to replace the cloudy pixel (when there is a tie, the cloudy pixel is replaced by snow; Parajka and Blöschl, 2008; Tong et al., 2009a, b). (3) Gafurov and Bárdossy (2009) proposed that when the eight neighboring pixels with a lower elevation than the cloudy pixel are covered by snow, the cloudy pixel should be labeled as being snow covered. (4) López-Burgos et al. (2013) replaced the cloudy pixel based on the elevation and aspect; that is, if any neighboring pixel has snow cover and has the same aspect and a lower elevation, then the cloudy pixel is classified as snow. In some cases, the same cloudy pixel may be labeled as snow or no snow by applying different rules. Therefore, the most suitable rule should be chosen according to the regional snow cover change rule.

The SF is usually effective for small-area cloud cover, and cloud with a proportion of no more than 10 % can be removed. For example, the SF in an eight-neighborhood pixel decreases the cloud coverage of MODIS products by 7 %, and the decrease in OA is just 0.7 % (Parajka and Blöschl, 2008). According to practical application, the SF is not very sensitive to the size of filter window (Zhou et al., 2005). Additionally, in mountainous regions, the elevation is assumed to be the dominant factor affecting the snow cover distribution. Due to the complicated topography, the accuracy of the SF usually declines with rising elevation (Tong et al., 2009b).

2.2 The snow line mapping approach (SNOWL)

The SNOWL method (Parajka et al., 2010; Dietz et al., 2013), which is also called the snow transition elevation method (Gafurov and Bárdossy, 2009; Shea et al., 2013), reclassifies the cloudy pixels as snow or land according to the elevation distribution characteristics of the snowpack. The cornerstone of SNOWL is the snow line and the land line. The snow line is the snow-covered elevation where all pixels above it are covered by snow, and the land line is the minimum elevation where snow exists (Krajčí et al., 2014, 2016; Lei et al., 2012). On that account, all the cloudy pixels above the snow line are labeled as snow, and all the cloudy pixels below the land line are labeled as land by SNOWL. The cloudy pixels between the snow line and land line are labeled as partial snow. However, the partial snow brings some uncertainties to the monitoring of snow cover variation.

SNOWL is relatively simple and easy to implement and performs well in both high and low elevations. In order to make sure that the snow line and land line are accurately labeled, the product should be at least 70 % cloud free (Gafurov and Bárdossy, 2009). Generally speaking, as the cloud fraction increases, the accuracy of cloud removal decreases. Some scholars considered aspect, topography, and land cover classes to improve the accuracy of SNOWL (Paudel and Andersen, 2011; Da Ronco and De Michele, 2014a). In a few special cases, the snow line for the whole area is hard to find, so the regional snow line for the local area is labeled (Parajka et al., 2010). In the Trans-Himalayan region, the improved regional SNOWL method removes 38.28 % of the cloud, and the misclassification error is less than 2 % (Paudel and Andersen, 2011).

2.3 Locally weighted logistic regression (LWLR)

For the LWLR method (López-Burgos et al., 2013), the snow cover probability of the cloudy pixel is predicted by the topographic and spatial relations with its neighboring cloudless pixels. LWLR uses two explanatory variables of elevation and aspect for snow occurrence probability. The information of the neighboring pixels is inversely weighted with distance and is fit to a logistic curve. The pixels closer to the cloudy pixel are assigned with a heavier weight than those that are farther apart. Finally, the estimated snow probability by LWLR is converted to a binary result according to a selected threshold. However, the method for choosing a better threshold to obtain a binary result needs more tests, and its high cost (>20 h in the Salt River basin in Arizona) limits its application to some degree. Experiments demonstrated that the LWLR method could reduce the cloud obscuration by 93.8 % in the Salt River basin in Arizona (López-Burgos et al., 2013).

Among all the spatial methods, the computational complexity of LWLR is the highest. The spatial methods mainly depend on the spatial patterns of the snowpack to reclassify the cloudy pixels. For the majority of the spatial methods, the core idea is to utilize neighboring cloud-free pixels. However, when the cloud fraction is high, the accuracy will decrease.

3 Temporal methods

The temporal methods reclassify the cloudy pixels according to the temporal correlation and change rule of the snow cover. According to the time span of the product used, the temporal methods are mainly classified as the Terra and Aqua combined (TAC) method and temporal filters. The reason why the TAC method is considered to be a kind of temporal method, in our opinion, other than the multi-source fusion, is that the Aqua and Terra satellites are both equipped with MODIS (of nearly the same design). Their combination amounts to a composite of MODIS at different times.

3.1 Terra and Aqua combined (TAC)

Among the temporal methods, the TAC method is the simplest and most straightforward method. Since cloud coverage is very variable in a few hours, the TAC method (Parajka and Blöschl, 2008; Xie et al., 2009; Wang and Xie, 2009; Mazari et al., 2013) has become the most popular temporal method. The TAC method can decrease the cloud coverage ratio by 5 %–20 % without significantly sacrificing accuracy of the products. This method blends the same-day MODIS snow products on a pixel basis. If a pixel is cloudy in one product and cloud-free in another product, the cloudy pixel will be updated by the classification of the cloud-free pixel. The combination usually has the following priority scheme: snow > no snow > cloud.

The basic assumption of the TAC method is that complete snowmelt and snowfall did not occur during the time interval. In the process of merging the two products, they are considered to be identical. In fact, the Terra and Aqua products still have some small differences, for the following reasons. The first is that they are acquired at different times (3 h). The second is that most of the detectors in Aqua MODIS band 6 failed (Shen et al., 2014). In the early days (before 2012), the snow mapping method for Aqua used band 7 as a substitute for band 6; later on, Aqua MODIS band 6 was restored by the quantitative image restoration method (Gladkova et al., 2012) with a high degree of accuracy. Even so, the two products still have slight differences. According to the ground observations, the TAC method inherits the reclassification errors of the input data of Terra and Aqua (Xia et al., 2012). Gao et al. (2010b) confirmed that the TAC method can reduce the cloud cover by 5 %–14 % and 8 %–12 % at the yearly and monthly scale, respectively, and the OA is 89.7 %, which is lower than MOD10A1 by 0.7 % and higher than MYD10A1 by 1.4 %.

3.2 Temporal filters

Another popular kind of temporal method is the temporal filter methods (Parajka et al., 2012; Hori et al., 2017), which mainly include adjacent temporal deduction (ATD), multi-day combination (MDC), season filter (SFil), and temporal interpolation using a mathematical function. The first three kinds of methods are applied to the SCA product (as shown in Fig. 2), and the last method is suitable for the FSC product. It should be noted that temporal filters are also referred to as temporal interpolation methods in some literature (López-Burgos et al., 2013; Hüsler et al., 2014).

Figure 2 Temporal filters for the MODIS SCA products.

3.2.1 Adjacent temporal deduction (ATD)

ATD (Paudel and Andersen, 2011; Lindsay et al., 2015; Dietz et al., 2013) is an effective way to deduce the surface conditions via the same pixel in the previous and subsequent days without reducing the spatial and temporal resolution. It is assumed that if the preceding and the following day of the cloudy pixel remain the same (land or snow), the condition of the cloudy pixel will remain unchanged (Gafurov and Bárdossy, 2009). When the previous and the next day are different, the cloudy pixel is still labeled as cloud. In fact, ATD uses ±1 d information and is similar to the works of Li et al. (2008) and Gafurov and Bárdossy (2009), who used ±1 or ±2 d information. However, ATD is different from direct temporal replacement (Zhao and Fernandes, 2009), with the sacrifice of the 1 d temporal resolution.

Owing to the useful information from the previous day and the subsequent day, the accuracy of cloudy pixel reclassification is high. On the basis of the TAC method, ATD can not only decrease cloud fraction by 25 % but can also achieve an accuracy of 96.3 % (Gafurov and Bárdossy, 2009). However, since the snow cover is assumed to remain unchanged in the given temporal interval, the accuracy in snow-transitional periods is obviously lower than that in snow-stable periods (Gao et al., 2010b). In other words, ATD is not suitable for a context with variable snow covers.

3.2.2 Multi-day combination (MDC)

Among the temporal filter methods, MDC (Parajka and Blöschl, 2008; Dietz et al., 2012b; Zhang et al., 2012; Gao et al., 2011a) is the most widely used method. The cloudy pixels are replaced by the cloudless pixels in a temporal interval ranging from 1 d to more than 10 d with a constant or flexible way (Chen et al., 2014). For example, given a temporal window of 10 d, the constant way of MDC means that the combination is implemented in 10 d, and the flexible way represents how the combination can be implemented in varying days (≤10 d). To some degree, the above-mentioned ATD can be regarded as a special case of MDC. As far as MDC is concerned, the temporal window can be preceding (Wang et al., 2014) or preceding and backward (Sharma et al., 2014; Foppa and Seiz, 2012; Coll and Li, 2018). Usually, priority is given to the preceding days (Marchane et al., 2015). For the preceding temporal window, it amounts to a backward temporal filter in which the cloudy pixel is replaced by the previous cloud-free pixel in turn. The cloud-gap-filled method (Hall et al., 2010) also belongs to this scenario. When the land cover is considered to be unchanged, the information of the next 2 d and the previous 2 d can be used (Gafurov and Bárdossy, 2009; Da Ronco and De Michele, 2014a). To improve the accuracy of the cloud removal, certain thresholds (e.g., cloud percentage and the number of composite days) should be imposed in MDC (Gao et al., 2010b).

MDC can reduce a high fraction of cloud coverage. On the one hand, as the temporal window increases, the temporal resolution and accuracy of the result will decrease (blur). For example, it has been demonstrated that the accuracies of the 2, 4, 6, and 8 d combinations are 89.5 %, 89.0 %, 88.2 %, and 87.8 %, respectively (Gao et al., 2010b). On the other hand, the remaining cloud cover will decrease with the increase of the temporal window span. Hence, it is a trade-off between the remaining cloud and the blurring of the snow information (Hüsler et al., 2014). The longer temporal window also corresponds to a larger area of snow cover. Thus, a balance between temporal resolution and SCA should be considered. When cloud covers the pixel for more than the temporal window, MDC does not work. In this situation, the remaining cloudy pixels should be processed by other methods.

3.2.3 Season filter (SFil)

Among the temporal methods, SFil (Gafurov and Bárdossy, 2009; Gafurov et al., 2013; Lindsay et al., 2015) uses the longest time-series information to reclassify the cloudy pixels. This method needs two thresholds in one unbroken hydrological year: a complete snowmelt day and snow accumulation start day. As shown in Fig. 3, the complete snowmelt day is the day when the pixel is no longer covered by snow, and the snow accumulation start day is the start day when the snow accumulates. The hydrological year is divided into a land season and a snow season by SFil (Da Ronco and De Michele, 2014a). For example, the cloudy pixels before the snow accumulation start day and after the complete snowmelt day are in the land season, and they will be reclassified as land.

SFil can remove all the cloud cover; however, it does not take the phenomenon of more than one snow cycle occurring in 1 hydrological year into consideration. Thus, three thresholds are introduced in each hydrological year (Paudel and Andersen, 2011): the maximum snow extent day, the minimum snow extent day, and snow accumulation start day. In this way, the improved SFil effectively increases the accuracy. Similarly, Lindsay et al. (2015) defined two snow seasons: the full snow season and the continuous snow season. The full snow season is the period between the first day and the last day of snow cover, and the continuous snow season has snow cover for at least 14 d, with intervening snow-free periods of no more than 2 d. Since the short-term snowfall or snowmelt is not considered by SFil, it is often applied after other spatial or temporal methods. On the whole, SFil obtains a slightly lower accuracy than other temporal filter methods.

Figure 3 Snow cover depletion curve and the threshold days.

3.2.4 Temporal interpolation using a mathematical function

As stated previously, ATD, MDC, and SFil are applied to the binary MODIS product, and very little attention has been paid to the FSC product. For example, Tang et al. (2013) filled the cloud-contaminated pixels by cubic spline interpolation with the cloud-free observation pixels, and the mean absolute error was less than 0.1. This method is based on the the features of cloud that are rapidly changing and shift daily. In order to distinguish it from the temporal interpolation of SCA products, the temporal interpolation for FSC products is called temporal interpolation using a mathematical function, which involves interpolating the cloudy information of the FSC product along the temporal dimension of the same pixel (Tang et al., 2013, 2017; Xu et al., 2017).

In summary, the temporal methods make use of the cloud instability and the snow correlation in time. They can effectively reduce the cloud cover, partly or completely, and the accuracy is high. ATD and MDC (with time series that are not long enough) cannot reduce the cloud cover completely and may neglect short snowfall events. SFil can remove the cloud cover completely. However, when the cloud cover exists persistently in a region, the accuracy of the temporal methods will decrease. For the FSC products, temporal interpolation using a mathematical function is effective (Tang et al., 2017).

4 Spatio-temporal methods

Although the spatial methods and temporal methods remove the cloud with a high accuracy, the majority of them cannot reduce the cloud completely. In order to minimize the extent of the cloud cover, many scholars have come up with spatio-temporal methods to use the complementary advantages of the temporal methods and spatial methods. This type of method relies on the correlations of snow cover in space and time with two basic forms. One is to utilize the spatial and temporal information step by step, usually as a multi-step combination method (Parajka and Blöschl, 2008; Da Ronco and De Michele, 2014b; Gurung et al., 2011; Zhou et al., 2013; Şorman et al., 2019). The other is to utilize the spatial and temporal information simultaneously (Li et al., 2017; Xia et al., 2012; Huang et al., 2018; Poggio and Gimona, 2015), which we call one-step utilization.

4.1 Multi-step combination

Multi-step combination combines the spatial methods and temporal methods alternately (Dariane et al., 2017). Although the spatial methods or temporal methods are not able to remove cloud absolutely, their successive combination can make a difference. By removing cloud progressively in each step, an accumulated cloud-free result can be obtained by multi-step combination. In other words, attempts at further reducing the residual cloud in the result of the previous step are made in the next step. Among the combined multiple steps, the TAC method is the first step in most cases. SF, SNOWL, and temporal filters are then applied in a variety of ways.

For example, a three-step method of TAC, the SF, and the temporal filter in turn was proposed (Parajka and Blöschl, 2008). Gafurov and Bárdossy (2009) proposed six steps of TAC, MDC, SNOWL, the SF with direct side-border neighboring pixels, the SF with eight neighboring pixels, and SFil. Paudel and Andersen (2011) proposed a five-step approach of TAC, ATD, the SF with four neighboring pixels, SNOWL, and SFil. López-Burgos et al. (2013) proposed a four-step combination of TAC, temporal filter, the SF, and LWLR. Da Ronco and De Michele (2014a) proposed the five-step method of TAC, the conservative temporal filter, SNOWL, the 6-day backward temporal filter, and SFil.

Whatever the combination of spatial methods and temporal methods, multi-step combination independently utilizes snow correlations in space and time, and the result from the previous step directly influences the next step. The combination order of the multiple steps is determined by their respective characteristics. For example, SNOWL requires the cloud fraction to be less than 30 %, which can be satisfied after the use of the TAC method and other temporal methods. As a result, SNOWL is often applied immediately after these methods. No matter what the combination, multi-step combination needs to adopt a feasible strategy based on the topographic features, temporal variation, and spatial heterogeneity of the snow cover. At the same time, a trade-off exists between the accuracy and the cloud fraction for the various steps. These simple multi-step combination methods have been proven to be effective and efficient in cloud reduction and agree well with station observations (Paudel and Andersen, 2011). However, this independent and successive utilization cannot take full consideration of the spatio-temporal information.

4.2 One-step utilization

In contrast, one-step utilization utilizes the spatio-temporal information simultaneously, rather than step by step. As we know, to date, most attention has been paid to the multi-step combination methods. It is only in the last few years that efforts have been made with one-step utilization methods. The one-step utilization methods are introduced in the following.

Xia et al. (2012) first introduced variational interpolation (Shen and Zhang, 2009) from the image processing field to construct a three-dimensional implicit function with five consecutive daily images (space–time manifold), which has an advantage over representing the complicated surfaces in high-dimensional spaces. The shape of the snow boundary can be easily obtained by the implicit function. It is a good case of utilizing the snow cover evolution with continuity in space and time. The results indicate that variational interpolation can maintain a close accuracy to the original product. However, its computational efficiency needs to be improved, since it intends to retrieve the space–time surface of snow cover on consecutive days.

The cloud-filling method (Poggio and Gimona, 2015) is a hybrid of the generalized additive model (GAM) and the geostatistical space–time model. The multi-dimensional spatio-temporal GAM models the binary variables, and geostatistical kriging accounts for the spatial details. The space–time correlations of snow cover are well exploited by the cloud-filling method. Even for a high fraction of cloud cover, the cloud-filling method still achieves a satisfactory accuracy and is suitable for the seasonal variation of snow cover. However, the requirement for a lot of ancillary data (e.g., land surface temperature, land cover, and soil pattern data) limits its application to some degree.

The adaptive spatio-temporal weighted method (ASTWM; Li et al., 2017) estimates the cloudy pixel according to the probability of snow cover, which is the adaptively weighted combination of the snow probabilities in space and time. Experiments demonstrated that the ASTWM not only removes the cloud completely but also achieves a high OA of above 93 % under different cloud fractions. However, the ASTWM resorts to the optimal weight of spatial and temporal probability with a high cost, and it may be able to reclassify snow cover as not being snow cover under darker conditions (Krajčí et al., 2014).

Combining the spatio-temporal–spectral information and environmental relation, Huang et al. (2018) proposed a spatio-temporal model based on the hidden Markov random field (HMRF). The information of the cubic spatio-temporal neighborhood is effectively utilized by the HMRF-based method. The snow mapping accuracy of HMRF-based method can reach 88 % and decrease the cloud cover to 1 %, which improves the overall snow mapping accuracy and reduces the omission error of the original product. For the challenging task of snow cover mapping in the transition periods of snow and in forest areas, it also obtains promising results.

Recently, increasing spatio-temporal methods have been proposed to alleviate cloud cover. In fact, multi-step combination has a numerical advantage over one-step utilization. Multi-step combination usually makes an assurance of removing the whole cloud by successive multiple steps. In contrast, one-step utilization achieves this goal in one step. One-step utilization jointly utilizes the spatio-temporal information from the snow coverage, which contributes to more promising cloud removal results.

5 Multi-source fusion methods

The above-mentioned three types of methods mainly use the spatial and temporal information of the snow cover from the same optical remote sensing sensor. In contrast, multi-source fusion methods (Romanov et al., 2000; Gao et al., 2010a; Yu et al., 2012; Gafurov et al., 2015; Wang et al., 2015; Dong and Menzel, 2016) utilize the complementary information among different sources (Shen et al., 2016), such as optical observations, microwave observations, and station observations.

5.1 Optical and microwave observations

In general, optically based products have a high spatial resolution and are influenced by cloud cover, whereas microwave-based products have a low resolution and a good cloud-penetrating capacity. Therefore, the fusion of optical and microwave data has been the most representative multi-source fusion method of cloud removal, e.g., MODIS and AMSR-E (Liang et al., 2008a; Gao et al., 2011b; Akyurek et al., 2010; Huang et al., 2014, 2016; Deng et al., 2015; H. Li et al., 2014; Bergeron et al., 2014), GOES and SSM/I (Romanov et al., 2000), and the visible/infrared spin–scan radiometer (VISSR) and microwave radiation imager (MWRI; Yang et al., 2014).

The fusion of MODIS and AMSR-E is the most frequently used method. For the AMSR-E SWE product, a value of zero represents a land surface, and values of 1–240 mm represent a snow-covered surface. The cloudy pixels of the MODIS product are reclassified as land if SWE = 0 or as snow if SWE ranges from 1 to 240 (Wang et al., 2015). For more fusion rules, please refer to Liang et al. (2008a). After recoding the MODIS and AMSR-E products, some scholars have adopted a rule of the maximum value composite (Yu et al., 2012; Wang et al., 2018). It is noted that the spatial resolution (25 km) of the AMSR-E SWE should be resampled to the 500 m of the MODIS snow cover product before fusion. In addition, for the AMSR-E SWE, the gaps with daily changing locations are filled by the ATD method (Wang et al., 2015). The fusion accuracy is much higher than the MODIS daily and 8 d products (Liang et al., 2008a). However, the fusion decreases the spatial resolution and results in uncertainties because of the low spatial resolution of AMSR-E (Gao et al., 2010a, b).

In addition, to remove the cloud cover in the MODIS product, Yu et al. (2016) developed the method of fusing MODIS and the combined IMS product, where the cloudy pixels of MODIS are replaced by the values of the IMS pixels. As stated previously, the IMS product itself is a combination of optical and microwave data. Compared with the AMSR-E spatial resolution of 25 km, the cloud-free IMS product has a much higher spatial resolution of ∼1 and ∼4 km, which is also resampled to the same resolution as MODIS. Therefore, the fusion of MODIS and the IMS product has an advantage over the fusion of MODIS and AMSR-E. Moreover, Wang et al. (2018) proposed fusing MODIS and the AMSR-E daily with the SWE product or IMS product. They pointed out that the accuracy of the cloudy pixel reclassification is dependent on AMSR-E and IMS.

5.2 Optical and station observations

The observations of the existing meteorological stations are long-term, high-precision, and point-based observations. In contrast, remote sensing observations are spatially continuous. Based on the correlations between station observations and spatial snow cover patterns, the snow cover can be reconstructed by fusing station data and high-spatial-resolution remote sensing data. For example, based on station observations and optical observation data, the conditional probability of each cloudy pixel being snow can be calculated to reclassify the residual cloudy pixels (Dong and Menzel, 2016; Gafurov et al., 2015, 2016). The conditional probability represents the probability of a pixel being snow cover, on the condition of the SD being higher than zero at the station. The conditional probability also implies the similarity between different meteorological stations. In other words, the snow presence in one station is predicted by the information about the presence of snow at another station. This method has been confirmed as being effective, especially during the snow season. In the Zerafshan River basin, the accuracy is only slightly lower than the original MODIS product (Gafurov et al., 2015), according to Landsat-derived snow cover (Gafurov et al., 2013). However, to some degree, the distribution and number of meteorological stations limit the predictive ability of snow cover reconstruction.

Gafurov et al. (2016) developed an all-in-one software package called MODSNOW-Tool, with advanced cloud removal algorithms for MODIS snow cover products. The integrated algorithms include the six-step method in Gafurov and Bárdossy (2009) and the conditional probability method in Gafurov et al. (2015). MODSNOW-Tool is equipped with operational and non-operational modes, which consist of seven processing modules. The operational mode generates a daily cloudless snow cover map without user interaction, and the non-operational mode generates a historical snow cover map. This tool can remove the complete cloud cover, which is a major breakthrough.

5.3 Optical, microwave, and station observations

In terms of multi-source fusion methods, the fusion of optical and microwave data is the most common approach, and the fusion of optical and station observations has attracted relatively little attention. Moderate attention has been paid to the fusion of optical, microwave, and station observations. Brown et al. (2010) utilized 10 data sources, including optical, microwave, and station observations, to estimate the cloud-free snow cover. Before the fusion of the 10 data sources, the consistency of each dataset was assessed by their correlations. The fusion result was then obtained by converting the average anomaly series to the first differences then joining the difference series, in which the average anomaly series was computed from each reference period (Brown et al., 2010). This method has been well validated for monitoring the SCE variation in the Arctic region (Brown et al., 2010).

6 Validation and evaluation

6.1 In situ and remote sensing data based evaluation

The evaluation of the accuracy and effectiveness of cloud removal is also very important to the cloud-removed MODIS snow cover products. When the SD data of climate stations are available, a time series of in situ observations can be used to validate the temporal effectiveness of the cloud-removed products in real-data experiments. The SD data of the in situ observations are usually set as the standard data. The nearest pixel to the site is classified as snow when the SD exceeds a threshold value; otherwise, it is classified as no snow (Parajka and Blöschl, 2008). Through the comparison of the in situ SD and the reclassified snow cover product, the effect can be evaluated.

In situ observation based evaluation is a direct and valid method. However, in most cases, because of data privacy policies or the absence of climate stations, researchers cannot acquire in situ observations of snow cover. To conduct the validation, the majority of researchers resort to evaluation based on remote sensing data via simulated experiments, following the work of Gafurov and Bárdossy (2009). Firstly, the Aqua and Terra snow cover products were combined by the TAC method. A number of images with a cloud fraction of less than 10 % were selected as “truth” products, and the cloud masks of the other images with a larger cloud fraction were applied to cover the truth products and get the “observation” products. Next, the observation products were reclassified. Finally, the results were compared with the truth products. In addition, some other SCA products with a higher resolution can also be considered as the true ground data. For example, the Landsat TM/ETM+/OLI FSC products have been used to validate the MODIS FSC products (Liang et al., 2017; Crawford, 2015; Kuter et al., 2018). However, when the MODIS snow products are covered by cloud, the optically based products with a higher spatial resolution will also be covered by cloud. This method is thus only effective for the validation of clear-sky accuracy.

6.2 Quantitative evaluation indicators

There have been many kinds of indicators used to evaluate the result of cloud-removed MODIS snow products. Usually, the effectiveness is described by cloud fraction, while the accuracy is evaluated by the OA, overall clear-sky accuracy (OC), overestimation error (OE), and underestimation error (UE; Gafurov and Bárdossy, 2009), which are calculated by the following expressions:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{OA}={\displaystyle \frac{N_{\mathrm{s}}^{\mathrm{s}}+N_{\mathrm{ns}}^{\mathrm{ns}}}{N_{\mathrm{a}}}},\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{OC}={\displaystyle \frac{N_{\mathrm{s}}^{\mathrm{s}}+N_{\mathrm{ns}}^{\mathrm{ns}}}{N_{\mathrm{c}}}},\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{OE}={\displaystyle \frac{N_{\mathrm{s}}^{\mathrm{ns}}}{N_{\mathrm{c}}}},\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{UE}={\displaystyle \frac{N_{\mathrm{ns}}^{\mathrm{s}}}{N_{\mathrm{c}}}},\end{array}$$

where N_a denotes the amount of cloudy pixels from the observation product, N_c denotes the amount of cloudy pixels in the observation product, $N_{\mathrm{s}}^{\mathrm{s}}$ denotes the amount of cloudy pixels which are reclassified as snow pixels in the observation product and are snow pixels in the true product, $N_{\mathrm{ns}}^{\mathrm{ns}}$ denotes the amount of cloudy pixels which are reclassified as not being snow pixels in the observation product and are not snow pixels in the true product, $N_{\mathrm{s}}^{\mathrm{ns}}$ denotes the amount of cloudy pixels which are reclassified as snow pixels in the observation product and are not snow pixels in the true product, and $N_{\mathrm{ns}}^{\mathrm{s}}$ denotes the amount of cloudy pixels which are reclassified as not being snow pixels in the observation product and are snow pixels in the true product.

7 Future directions

Firstly, multi-source fusion is still a promising direction for the cloud removal of the MODIS snow cover product. Although optical observations, microwave observations, and station observations have been fused to derive the spatially seamless snow cover product, the complementary information of the multi-source observations is not utilized best. For example, the microwave-based snow cover product is just used to replace the optically based snow cover product under cloud coverage simply. On one hand, the correlations of heterogeneous products are not well modeled for a better use. On the other hand, most of the existing fusion is to replace MODIS with the coarser AMSR-E; the difference between the spatial resolution of microwave-based and optically based products is very great, which is usually neglected in the replacement process. Additionally, the station observations, which are sparsely distributed in space, usually provide the temporal variation rule of snow cover. Dong and Menzel (2016) provided an available approach to fuse the optical and station observation with the conditional probability interpolation. However, the overestimation error of the fused snow cover products is still high. In the future, the mathematical relation among the optical observations, microwave observations, and station observations should be modeled comprehensively.

Secondly, the cloud removal algorithms of MODIS snow cover products will be a benefit of the platforms with higher spatial resolution more easily and frequently. With the development of the technology, the spatial resolution of the sensor becomes higher and higher. In the framework of the multi-source fusion, the microwave-based observation with a higher spatial resolution than AMSR-E should make a difference, especially with the Sentinel series. For example, Sentinel-1 SAR has the spatial resolution of 20 m (Snapir et al., 2019; Nagler et al., 2016), which will significantly improve the fusion accuracy of the MODIS snow cover product. Additionally, the optical observations of Sentinel series, e.g., the Sentinel-2 Multispectral Instrument (MSI) and Sentinel-3 Sea Land Surface Temperature Radiometer (SLSTR; Nagler et al., 2018; Zhu et al., 2015), also have the potential to provide a snow cover product with higher spatial resolution in the future. However, the spatial resolution of Sentinel series is higher than MODIS, which results in the problem with a smaller image swath and a longer revisit period. In addition, the high-spatial-resolution data will not only contribute to snow mapping, e.g., unmanned aerial vehicles (UAVs) acting as an effective supplement for snow mapping (Liang et al., 2017), but they will also play a significant part in the accuracy validation of the cloud-removed MODIS snow cover products in the near future.

Thirdly, the new algorithms for MODIS Collection 6 (C6) products should be developed correspondingly. As we know, the most common algorithms of cloud removal for MODIS snow products have been aimed at the binary product (V005). In the spring of 2016, MODIS C6 products were published (Malmros et al., 2018). In the MODIS C6 products, the binary SCA products have been substituted by the NDSI, and the FSC product is not supplied at all (Hall and Riggs, 2016a, b; Riggs et al., 2017). Research work has demonstrated that the has a strong spatial and temporal agreement with Landsat TM/ETM+ (Kuter et al., 2018). Increasing research is thus moving toward MODIS C6 data (Malmros et al., 2018; Huang et al., 2018). On the one hand, the aforementioned spatial methods, temporal methods, spatio-temporal methods, and multi-source fusion methods will still work well, on the condition that the newly released MODIS C6 data are converted into a binary SCA product by the MODIS snow mapping algorithm SNOMAP. On the other hand, the cloud removal algorithms for the NDSI products should also be paid attention to. For the snow product transformation from binary product to NDSI, the cloud removal methods should be changed accordingly. The NDSI has the property of seasonal periodicity and is similar to a vegetation index to some degree, e.g., the normalized difference vegetation index (NDVI) or enhanced vegetation index (EVI). In this sense, the methods of reconstructing a gap-free vegetation index (Yang et al., 2015; Chen et al., 2004; Lovell and Graetz, 2001; Roerink et al., 2000; Poggio et al., 2012) will provide a reference for the cloud-free NDSI. Without doubt, the particular features of the NDSI should be considered. In terms of the cloud removal algorithms for the NDSI, machine learning (Tahsin et al., 2017), with its advantages in processing high-dimensional data, will also be a promising direction. Owing to the long time series of observed snow cover, the variation rule of snow cover could be discovered more easily by machine learning.

8 Conclusions

In this paper, the existing methods of generating cloud-free MODIS snow cover products have been summarized from the four aspects of spatial algorithms, temporal algorithms, spatio-temporal algorithms, and multi-source fusion algorithms. These methods utilize the spatial and temporal variation characteristics and the complementary properties of different observation approaches. Thanks to the spatially correlated relations of snow cover, the spatial methods are relatively effective in the removal of neighboring cloudy pixels but are usually ineffective for large-area cloud cover. The temporal methods remove the cloud cover of the products using the temporal variation rule of snow cover. In addition to a high accuracy of cloud removal, the temporal methods have the ability to remove all the cloud, on the condition that the time series is long enough. As their name implies, the spatio-temporal methods take advantage of the spatial methods and temporal methods by successive or one-step utilization of them. The multi-source fusion methods are based on the complementary observations of different types. The fusion of optical, microwave, and station observations contributes to a promising cloud removal result. Algorithms oriented towards the MODIS C6 product will be developed in the near future.