Next Article in Journal
A Bibliometric Profile of the Remote Sensing Open Access Journal Published by MDPI between 2009 and 2018
Previous Article in Journal
SAR Interferometry for Sinkhole Early Warning and Susceptibility Assessment along the Dead Sea, Israel
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Key Laboratory of Remote Sensing of Gansu Province, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou 730000, China
2
Heihe Remote Sensing Experimental Research Station, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou 730000, China
3
Key Laboratory of Western China’s Environmental Systems (Ministry of Education), Lanzhou University, Lanzhou 730000, China
*
Author to whom correspondence should be addressed.
Submission received: 23 November 2018 / Revised: 24 December 2018 / Accepted: 31 December 2018 / Published: 7 January 2019

Abstract

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

Graphical Abstract

1. Introduction

Snow cover is one of the fastest changing objects on the Earth’s surface. It strongly affects the climate system, radiation, and energy balance between the atmosphere and the Earth, hydrological circulation, biogeochemical cycling, and even human activities [1,2]. In the Northern Hemisphere, the snow coverage ranges from approximately 45.2 to 1.9 million square kilometers, which occurs in January and August, respectively. The drastic variation of snow cover spatial extent is a sensitive indicator of climate change [3]. The melting of seasonal snow cover supplies the dominant water source for river runoff and ground water discharge, particularly for arid and semi-arid regions [4]. Approximately 1/6 of the world population and more than 1/5 of the Chinese population depend on water resources supplied by snow cover or glaciers [5]. Snow cover is also a critical parameter in hydrological and ecological process models [6]. There is an urgent requirement to accurately map the spatial distribution of snow cover. This is not only because of a rising demand for fresh water management but also due to the concerns about the effects of climate change.
Several remotely sensed snow cover products provide the long series of global/regional snow cover information [7]. The suite of MODIS snow cover products, available since 2000, has been applied widely in mapping snow cover area because of its moderate spatial resolution and high temporal resolution [8,9]. Extensive studies have demonstrated that the standard MODIS snow products have high snow classification accuracy (>90%) under a clear sky term. The evaluation results of MODIS snow cover product report that it has an average classification accuracy of 93% for MOD10A1 in the Upper Rio Grande Basin [10]; 95% accuracy for eight-days composited MOD10A2 in Northern Xinjiang, Northwest China [11], an accuracy exceeding 90% for MOD10A1 in the Qinghai-Tibet Plateau [12,13], an accuracy of 93% for MOD10A1 in Canada [14], an accuracy of 95% for MOD10A1 in Europe [15] and an accuracy exceeding 92% for MOD10A1 and MYD10A1 in Central Asia [16].
Nevertheless, cloud contamination in MODIS snow cover products has inevitably caused great data gaps which have seriously hindered its application in total environment research. Lin et al. (2013) suggested that the global annual mean cloud cover was approximately 66% from the project of International Satellite Cloud Climatology [17]. Parajka et al. (2006) demonstrated that the clouds have obscured, on average, 63% of MODIS snow cover maps in Austria from 2000 to 2005 [15]. Similar average cloud coverages of approximately 70% in the Quesnel River Basin and 50–60% in Alaska are also reported [18,19]. The mean cloud coverage of MODIS snow cover products was 46.86% (Terra) and 48.47% (Aqua) in Europe from 2000 to 2011 [20], 39% in the Upper Salt River Basin from 1 October 2004 to 31 May 2005 [21]. The mean overall classification accuracy of MODIS snow cover products was approximately 34–50% in all weather conditions [22,23,24]. Wang et al. (2009) verified that the mean annual cloud coverage of MODIS snow cover products was 44% (Terra) and 47% (Aqua) in Northern Xinjiang, China. The overall classification accuracy was below 40% in all-sky conditions [25].
A series of efforts have been made to fill the cloud gaps in the MODIS snow cover maps. The combination of MODIS Terra and MODIS Aqua snow cover products can reduce approximately 10–20% of cloud obscuration and achieve an accuracy of 90% [21,25]. Temporal filtering directly replaces cloud pixels by using ground information of a pixel in the previous few days or/and subsequent days [26]. Previous research has demonstrated that seven days of temporal filtering can eliminate more than 95% of clouds while maintaining an overall accuracy exceeding 92% over Austria [27]. Spatial filtering is designed to replace cloud-contaminated pixels by using information on surrounding adjacent pixels that are not obscured by clouds on the same date, which can fill gaps by approximately 6–13% [23,27,28,29,30]. SNOWL is also a spatial filtering method based on snow transition elevation, which is capable of reducing cloud coverage in standard MODIS snow cover products from 60 to 20% [30,31]. In practice, the combined use of temporal and spatial filtering is the most common way for MODIS cloud mitigation, which can produce cloudless images as well as maintain the accuracy of the source products under clear-sky conditions. The combination of MODIS and passive microwave (e.g., AMSR-E) snow cover products take advantage of MODIS’s high spatial resolution and microwave penetration into the cloud, which is capable of eliminating the cloud-contaminated pixels but has a reduced accuracy due to the low spatial resolution of the passive microwave [23,24,32].
The additional limitation of existing methods often leaves out a percentage of cloud pixels or has low effectiveness when cloud obstruction fills the same area for several continuous days or is widely distributed. Furthermore, these studies were only applicable to MODIS binary snow products, and incapable of delineating the fractional snow cover (SCF) of a pixel, which might be more valuable for hydrological and ecological process models [33]. Dozier et al. (2008) first proposed the notion of treating the time-varying snow cover as a space-time cube. They filled the gaps caused by clouds via cubic spline temporal interpolation at each grid [34]. Tang et al. (2013) also used this cubic spline interpolation algorithm to eliminate the cloud-contaminated pixels in MODIS SCF products of the Tibetan Plateau [12]. This one-dimensional temporal interpolation generates a large number of outliers in the case of continuous cloud coverage whereas the spatial distribution characteristics of snow cover was not considered.
Cheng et al. (2014) developed a new approach based on the idea of similar pixel replacement [35]. The spatio-temporal Markov random fields (STMRF) function was constructed to search for the most appropriate similar non-cloud pixels to replace a cloud-contaminated pixel in MODIS land surface temperature products. Inspired by this work, we presented an innovative gap-filling method, based on the concept of the non-local spatio-temporal filter (NSTF), to eliminate the cloud contamination in daily MODIS SCF products. The ground information of a cloud pixel is estimated through selected similar pixels in the remaining known region. We anticipate that twofold goals can be achieved: (1) reconstructing the ground information of MODIS SCF products obscured by cloud cover and (2) achieving the overall accuracy of MODIS SCF source products under clear a sky term after filling the gaps.

2. Study Area and Data

2.1. Study Area

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

2.2. Data

2.2.1. MODIS SCF Products

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

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

The daily snow depth (SD) recorded by 50 meteorological stations were provided by the Chinese Meteorological Data Sharing Service System (CMDS) (http://cdc.cma.gov.cn/). Most of these meteorological stations are located in the inhabited and low elevation valleys. Only three stations are located at elevations above 2000 m (Figure 1). The SD observations were recorded in centimeters, and the fractional part is rounded up.
The land cover data was extracted from land cover products of China, which was provided by the Cold and Arid Regions Science Data Center at Lanzhou (http://westdc.westgis.ac.cn). The overall classification accuracy of this land cover map was 71%, which was higher than the accuracies of other land cover maps (i.e., global land cover dataset of IGBP, MODIS land cover map, a global land cover map produced by the University of Maryland, and the land cover map produced by the global land cover for the year 2000 (GLC 2000) project. The overall classification accuracy of these four land cover maps is between 54.17–59.28%) [37,38]. The land cover map uses a unified 1000 m resolution and ALBERS projection. We reprojected it into UTM projection with the WGS84 datum and resampled it into 500 m resolution so that it has consistent projection and resolution with MODIS SCF data.
The NASA Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) data (SRTM3, 90 m) was provided by the Consortium for Spatial Information (CGIAR-CSI) (http://srtm.csi.cgiar.org). The eighteen DEM tiles are mosaicked, reprojected, and resampled to achieve consistency with the MODIS SCF images. The elevation, aspect, and slope maps were extracted from processed DEM data.

3. Methodology

3.1. Investigation on Cloud Gaps in MODIS SCF Products

How much cloud gaps are in the MODIS SCF products? How long is the cloud cover duration? To answer these two questions, we first reclassified the MODIS SCF products into two categories: fractional snow (0 to 100) and cloud (250). After reclassification, annual cloud coverage is calculated to analyze the influence of cloud gaps on MODIS SCF products. Furthermore, to explore the relationship between elevation and cloud cover frequency, we also take elevation into account in our calculations. For example, the daily cloud coverage for elevations above the specified threshold Z T of day d can be calculated using the following formula.
C d ( Z > Z T ) = N d ( Z > Z T ) N T o l ( Z > Z T )
where N d ( Z > Z T ) and N T o l ( Z > Z T ) are the numbers of cloud pixels and total number of pixels in the elevation above Z T for date d in the study area. If Z T is taken as the minimum elevation ( Z m i n ) of study area, the study area is considered in its entirety.

3.2. Gap-Filling Procedure

3.2.1. Theoretical Basis

Existing studies have found that there are always numerous similar pixels in a remote sensing image. Similar pixels do not exist only in an adjacent location but also appear in non-local neighboring regions [35,39]. Similar pixels have the following characteristics. (1) They have approximately equal spectral characteristics and, (2) despite images obtained at different times having different spectral characteristics, similar pixels have approximately the same reflectance changes over time and relative coincident positions in multi-temporal images [40]. For example, Figure 2 shows the snow distribution of MOD10A1 in a small region (the size is 31 × 31 pixels, with the upper left corner is located at (45.05°N, 86.00°E)) of study area on (a) 30 January 2003 and (b) 18 February 2003. The pixels displayed in the red box in Figure 2a are the similar pixels of the target pixel (black box), which are also similar pixels of target pixels in Figure 2b. The land cover types of these pixels are cropland, and these pixels have close reflectance within an image and consistent changes with time. Taking this basic concept, the ground information of a cloud pixel within a MODIS SCF image can be estimated by using the appropriate similar pixels in the remaining known region with the help of multitemporal reference images.
However, some other studies point out that a pixel would also show considerable reflectance change where land cover type changes, therefore causing the phenomenon of temporal difference measurement [41,42]. For example, regarding the pixel we picked out with a green square in Figure 2a, its land cover type is bare land. It is also the similar pixel of the target pixel before snow melts in Figure 2a, but it is no longer the similar pixel of target pixels in Figure 2b after partial snow melting. This may be caused by different snow melting rates under changed underlying surface conditions. To tackle this problem, the pixel which had the same land cover type with the target pixel was regarded as a similar pixel [41,42,43].

3.2.2. Gap-Filling Method

An innovative gap-filling method via NSTF is proposed to reconstruct the ground information of the cloud gaps in the daily MODIS SCF products. As there are significant gaps in large cloud regions, it is unreliable and unrealistic to rely on the remaining known regions of the image to match similar pixels. Therefore, a combination of Terra and Aqua, and temporal filtering are executed first to fill some gaps before NSTF. However, the NSTF method can be used independently. The detailed schematic of the MODIS SCF products gap-filling procedure is presented in Figure 3.
1. Daily MOD10A1 and MYD10A1 Combination
Daily MOD10A1 and MYD10A1 products were combined first by using the composite rules where priority was assigned to MOD10A1, as the validation studies demonstrated that it had more accurate retrievals than MYD10A1 [13]. Hereafter we refer to MOD10A1, MYD10A1, and the combination products to MOD, MYD, and MOYD, respectively. The specific combination strategy is the following. When a pixel is cloud-free in both MOD and MYD, the value in MOD is taken. When a pixel is cloud-free in only one of the products, the cloud free observation is taken. In addition, because the MYD product was not available until July 2002, we only used the MOD product prior to that date.
2. Adjacent Temporal Filtering
Adjacent temporal filtering (ATF) rests on the fundamental assumption that snow will persist on the surface for a certain period of time. Previous studies of one-day to eight-day temporal filters have confirmed that the overall accuracy of multi-days combined snow cover products improves with the increase of composition days [26,27,32,44,45]. However, Gao et al., (2010) found that the overall accuracy of multi-day combined snow cover products decreases with increasing composition days [19]. To determine the optimum composition days, we calculated a cloud persistence index (CPI) for each pixel via finding the average of the continuously cloud covered days [25]. The result indicates that three days is the optimum number of composition days, as the mean CPI for the study area from 2001 to 2016 is 2.72 days. In addition, to guarantee the real-time application ability of the proposed method, we only use images from the previous date. Thus, we determined the three-days backward temporal filtering method used in this study. For a cloudy pixel, we searched images from the previous three days and, if a cloud-free observation was obtained from two or more days, the last valid SCF at this location was assigned to this cloud pixel.
3. NSTF
To match the similar pixels in the remaining non-cloud covered regions to fill the gap pixel, the following three criteria need to be considered. (1) The similar pixels have similar SCF values and have an approximately similar SCF variation in the time domain; (2) Considering the continuity of snow distribution in the spatial domain, the similar pixels have an approximately similar SCF proximity with their respective surrounding neighbors; (3) The similar pixels have approximately similar topographic features. According to the above three criteria, three objective functions are defined. The similar pixels are the solution of this three-objective optimization problem which requires the simultaneous minimization of three objective functions under equality constraints. However, the three objective functions are almost impossible to minimize simultaneously because there is often conflict between multiple objectives. To provide a means to assess trade-off between three conflicting objectives, an automatic machine learning technique of the fast elitist and non-dominated sorting multi-objective genetic algorithm (NSGA-II) is adopted [46,47].
First, for a cloud pixel P on target image, it is assumed that P is the pixel at the same position as P on an auxiliary reference image. The historical image closest to the target image with B ( P ) is a 3 × 3 non-cloud block centered at P (Figure 4). This was selected as the most appropriate auxiliary reference image. We then extracted the M non-cloud common pixels of the target image and reference image within a 51 × 51 pixels search window centered at location P , as is too time-consuming to search for similar pixels in the entire study area [41]. To guarantee the precision of the NSTF, we let M account for at least 30% of the total number of pixels in searching window. If this condition was not satisfied, the search window was expanded to 101 × 101 pixels then 151 × 151 pixels, etc., until the entire study area is encompassed.
Second, the three-objective MODIS SCF products gap-filling problem can be expressed as:
M i n i m i z e : [ f 1 ( P S ) , f 2 ( P S ) , f 3 ( P S ) ]
S u b j e c t t o : L C ( P S ) = L C ( P )
where P s is a candidate similar pixel of P . L C is the land cover type.
A: Minimize Variation in Time Domain
The objective considers that similar pixels having similar change trends over time, and relative coincident positions in multi-temporal images can be mathematically defined as
f 1 ( P S ) = S T ( P S ) S T 2
where the subscript T represents the target image, S T ( P S ) is the SCF value of P S , S T is the average SCF value of pixels in a corresponding location P S in the target image, and P S is a candidate similar pixel of P which was calculated in advance (Figure 4). The calculation procedure of P S is
S R ( B ( P ) ) S R ( B ( P S ) ) 2 d × 2 S R ( P )
S u b j e c t   t o :   L C ( P S ) = L C ( P )
where B ( P S ) is a 3 × 3 non-cloud block centered at location P S ; subscript R represents the reference image; S R ( P ) is the SCF value of P ; and S R ( B ( P ) ) and S R ( B ( P S ) ) represent the average SCF value of blocks B ( P ) and B ( P S ) , respectively. d is a free parameter related to the sensor [41]; it is set as 0.01 for the MODIS case. Equation (4) is the similarity measurement to ensure that the difference between the two blocks is extremely slight.
B: Minimize Discontinuity in Spatial Domain
The objective function f2 considers the spatial continuity of snow cover distribution. The continuity between a cloud pixel and surrounding spatial neighbors can be presented as
f 2 ( P S ) = S T ( P S ) S T ( P ) ¯ 2 + S R ( P S ) S R ( P ) ¯ 2
where P is the spatial non-cloud pixel of first order adjacency to P (Figure 4), and S T ( P ) ¯ and S R ( P ) ¯ represent the average SCF values of P in the target and reference image, respectively.
C: Minimize the Topographic Differences
Similar to [21], we also assume that there is equivalent exposure and the same amount of snow cover under the same terrain conditions. Thus, the objective function f3 is used to enforce similar pixels having similar topographic features.
f 3 ( P S ) = E ( P S ) E ( P S ) 2 + S ( P S ) S ( P S ) 2 + A ( P S ) A ( P S ) 2
where E , S , and A represent normalized elevation, slope and aspect, respectively. Among them, E and S adopt the Min-Max normalization method, and the normalization method of A was: A = | A 0 180 | / 180 ( A 0 was the original aspect value).
Third, the NSGA-II algorithm was executed to search for N (the value is 10 in the present study) similar pixels among the M non-cloud common pixels. We encoded the chromosome ( C M , a solution to this three-objective minimal optimization problem) as the combination of the location of the N similar pixels.
C M = [ x 1 y 1 x 2 y 2 x 3 y 3 x N y N ]
where ( x i , y i ) is the location of i-th similar pixels, i = 1, 2, …, N.
The detailed execution process of NSGA-II can be referred to in [46].
Finally, the arithmetic mean SCF value of N similar pixels (that is, the average SCF value of 10 pixels in corresponding locations of ( x i , y i ) ) is arranged for cloud gap pixel P .

3.3. Validation Methodology

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

3.3.1. Validation Based on “Cloud Assumption”

Considering that percentages of cloud obstruction and the duration of the clouds will affect the accuracy assessment results directly, two strategies are adopted to carry out validation [25].
1. One-Day Cloud Mask
The MODIS SCF image with the lowest cloud coverage per month from 2001 to 2016 is considered the “cloud-free” image. Then, to ensure that the artificial fictitious cloud gaps would occur as close to real life as possible, we calculated the cloud coverage level of 25, 50, and 75% probability of occurrence (PO) via the empirical cumulative distribution function (ECDF) to describe the cloud pattern of low, moderate, and high occurrence in each month. The closest images in time to 25% PO, 50% PO, and 75% PO (abbreviated to P25, P50, and P75 here after) in each month were identified, and the corresponding cloud gaps were extracted to mask the selected “cloud-free” image of this month. Consequently, each selected “cloud-free” SCF image was related to three mask images, and it co-produced a total of 576 validation images during the study period of 16 years.
2. Multi-Day Cloud Mask
To imitate the real natural conditions of consecutive overcast days, the spatio-temporal additional cloud masked images were generated. Because the mean CPI of the study area during the period from 2001 to 2016 is 2.72 days, we chose two groups of three consecutive days characterized by minimum cloud cover every snow season (i.e., November to March next year as a snow season) to be the “cloud-free” images. Additionally, the cloud gaps in the corresponding images of the previous snow cover season records were extracted to overlap onto the selected “cloud-free” images. Moreover, we chose additional groups of two, four, five, and seven consecutive days every snow season to launch more extensive validation. This co-produced a total of 90 validation groups.
Four performance evaluation criteria (i.e., root mean square error ( R M S E ), coefficient of determination ( R 2 ), overestimated error ( O E ), and underestimated error ( U E )) are defined as follows.
R M S E = 1 N i = 1 N ( x i s x i T ) 2
R 2 = ( 1 N 1 i = 1 N ( x i s x i s ¯ σ s ) ( x i T x i T ¯ σ T ) ) 2
O E = ( i = 1 N α i / N ) × 100 % α i = { 1 , x i s x i T 5 % 0 , x i s x i T < 5 %
U E = ( i = 1 N β i / N ) × 100 % β i = { 1 , x i s x i T - 5 % 0 , x i s x i T > - 5 %
where x i s is the estimated SCF value of the i-th fictitious cloud pixel, x i T is the corresponding ground truth MODIS SCF value, and N is the total number of fictitious cloud masked pixels.
In addition to the four above traditional evaluation criteria, a new bias-insensitive metric called spatial efficiency (SPAEF) was used for snow cover spatial pattern validation, the SPAEF is defined as follows [48]:
S P A E F = 1 ( A 1 ) 2 + ( B 1 ) 2 + ( C 1 ) 2
where A is the correlation coefficient between the gap-filled and ground truth MODIS SCF image, B is the fraction of coefficient of variation, and C is the percentage of histogram intersection. The optimal value of S P A E F is 1, which means that the gap-filled fictitious MODIS SCF image is exactly the same as the ground truth MODIS SCF image.

3.3.2. Validation Based on In Situ SD Observations

The original and gap-filled MODIS SCF images were compared with the in situ SD observations corresponding to the 50 meteorological stations. We applied the validation strategy based on the confusion matrix (Table 1), which is recommended by extensive studies such as [49].
a , b , c , and d are the number of pixels correctly hit, commission alarmed, omission alarmed, and correctly rejected, respectively. e and f are the cloud pixels in the MODIS SCF image when the meteorological station reports snow and no snow, respectively. The ε 1 and ε 2 are the defined threshold values of observed SD and MODIS SCF to determine whether a station or a MODIS pixel is covered by snow, respectively.
The four validation indices are defined as follows.
OA = ( ( a + d ) / ( a + b + c + d + e + f ) ) × 100 %
MU = ( c / ( a + b + c + d ) ) × 100 %
MO = ( b / ( a + b + c + d ) ) × 100 %
Bias = ( ( a + b ) / ( a + c + e ) ) × 100 %
Overall accuracy ( OA ) is the probability a MODIS pixel that correctly classified; MU and MO represent the underestimated and overestimated snow event, respectively; and Bias refers to the probability of snow event being correctly identified by MODIS. The perfect estimation of MODIS is OA =   1 , Bias = 1 (we also set Bias to 1 when a = b = c = e = 0 ), and MU = Μ O = 0 . The imperfect estimation of MODIS is OA =   0 , MU = MO = 1 , Bias = 0 or (we set Bias to 0 when a = c = e = 0 and b 0 , which means no snow is observed whereas MODIS indicates snow cover).

4. Experimental Results

4.1. Cloud Gaps in MODIS SCF Products Over Northern Xinjiang

After investigating the cloud gaps in MODIS SCF products of the study area, we found that the extent of cloud obstruction is even worse than we had desired. Annual average cloud coverage of the study area is 43.53–49.85% for Terra and 42.43–52.31% for Aqua from 2001 to 2016, with a mean cloud coverage of 46.48 and 47.78%, respectively (Table 2). The results are effectively consistent with previous studies of [24]. The cloud cover duration maps demonstrate that the mean number of cloud observations is 94–275 days for MODIS Terra SCF images and 96–316 days for MODIS Aqua SCF images, with mean cloud days of 169 and 175, respectively (Figure 5).
We took MODIS Terra SCF products as an example (though there were similar patterns for MODIS Aqua). The calculated annual cloud coverage in different elevation zones are shown in Figure 6. From Figure 6, the annual average cloud coverage for each elevation zone demonstrates that the cloud obstruction probability increases with elevation. The annual average cloud percentages are 42–48.94, 43.05–49.32, 45.07–52.12, 50.26–57.31, and 43.53–49.91 for elevations of <1000, 1000–2000, 2000–3000, >3000 m, and the entire study area, with average values of 45.31, 46.30, 48.61, 53.09, and 46.52%, respectively. The systematic trends of cloud obstruction for the entire study area are the most consistent with those of the elevation zone of 2000–3000 m. The influence of topography on cloud cover is also confirmed in Figure 5, in which the number of cloud cover days have shown a positive pattern in relation to elevation. This result is partly attributable to the water vapor uplift caused by orography, making its condensation into clouds more likely. Meanwhile, atmospheric circulation may be contributing to clearing the clouds at low-altitude regions [37].

4.2. The Effectiveness of Gap-Filling Methodology

The proposed gap-filling methodology generated almost continuous spatio-temporal daily MODIS SCF products, leaving only 0.52% of cloud gaps on long-term average. The maximum gaps remaining is 5.92%. The mean cloud coverage remaining in different elevation zones after the execution of each gap-filling step is presented in Figure 7. We found an extremely consistent change trend between the single year of 2004 and the average after 16 years. Similar conditions were found in other years, which will not be explored here. MOYD, ATF, and NSTF reduced cloud cover by an average of 8.19% (compared to the original Terra data) or 9.72% (compared to the original Aqua data), 21.4%, and 16.43%, respectively.
Figure 8 shows an example of the procedure applied to gap filling over four days for 1 January 2001, 1 January 2005, 1 January 2009 and 1 January 2013. For these four days, Terra and Aqua maps had a cloud coverage of 85.6, 86.92, 70, and 47.58% and 84.72, 91.87, 76.1, and 61.43%, respectively. MOYD left out a cloud coverage of 84.94, 83.09, 62.73, and 41.85%. The intervention of ATF affects 25.39, 70.87, 45.25, and 29.78% of the cloud coverage on those days, respectively. NSTF removes the remaining 59.55, 12.22, 17.48, and 12.07% of the clouds on those days, respectively, and closes the gap-filling procedure. By consecutively performing three gap-filling steps, even if more than 90% of the ground information is obscured by clouds, extensive cloud gaps can be reconstructed and more reasonable ground information is obtained. This agrees with the general knowledge of snow spatial distribution characteristics in Northern XinJiang.

4.3. Accuracy Validation Results

4.3.1. Validation Based on “Cloud Assumption”

1. Accuracy: An Example of One-Day Cloud Mask
We take the special case of January 2004 as an example. The clearest MODIS Terra SCF image appeared on 21 January with 18.18% cloud gaps. This was regarded as the “cloud-free” MODIS SCF image (Table 3). The ECDF of cloud coverage in each month for the MODIS Terra SCF product in 2004 is shown in Figure 9. The closest images to P25, P50, and P75 in January were the images from 1 January 2004, 20 January 2004, and 8 January 2004, with cloud coverages of 47.55, 55.08 and 69%, respectively. After these three extracted images were assigned to the clearest SCF image on 21 January, the percentage of cloud cover in the three artificially cloud-filled images were 53.59, 59.78 and 71.87%, respectively. Note that we also performed this cloud-filled procedure for the corresponding MODIS Aqua SCF image.
After the gap-filling procedure was finished, the ground information was perfectly reconstructed visually. Despite the clouds obscuring more than half of the study area, the SCF values of the reconstructed image were extremely close to the true value of the “cloud-free” MODIS SCF image in all three cloud patterns (Figure 10). MOYD reduced 14.75, 19.85, and 17.31% of the cloud cover and ATF removed 33.7, 45.02, and 49.47% of the cloud cover for P25, P50, and P75 cloud-masked conditions, respectively. NSTF eliminates the remaining clouds and completely cloudless images were produced. As seen from Table 4, the proposed gap-filling approach can achieve the preferred accuracy. The R M S E was 0.07–0.13, O E was 0.59%–2.07%, U E was below 1%, and R 2 was 0.86–0.96. The spatial pattern validation metric of S P A E F for the MODY, ATF, and NSTF steps, and final images exceeded 0.8, 0.65, 0.7, and 0.75, respectively. This illustrates that the gap-filled fictitious MODIS SCF image is very consistent with assumed ground truth MODIS SCF image. As expected, there was better gap-filling performance under lower cloud coverage level conditions. In addition, the performance of MOYD was superior to ATF, while ATF was superior to NSTF in most cases. This was because the implementation of the second step was based on the results of the previous step. As such, the imprecision of the previous step was passed through the input information of the second step and error was accumulated.
2. Accuracy: An Example of Multi-Day Cloud Mask
In view of the fact that the reconstruction of SCF information under consecutive cloud cover in the snow season has more significant importance for practical applications, we take the special case of the 2015 to 2016 snow season as an example. Three consecutive days characterized by lowest average cloud obstruction for MODIS Terra SCF images on February 26, February 27, and February 28 of 2016 are selected as the “cloud-free” true SCF images. The cloud coverages during these three days were 17.50, 13.12, and 64.68%, respectively (Table 5). The corresponding images on the same dates in 2015 had cloud coverages of 66.03, 44.06, and 45.50% were used to mask the selected “cloud-free” images. The three artificially cloud-filled images were produced with cloud coverages of 70.43, 50.44, and 79.59%. The same was also done for the corresponding Aqua images.
From Table 6, we can see that the performance of this consecutive cloud cover case is slightly worse, compared to the one-day cloud mask case. MOYD reduced the cloud cover by 8.8, 16.92, and 16.73% on the three days, respectively. ATF reduced the cloud cover by 34.71, 11.58, and 35.59% on the three days, respectively. Although the NSTF step eliminated a lot of cloud cover (26.02, 20.94, and 26.42% on the three days, respectively), it does not completely eliminate the cloud cover. Tiny amounts (i.e., 0.9, 1, and 0.95% on the three days, respectively) of cloud gaps were still retained. Higher accuracy was also achieved with lower overall R M S E values of 0.09–0.11, O E values of 1.15–1.58%, U E values of 0.15–2.41%, higher R 2 values of 0.93–0.96, and higher S P A E F values of 0.8–0.88. Despite the significant differences in terms of both cloud percentage and cloud patterns, there was little performance difference among these three days.
3. Validation Results Based on “Cloud Assumption”
The generalization validation results of 576 one-day and 90 groups of multi-day (2–7 days) cloud masked images are shown in Figure 11. The comparison of the result of each gap-filing step with the selected “cloud-free” true MODIS SCF image confirms the efficacy of the proposed gap-filling method. For the one-day case, most cloud-free images on each month can be found with an average cloud coverage less than 1%. The average cloud coverage of the artificially cloud-filled images was 30, 47, and 66%, respectively, when P25, P50 and P75 additional cloud cover was introduced. The average cloud coverage of the selected 2–7 consecutive “cloud-free” days was approximately 11–27% and 50–64% after artificially spatiotemporal additional cloud cover was incorporated, as shown in Figure 11a. The gap-filling efficiency is further confirmed by Figure 11b. More than 99% of gaps were reconstructed. The average amount of cloud gaps reduced by MODY, ATF, and NSTF are 9.21, 25.19, and 17.28%, which are basically similar with the results averaged from 16 years given in Section 4.2. Based on Figure 11c–f and Figure 12, although there are significant differences in both cloud percentage and cloud patterns among different validation strategies, all of them are characterized by lower values of R M S E , O E , and U E , as well as higher values of R 2 and S P A E F . The average R M S E was about 0.1, R 2 was greater than 0.8, O E was 1.13%, U E was 1.4%, and S P A E F was 0.78 for the final images. The performance of MOYD and ATF was very prominent, with the S P A E F below 0.1, the O E and U E below 1%, and S P A E F exceeding 0.8 in all validation strategies. However, the performance of NSTF was slightly worse, particularly for the validation strategy of four to seven consecutive cloud covered days. The R M S E reached up to 0.15, some disagreement errors exceeded 2%, and the average S P A E F was approximately 0.65. This is because, after three days of backward temporal filtering in ATF, it will inevitably leave out a larger percentage of gaps. There will also be less available non-cloud pixels for NSTF. Together with the accumulative effect of error, these factors lead to a slightly inferior performance of NSTF.

4.3.2. Validation Based on In-Situ SD Observations

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

5. Discussion

5.1. The Sensitivity of Different SCF and SD Threshold Combinations

When comparing the MODIS SCF value with in situ SD observations, the most suitable SD threshold value ( ε 1 ) and SCF threshold value ( ε 2 ) must be chosen. ε 1 = 1 [26,51,52], ε 1 = 3 [21,28], ε 1 = 4 [22,24,30,31], ε 2 = 15 % [53,54], and ε 2 = 50 % [12,55] have been used in the literature. We use these thresholds to analyze the sensitivity of different thresholds to validate results, as shown in Figure 14.
The results indicated that a higher OA (above 90%) and lower MU and MO (below 2%) are achieved for all the different threshold combinations in the non-snow season, but Bias is closest to 1 when the threshold value combination of ε 1 = 3 and ε 2 = 50 % is used. This illustrates that the inconsistency between the final gap-filled MODIS SCF product and in situ SD observations is minimal. In snow seasons, OA is over 80% and Bias is also extremely close to 1 for all different threshold combinations, which further confirms that there is favorable consistency between the final gap-filled MODIS SCF product and in situ SD observations. The MODIS SCF threshold value of ε 2 = 50 % produced a relatively a higher MU and lower MO . When the ε 2 is fixed, a large value of ε 1 produces relatively lower MU and higher MO values. For example, MU is 12% and MO is 0 for the threshold combination of ε 1 = 1 and ε 2 = 50 % , and MU is 8.2% and MO is 9.08% for the threshold combination of ε 1 = 4 and ε 2 = 50 % in March. The above results demonstrate that the most suitable ε 1 should be between 1 and 4cm and ε 2 should be between 15 and 50%. Therefore, we drew the inference that the threshold combination of ε 1 = 3 and ε 2 = 50 % is the best choice for the MODIS accuracy assessment in the study area.

5.2. Comparison with Another Gap-Filling Method

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

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

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

6. Conclusions

Due to atmospheric disturbances, much ground information usually cannot be acquired. This is what mainly limit MODIS SCF products from being widely applied. Addressing this challenge, an innovative gap-filling method for MODIS SCF products based on the concept of NSTF is proposed. Ground information of a gap pixel is estimated through selected similar pixels in the remaining known region via an automatic machine learning technique of NAGA-II. As there are significant gaps in large cloud regions, it is unreliable and unrealistic to rely on the remaining known regions of the image to match similar pixels. The widely confirmed method of MOYD and ATF are first executed sequentially to reconstruct some gap pixels. However, the NSTF method can be used independently.
In this study, we first answered the two questions on how much cloud gaps are in the MODIS SCF products and how long is the cloud cover duration. The influence of topography on cloud cover is also evident. Almost cloud-free daily MODIS SCF images can be produced after consecutive implementation of three cloud removal steps. The long-term mean cloud coverage is 0.52% on average. The generalization validation results of 576 one-day cloud masked images and 90 groups of multi-day (2–7 days) spatiotemporal additional cloud masked images prove the reliability and feasibility of the proposed methodology. The average R M S E is approximately 0.1, R 2 exceeds 0.8, O E error is 1.13%, U E error is 1.4%, and S P A E F is 0.78, even if more than 80% of the surface covering information is missing. The validation based on 50 in situ SD observations also demonstrates the superiority of the method in terms of accuracy and consistency. The average MU and MO error have increased about 1.16 and 0.53% for the study area in 16 years. In addition, the validation test in terms of land cover types also proves the efficiency and accuracy of the proposed gap-filling method. In general, existing MODIS cloud gap-filling methods have satisfactory accuracy in the stable snow cover phase, but they cannot achieve a satisfying accuracy in the transition phase of snow cover accumulation and ablation [10,14]. Different from the existing methods, the proposed method has a no significant performance difference in seasonal or different elevation zones. NSTF can consider both the spatio-temporal correlation of multitemporal MODIS SCF products and auxiliary information, such as topography and land cover, in an integrated manner, instead of analyzing each individually. Thus, it can more fully exploit the space-time correlation within a pixel. Nevertheless, there are still some limitations that should be improved. For example, we simply used a fixed value of 10 as the number of similar pixels. This is obviously defective. We will use adaptive variable values according to the cloud coverage level in subsequent studies.

Author Contributions

C.H. designed the whole study. J.H. developed the algorithm, conducted the data analysis, and crafted the manuscript. Y.Z., J.G. (Jifu Guo) and J.G. (Juan Gu) significantly edited and revised the manuscript.

Funding

This research was funded by the Strategic Priority Research Program of Chinese Academy of Sciences (grant No. XDA19040504) and the National Natural Science Foundation of China (Grant No. 41501412 and 41671375).

Acknowledgments

The MODIS data are available from the NASA Distributed Active Archive Center at NSIDC. The snow depth data are available from Chinese Meteorological Data Sharing Service System (CMDS), The land cover data are available from the Cold and Arid Regions Science Data Center and the DEM data are available from the Consortium for Spatial Information (CGIAR-CSI), who we would like to thank.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript and in the decision to publish the results.

References

  1. Robinson, D.A.; Dewey, K.F.; Heim, R.R. Global Snow Cover Monitoring—An Update. Bull. Am. Meteorol. Soc. 1993, 74, 1689–1696. [Google Scholar] [CrossRef]
  2. Brown, R.D. Northern hemisphere snow cover variability and change, 1915-97. J. Clim. 2000, 13, 2339–2355. [Google Scholar] [CrossRef]
  3. Barry, R.G. The role of snow and ice in the global climate system: A review. Polar Geogr. 2002, 26, 235–246. [Google Scholar] [CrossRef]
  4. Musselman, K.N.; Clark, M.P.; Liu, C.H.; Ikeda, K.; Rasmussen, R. Slower snowmelt in a warmer world. Nat. Clim. Chang. 2017, 7, 214–219. [Google Scholar] [CrossRef]
  5. Barnett, T.P.; Adam, J.C.; Lettenmaier, D.P. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 2005, 438, 303–309. [Google Scholar] [CrossRef] [PubMed]
  6. Brown, L.; Thorne, R.; Woo, M.K. Using satellite imagery to validate snow distribution simulated by a hydrological model in large northern basins. Hydrol. Process. 2008, 22, 2777–2787. [Google Scholar] [CrossRef]
  7. Estilow, T.W.; Young, A.H.; Robinson, D.A. A long-term Northern Hemisphere snow cover extent data record for climate studies and monitoring. Earth Syst. Sci. Data 2015, 7, 137–142. [Google Scholar] [CrossRef] [Green Version]
  8. Hall, D.K.; Riggs, G.A.; Salomonson, V.V.; DiGirolamo, N.E.; Bayr, K.J. MODIS snow-cover products. Remote Sens. Environ. 2002, 83, 181–194. [Google Scholar] [CrossRef] [Green Version]
  9. Frei, A.; Tedesco, M.; Lee, S.; Foster, J.; Hall, D.K.; Kelly, R.; Robinson, D.A. A review of global satellite-derived snow products. Adv. Sp. Res. 2012, 50, 1007–1029. [Google Scholar] [CrossRef] [Green Version]
  10. Klein, A.G.; 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. 2003, 86, 162–176. [Google Scholar] [CrossRef]
  11. Huang, X.D.; Liang, T.G.; Zhang, X.T.; Guo, Z.G. Validation of MODIS snow cover products using Landsat and ground measurements during the 2001–2005 snow seasons over northern Xinjiang, China. Int. J. Remote Sens. 2011, 32, 133–152. [Google Scholar] [CrossRef]
  12. Tang, Z.; Wang, J.; Li, H.; Yan, L. Spatiotemporal changes of snow cover over the Tibetan plateau based on cloud-removed moderate resolution imaging spectroradiometer fractional snow cover product from 2001 to 2011. J. Appl. Remote Sens. 2013, 7. [Google Scholar] [CrossRef]
  13. Yang, J.T.; Jiang, L.M.; Menard, C.B.; Luojus, K.; Lemmetyinen, J.; Pulliainen, J. Evaluation of snow products over the Tibetan Plateau. Hydrol. Process. 2015, 29, 3247–3260. [Google Scholar] [CrossRef]
  14. Simic, A.; Fernandes, R.; Brown, R.; Romanov, P.; Park, W. Validation of VEGETATION, MODIS, and GOES plus SSM/I snow-cover products over Canada based on surface snow depth observations. Hydrol. Process. 2004, 18, 1089–1104. [Google Scholar] [CrossRef]
  15. Parajka, J.; Bloschl, G. Validation of MODIS snow cover images over Austria. Hydrol. Earth Syst. Sci. 2006, 10, 679–689. [Google Scholar] [CrossRef] [Green Version]
  16. Gafurov, A.; Kriegel, D.; Vorogushyn, S.; Merz, B. Evaluation of remotely sensed snow cover product in Central Asia. Hydrol. Res. 2013, 44, 506–522. [Google Scholar] [CrossRef]
  17. Lin, C.-H.; Tsai, P.-H.; Lai, K.-H.; Chen, J.-Y. Cloud Removal from Multitemporal Satellite Images Using Information Cloning. IEEE Trans. Geosci. Remote Sens. 2013, 51, 232–241. [Google Scholar] [CrossRef]
  18. Tong, J.; Dery, S.J.; Jackson, P.L. Interrelationships between MODIS/Terra remotely sensed snow cover and the hydrometeorology of the Quesnel River Basin, British Columbia, Canada. Hydrol. Earth Syst. Sci. 2009, 13, 1439–1452. [Google Scholar] [CrossRef] [Green Version]
  19. Gao, Y.; Xie, H.J.; Yao, T.D.; 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. 2010, 114, 1662–1675. [Google Scholar] [CrossRef]
  20. Dietz, A.J.; Kuenzer, C.; Conrad, C. Snow-cover variability in central Asia between 2000 and 2011 derived from improved MODIS daily snow-cover products. Int. J. Remote Sens. 2013, 34, 3879–3902. [Google Scholar] [CrossRef]
  21. López-Burgos, V.; Gupta, H.V.; Clark, M. Reducing cloud obscuration of MODIS snow cover area products by combining spatio-temporal techniques with a probability of snow approach. Hydrol. Earth Syst. Sci. 2013, 17, 1809–1823. [Google Scholar] [CrossRef] [Green Version]
  22. Parajka, J.; Pepe, M.; Rampini, A.; Rossi, S.; Blöschl, G. A regional snow-line method for estimating snow cover from MODIS during cloud cover. J. Hydrol. 2010, 381, 203–212. [Google Scholar] [CrossRef]
  23. Liang, T.G.; Zhang, X.T.; Xie, H.J.; Wu, C.X.; Feng, Q.S.; Huang, X.D.; Chen, Q.G. Toward improved daily snow cover mapping with advanced combination of MODIS and AMSR-E measurements. Remote Sens. Environ. 2008, 112, 3750–3761. [Google Scholar] [CrossRef]
  24. Gao, Y.; Xie, H.; Lu, N.; Yao, T.; Liang, T. Toward advanced daily cloud-free snow cover and snow water equivalent products from Terra–Aqua MODIS and Aqua AMSR-E measurements. J. Hydrol. 2010, 385, 23–35. [Google Scholar] [CrossRef]
  25. Wang, X.W.; Xie, H.J.; Liang, T.G.; Huang, X.D. Comparison and validation of MODIS standard and new combination of Terra and Aqua snow cover products in northern Xinjiang, China. Hydrol. Process. 2009, 23, 419–429. [Google Scholar] [CrossRef]
  26. Dariane, A.B.; Khoramian, A.; Santi, E. Investigating spatiotemporal snow cover variability via cloud-free MODIS snow cover product in Central Alborz Region. Remote Sens. Environ. 2017, 202, 152–165. [Google Scholar] [CrossRef]
  27. Parajka, J.; Blöschl, G. Spatio-temporal combination of MODIS images—Potential for snow cover mapping. Water Resour. Res. 2008, 44. [Google Scholar] [CrossRef]
  28. Gafurov, A.; Bardossy, A. Cloud removal methodology from MODIS snow cover product. Hydrol. Earth Syst. Sci. 2009, 13, 1361–1373. [Google Scholar] [CrossRef] [Green Version]
  29. Dong, C.; Menzel, L. Improving the accuracy of MODIS 8-day snow products with in situ temperature and precipitation data. J. Hydrol. 2016, 534, 466–477. [Google Scholar] [CrossRef]
  30. Krajci, P.; Holko, L.; Perdigao, R.A.P.; Parajka, J. Estimation of regional snowline elevation (RSLE) from MODIS images for seasonally snow covered mountain basins. J. Hydrol. 2014, 519, 1769–1778. [Google Scholar] [CrossRef]
  31. Da Ronco, P.; De Michele, C. Cloud obstruction and snow cover in Alpine areas from MODIS products. Hydrol. Earth Syst. Sci. 2014, 18, 4579–4600. [Google Scholar] [CrossRef]
  32. Huang, X.; Deng, J.; Wang, W.; Feng, Q.; Liang, T. Impact of climate and elevation on snow cover using integrated remote sensing snow products in Tibetan Plateau. Remote Sens. Environ. 2017, 190, 274–288. [Google Scholar] [CrossRef]
  33. Hall, D.K.; Riggs, G.A.; Foster, J.L.; Kumar, S.V. Development and evaluation of a cloud-gap-filled MODIS daily snow-cover product. Remote Sens. Environ. 2010, 114, 496–503. [Google Scholar] [CrossRef]
  34. Dozier, J.; Painter, T.H.; Rittger, K.; Frew, J.E. Time–space continuity of daily maps of fractional snow cover and albedo from MODIS. Adv. Water Res. 2008, 31, 1515–1526. [Google Scholar] [CrossRef]
  35. Cheng, Q.; Shen, H.; Zhang, L.; Yuan, Q.; Zeng, C. Cloud removal for remotely sensed images by similar pixel replacement guided with a spatio-temporal MRF model. ISPRS J. Photogramm. Remote Sens. 2014, 92, 54–68. [Google Scholar] [CrossRef]
  36. Ding, Y.; Li, Y.; Li, L.; Yao, N.; Hu, W.; Yang, D.; Chen, C. Spatiotemporal variations of snow characteristics in Xinjiang, China over 1961–2013. Hydrol. Res. 2018, 49, 1578–1593. [Google Scholar] [CrossRef]
  37. Ran, Y.H.; Li, X.; Lu, L. Evaluation of four remote sensing based land cover products over China. Int. J. Remote Sens. 2010, 31, 391–401. [Google Scholar] [CrossRef]
  38. Ran, Y.H.; Li, X.; Lu, L.; Li, Z.Y. Large-scale land cover mapping with the integration of multi-source information based on the Dempster-Shafer theory. Int. J. Geogr. Inf. Sci. 2012, 26, 169–191. [Google Scholar] [CrossRef]
  39. Chen, J.; Zhu, X.; Vogelmann, J.E.; Gao, F.; Jin, S. A simple and effective method for filling gaps in Landsat ETM+ SLC-off images. Remote Sens. Environ. 2011, 115, 1053–1064. [Google Scholar] [CrossRef]
  40. Zhu, X.; Liu, D.; Chen, J. A new geostatistical approach for filling gaps in Landsat ETM+ SLC-off images. Remote Sens. Environ. 2012, 124, 49–60. [Google Scholar] [CrossRef]
  41. Cheng, Q.; Liu, H.; Shen, H.; Wu, P.; Zhang, L. A Spatial and Temporal Nonlocal Filter-Based Data Fusion Method. IEEE Trans. Geosci. Remote Sens. 2017, 55, 4476–4488. [Google Scholar] [CrossRef] [Green Version]
  42. Shen, H.F.; Wu, P.H.; Liu, Y.L.; Ai, T.H.; Wang, Y.; Liu, X.P. A spatial and temporal reflectance fusion model considering sensor observation differences. Int. J. Remote Sens. 2013, 34, 4367–4383. [Google Scholar] [CrossRef]
  43. Gao, F.; Masek, J.; Schwaller, M.; Hall, F. On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance. IEEE Trans. Geosci. Remote 2006, 44, 2207–2218. [Google Scholar] [CrossRef]
  44. Gao, Y.; Xie, H.J.; Yao, T.D. Developing Snow Cover Parameters Maps from MODIS, AMSR-E, and Blended Snow Products. Photogramm. Eng. Remote Sens. 2011, 77, 351–361. [Google Scholar] [CrossRef]
  45. Paudel, K.P.; Andersen, P. Monitoring snow cover variability in an agropastoral area in the Trans Himalayan region of Nepal using MODIS data with improved cloud removal methodology. Remote Sens. Environ. 2011, 115, 1234–1246. [Google Scholar] [CrossRef]
  46. Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput. 2002, 6, 182–197. [Google Scholar] [CrossRef] [Green Version]
  47. Ramesh, S.; Kannan, S.; Baskar, S. Application of modified NSGA-II algorithm to multi-objective reactive power planning. Appl. Soft Comput. 2012, 12, 741–753. [Google Scholar] [CrossRef]
  48. Demirel, M.C.; Mai, J.; Mendiguren, G.; Koch, J.; Samaniego, L.; Stisen, S. Combining satellite data and appropriate objective functions for improved spatial pattern performance of a distributed hydrologic model. Hydrol. Earth Syst. Sci. 2018, 22, 1299–1315. [Google Scholar] [CrossRef] [Green Version]
  49. Wilks, D.S. Statistical Methods in the Atmospheric Sciences (International Geophysics Series; V. 91); Academic Press: Cambridge, MA, USA, 2006. [Google Scholar]
  50. Riggs, G.A.; Hall, D.; Salomonson, V. MODIS Snow Products User Guide; NASA Goddard Space Flight Center Rep.: Laurel, MD, USA, 2006.
  51. Maurer, E.P.; Rhoads, J.D.; Dubayah, R.O.; Lettenmaier, D.P. Evaluation of the snow-covered area data product from MODIS. Hydrol. Process. 2003, 17, 59–71. [Google Scholar] [CrossRef]
  52. Xu, W.F.; Ma, H.Q.; Wu, D.H.; Yuan, W.P. Assessment of the Daily Cloud-Free MODIS Snow-Cover Product for Monitoring the Snow-Cover Phenology over the Qinghai-Tibetan Plateau. Remote Sens.-Basel 2017, 9. [Google Scholar] [CrossRef]
  53. Painter, T.H.; Rittger, K.; McKenzie, C.; Slaughter, P.; Davis, R.E.; Dozier, J. Retrieval of subpixel snow covered area, grain size, and albedo from MODIS. Remote Sens. Environ. 2009, 113, 868–879. [Google Scholar] [CrossRef] [Green Version]
  54. Rittger, K.; Painter, T.H.; Dozier, J. Assessment of methods for mapping snow cover from MODIS. Adv. Water Res. 2013, 51, 367–380. [Google Scholar] [CrossRef]
  55. Marchane, A.; Jarlan, L.; Hanich, L.; Boudhar, A.; Gascoin, S.; Tavernier, A.; Filali, N.; Le Page, M.; Hagolle, O.; Berjamy, B. Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range. Remote Sens. Environ. 2015, 160, 72–86. [Google Scholar] [CrossRef]
Figure 1. Topographic relief and the location of meteorological observation stations in Northern Xinjiang, Northwest China.
Figure 1. Topographic relief and the location of meteorological observation stations in Northern Xinjiang, Northwest China.
Remotesensing 11 00090 g001
Figure 2. Similar pixels in two different scenes: (a) 30 January 2003 and (b) 18 February 2003.
Figure 2. Similar pixels in two different scenes: (a) 30 January 2003 and (b) 18 February 2003.
Remotesensing 11 00090 g002
Figure 3. Schematic of the MODIS fractional snow cover (SCF) products gap-filling procedure.
Figure 3. Schematic of the MODIS fractional snow cover (SCF) products gap-filling procedure.
Remotesensing 11 00090 g003
Figure 4. Sketch map of the similar pixels. (a) Target image T , in which P is a cloud pixel, P (gray-marked pixel) is the spatial non-cloud pixel of first order adjacency to P , P s is a candidate similar pixel of P ; (b) Reference image R , in which P is non-cloud pixel in the same position relative to P , B ( P ) is a 3 × 3 non-cloud block (red box) centered at P , P s is a candidate similar pixels of P .
Figure 4. Sketch map of the similar pixels. (a) Target image T , in which P is a cloud pixel, P (gray-marked pixel) is the spatial non-cloud pixel of first order adjacency to P , P s is a candidate similar pixel of P ; (b) Reference image R , in which P is non-cloud pixel in the same position relative to P , B ( P ) is a 3 × 3 non-cloud block (red box) centered at P , P s is a candidate similar pixels of P .
Remotesensing 11 00090 g004
Figure 5. Cloud cover duration maps of the standard (a) MOD: MOD10A1 and (b) MYD: MYD10A1.
Figure 5. Cloud cover duration maps of the standard (a) MOD: MOD10A1 and (b) MYD: MYD10A1.
Remotesensing 11 00090 g005
Figure 6. Annual average cloud coverage of MODIS Terra SCF product in different elevation zones.
Figure 6. Annual average cloud coverage of MODIS Terra SCF product in different elevation zones.
Remotesensing 11 00090 g006
Figure 7. Mean cloud coverage remaining in different elevation zones after the execution of each gap-filling step. (a) 2004 and (b) the average of 16 years. MOD: MOD10A1; MYD: MYD10A1; MOYD: daily combinated MOD and MYD product; ATF: adjacent temporal filtering; NSTF: non-local spatio-temporal filter.
Figure 7. Mean cloud coverage remaining in different elevation zones after the execution of each gap-filling step. (a) 2004 and (b) the average of 16 years. MOD: MOD10A1; MYD: MYD10A1; MOYD: daily combinated MOD and MYD product; ATF: adjacent temporal filtering; NSTF: non-local spatio-temporal filter.
Remotesensing 11 00090 g007
Figure 8. Cloud coverage of the original MODIS SCF images after the execution of each gap-filling step. The results from four days (1 January 2001, 1 January 2005, 1 January 2009, and 1 January 2013) are displayed from top to bottom. The columns, from left to right, represent MOD, MYD, and the three consecutive applications of MOYD, ATF and NSTF, respectively.
Figure 8. Cloud coverage of the original MODIS SCF images after the execution of each gap-filling step. The results from four days (1 January 2001, 1 January 2005, 1 January 2009, and 1 January 2013) are displayed from top to bottom. The columns, from left to right, represent MOD, MYD, and the three consecutive applications of MOYD, ATF and NSTF, respectively.
Remotesensing 11 00090 g008
Figure 9. The empirical cumulative distribution function (ECDF) of cloud cover in each month for Terra MOD10A1 in 2004. The red, green, and blue points represent the cloud coverage of the closest images to P25, P50, and P75, respectively.
Figure 9. The empirical cumulative distribution function (ECDF) of cloud cover in each month for Terra MOD10A1 in 2004. The red, green, and blue points represent the cloud coverage of the closest images to P25, P50, and P75, respectively.
Remotesensing 11 00090 g009
Figure 10. Gap-filling procedure for January 2004. (a) Selected “cloud-free” true MODIS SCF image. (b1, c1, d1) Extracted mask images with cloud conditions of P25, P50, and P75. (b2, c2, d2) The corresponding artificially cloud-filled images. The third to fifth rows represent the consecutive applications of MOYD, ATF, and NSTF.
Figure 10. Gap-filling procedure for January 2004. (a) Selected “cloud-free” true MODIS SCF image. (b1, c1, d1) Extracted mask images with cloud conditions of P25, P50, and P75. (b2, c2, d2) The corresponding artificially cloud-filled images. The third to fifth rows represent the consecutive applications of MOYD, ATF, and NSTF.
Remotesensing 11 00090 g010
Figure 11. Performance evaluation results of validation based on “cloud assumption”. (a) Cloud coverage of selected “cloud-free” MODIS Terra SCF images (CT), additional cloud cover introduced by mask images (CA), and artificially cloud-filled images (CM). (bf) Cloud removal efficiency, R M S E , R 2 , O E , and U E of MODY, ATF, NSTF, and final images.
Figure 11. Performance evaluation results of validation based on “cloud assumption”. (a) Cloud coverage of selected “cloud-free” MODIS Terra SCF images (CT), additional cloud cover introduced by mask images (CA), and artificially cloud-filled images (CM). (bf) Cloud removal efficiency, R M S E , R 2 , O E , and U E of MODY, ATF, NSTF, and final images.
Remotesensing 11 00090 g011
Figure 12. Spatial pattern validation results for the spatial efficiency (SPAEF) of MODY, adjacent temporal filtering (ATF), non-local spatio-temporal filtering (NSTF), and final images.
Figure 12. Spatial pattern validation results for the spatial efficiency (SPAEF) of MODY, adjacent temporal filtering (ATF), non-local spatio-temporal filtering (NSTF), and final images.
Remotesensing 11 00090 g012
Figure 13. Validation results based on in situ SD observations in different elevation zones.
Figure 13. Validation results based on in situ SD observations in different elevation zones.
Remotesensing 11 00090 g013
Figure 14. Monthly/annual variation of four validation indices based on different SD threshold values ( ε 1 ) and SCF threshold value ( ε 2 ). Validation strategies 1–6 represent threshold value combinations of ε 1 = 4 and ε 2 = 50 % , ε 1 = 3 and ε 2 = 50 % , ε 1 = 1 and ε 2 = 50 % , ε 1 = 4 and ε 2 = 15 % , ε 1 = 3 and ε 2 = 15 % , and ε 1 = 4 and ε 2 = 15 % , respectively.
Figure 14. Monthly/annual variation of four validation indices based on different SD threshold values ( ε 1 ) and SCF threshold value ( ε 2 ). Validation strategies 1–6 represent threshold value combinations of ε 1 = 4 and ε 2 = 50 % , ε 1 = 3 and ε 2 = 50 % , ε 1 = 1 and ε 2 = 50 % , ε 1 = 4 and ε 2 = 15 % , ε 1 = 3 and ε 2 = 15 % , and ε 1 = 4 and ε 2 = 15 % , respectively.
Remotesensing 11 00090 g014
Figure 15. Monthly/annual variation of four validation indices for the proposed method and the existing cubic spline temporal interpolation approach (spline).
Figure 15. Monthly/annual variation of four validation indices for the proposed method and the existing cubic spline temporal interpolation approach (spline).
Remotesensing 11 00090 g015
Table 1. Confusion matrix comparing MODIS SCF value with snow depth (SD) observation.
Table 1. Confusion matrix comparing MODIS SCF value with snow depth (SD) observation.
Observed SD
Snow   ( ε 1   cm ) No   Snow   ( < ε 1   cm )
MODIS SCFSnow ( ε 2 )ab
No snow ( < ε 2 )cd
Cloudef
Table 2. Mean cloud coverage in the MODIS SCF products.
Table 2. Mean cloud coverage in the MODIS SCF products.
YearMean Cloud Coverage (%)
MOD10A1MYD10A1
200147.20-
200246.9542.43 *
200348.3945.08
200446.5248.37
200544.7045.97
200644.8146.81
200743.6047.31
200843.5346.06
200947.4650.11
201048.6949.88
201146.1548.39
201244.2146.40
201345.8448.49
201446.4848.60
201549.8552.31
201649.3350.54
Average46.4847.78
* represents only the data from 4 July 2002 to 31 December 2002 was used to calculate the annual average cloud coverage.
Table 3. The cloud coverage of selected “cloud-free” MODIS SCF images, three extracted mask images, three artificially cloud-filled images, and after the execution of each gap-filling step for the one-day cloud mask.
Table 3. The cloud coverage of selected “cloud-free” MODIS SCF images, three extracted mask images, three artificially cloud-filled images, and after the execution of each gap-filling step for the one-day cloud mask.
Cloud Coverage (%)
P25P50P75
Original terra18.18
Mask image47.5555.0869.00
Cloud-filled image53.5959.7871.87
MOYD38.8439.9354.56
ATF19.8914.7622.40
NSTF000
Original Terra represents the MODIS Terra SCF image with the lowest cloud cover in January 2004; Mask image represents three extracted images closest to P25, P50, and P75 in January, respectively; Cloud-filled image represents mask images were artificially assigned to the original terra.
Table 4. Quantitative performance evaluation results for the one-day cloud mask.
Table 4. Quantitative performance evaluation results for the one-day cloud mask.
R M S E R 2 O E   ( % ) U E   ( % ) S P A E F
P25P50P75P25P50P75P25P50P75P25P50P75P25P50P75
MOYD0.070.110.130.950.890.840.370.891.460.360.891.480.890.860.82
ATF0.070.130.140.960.860.850.381.992.570.270.670.770.660.660.72
NSTF0.080.110.100.950.900.911.22.361.660.280.450.510.700.760.83
Final images0.070.120.130.960.880.860.591.652.070.30.720.930.910.750.80
Table 5. Cloud coverage of selected “cloud-free” MODIS Terra SCF images, three extracted mask images, three artificially filled images, and after the implementation of each gap-filling step for a consecutive three-days cloud mask. D1–D3 refers to the selected first to third days.
Table 5. Cloud coverage of selected “cloud-free” MODIS Terra SCF images, three extracted mask images, three artificially filled images, and after the implementation of each gap-filling step for a consecutive three-days cloud mask. D1–D3 refers to the selected first to third days.
Cloud Coverage (%)
D1D2D3
Original Terra17.5013.1264.68
Mask image66.0344.0645.50
Filled image70.4350.4479.59
MOYD61.6333.5262.96
ATF26.9221.9427.37
NSTF0.910.95
Table 6. Quantitative performance evaluation results for a consecutive three-day cloud mask.
Table 6. Quantitative performance evaluation results for a consecutive three-day cloud mask.
R M S E R 2 O E   ( % ) U E   ( % ) S P A E F
D1D2D3D1D2D3D1D2D3D1D2D3D1D2D3
MOYD00.060.1910.970.9300.150.3500.347.820.940.900.95
ATF0.070.130.080.970.930.950.740.850.470.151.801.430.900.970.84
NSTF0.150.130.120.900.880.934.392.762.920.231.121.570.740.730.70
Final images0.090.100.110.960.930.941.581.151.270.150.932.410.870.880.80
Table 7. The effect of the land cover types on the four validation indices of OA , Bias , MU , and MO .
Table 7. The effect of the land cover types on the four validation indices of OA , Bias , MU , and MO .
Land Cover TypeSCF ImagesL-CL-LL-SS-CS-LS-SCloud (%) OA   ( % ) Bias MU   ( % ) MO   ( % )
Cropland
(22)
MOD33,85453,86594020,83124910,18345.6053.410.360.381.44
MYD35,81352,18066620,84322710,19347.2452.010.350.361.05
MOYD27,58959,971109917,49432813,44137.5961.220.470.441.47
ATF902277,985165210,71563819,91016.4681.630.690.641.65
NSTF085,52631330137829,885096.241.051.152.61
Grasslands and shrubs
(11)
MOD18,49923,806142110,6661274429548.6446.870.354.144.61
MYD19,92922,97482312,0341424277753.3142.950.225.092.95
MOYD15,32826,790160895641715495641.5152.940.44.894.59
ATF567735,727232261862664738519.7871.90.65.544.83
NSTF040,21635100423611,999087.080.966.965.85
Urban and built up
(9)
MOD14,31721,1747510,004142334749.5849.980.250.570.3
MYD15,10220,4343011,038248220753.2846.150.161.080.13
MOYD11,70823,765939134276408342.4856.760.310.980.33
ATF405831,3511576688506629921.9076.740.481.320.41
NSTF034,9146520139612,097095.830.952.851.33
Forests
(3)
MOD38617171253713134144946.3252.710.281.530.28
MYD41186923163831168129748.6150.270.2520.19
MOYD30707954333232193187138.5460.080.361.920.33
ATF73810,267522093335286817.3180.320.552.480.38
NSTF010,92912805244772096.010.933.20.78
Unused
(4)
MOD618394701194069515144847.0250.070.264.461.03
MYD65119181803938517157747.9249.340.274.550.7
MOYD488210,7511393443632195738.1858.260.354.691.03
ATF165413,9231952315919279818.2076.690.55.151.09
NSTF015,268504014214611091.170.856.522.31
Note: The figures in parentheses represent the number of in situ SD observation stations. S, L, and C represent snow, land, and cloud pixels respectively. the first letter indicates in situ SD observation and the second letter represents MODIS observation.

Share and Cite

MDPI and ACS Style

Hou, J.; Huang, C.; Zhang, Y.; Guo, J.; Gu, J. Gap-Filling of MODIS Fractional Snow Cover Products via Non-Local Spatio-Temporal Filtering Based on Machine Learning Techniques. Remote Sens. 2019, 11, 90. https://0-doi-org.brum.beds.ac.uk/10.3390/rs11010090

AMA Style

Hou J, Huang C, Zhang Y, Guo J, Gu J. Gap-Filling of MODIS Fractional Snow Cover Products via Non-Local Spatio-Temporal Filtering Based on Machine Learning Techniques. Remote Sensing. 2019; 11(1):90. https://0-doi-org.brum.beds.ac.uk/10.3390/rs11010090

Chicago/Turabian Style

Hou, Jinliang, Chunlin Huang, Ying Zhang, Jifu Guo, and Juan Gu. 2019. "Gap-Filling of MODIS Fractional Snow Cover Products via Non-Local Spatio-Temporal Filtering Based on Machine Learning Techniques" Remote Sensing 11, no. 1: 90. https://0-doi-org.brum.beds.ac.uk/10.3390/rs11010090

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop