Reducing cloud obscuration of MODIS snow cover area products by combining spatio-temporal techniques with a probability of s ow approach

Satellite remote sensing can be used to investigate spatially distributed hydrological states for use in modeling, assessment, and management. However, in the visual wavelengths, cloud cover can often obscure significant portions of the images. This study develops a rule-based, multistep method for removing clouds from MODIS snow cover area (SCA) images. The methods used include combining images from more than one satellite, time interpolation, spatial interpolation, and estimation of the probability of snow occurrence based on topographic information. Applied over the upper Salt River basin in Arizona, the method reduced the degree of cloud obscuration by 93.8 %, while maintaining a similar degree of image accuracy to that of the original images.


Introduction
Water in the Southwestern United States is a scarce resource, requiring efficient management to meet the growing demands of a rapidly growing population.Winter precipitation is particularly important because cooler temperatures allow the accumulation and persistence of snowpack at higher elevations (Jacobs et al., 2005), with spring snowmelt providing inflow to the reservoirs used for water supply and hydropower generation.In fact, even though only 39 % of winter precipitation in the Salt and Verde watersheds falls as snow (Serreze et al., 1999), snowmelt accounts for up to 85 % of the surface water supply for the Phoenix metropolitan area (Hawkins, 2006).
Similarly, snowmelt has been shown to provide 40-70 % of groundwater recharge at several study sites in the Southwestern United States (Earman et al., 2006).This is significant because ∼ 40 % of Arizona's water supply is taken from the underlying aquifers (Megdal, 2004).However, with temperatures in the Southwestern United States projected to rise by 1 • C or more over the next hundred years, the extent and persistence of snowpack is threatened.The ability to estimate and monitor the evolution of snowpack is therefore extremely important.
In this paper we focus on snow cover area (SCA) products derived from spectral imagery collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor aboard two satellite platforms (Terra and Aqua, respectively).These satellite platforms follow the same orbit within 3 h of each other (therefore providing two observations of SCA per day), and have a higher spatial resolution (500 m, 1 day) than other available products such as the Advanced Very High Resolution Radiometer (AVHRR), the Geostationary Orbiting Earth Satellite (GOES) and higher temporal resolution than Landsat.However, the accuracy of the MODIS SCA product is affected by several factors, including land cover type, snow depth conditions, and the presence of cloud cover (Justice et al., 1998;Hall and Riggs, 2007;Zhou et al., 2005;Klein and Barnett, 2003;Bitner et al., 2002;Simic et al., 2004;Tekeli et al., 2005).
Several techniques developed for reducing the cloud cover noise in SCA images have been reported in the literature (see Sect. 2).In this work we investigate the performance V. López-Burgos et al.: Reducing cloud obscuration of MODIS snow cover area products of several such techniques over the Salt River basin in Arizona, and develop a hybrid rule-based, multistep method that computes the probability that a cloudy pixel is underlain by snow.Section 2 reviews the literature regarding snow cover imagery and cloud obscuration, Sect. 3 describes the study region and data used, Sect. 4 discusses three existing and one new method for cloud removal from snow cover images, Sect. 5 reports the results of testing these methods over the Salt River watershed in Arizona, and Sect.6 evaluates the results using ground truth data.Finally, in Sect.7 we discuss the results and make suggestions for future work.
2 Review of the literature

Sensors providing snow cover area imagery
Snow cover maps have been produced since 1966 using remotely sensed data.The most widely used satellites/sensors have been the Landsat thematic mapper (TM) and enhanced thematic mapper (ETM+), AVHRR, the GOES, and the Terra and Aqua MODIS sensors.The Landsat TM (launched in 1982) and the enhanced TM (launched in 1999) provides images at 30 m spatial resolution and a revisit time of 16 days.Although providing superior spatial resolution, the limited temporal resolution limits the applicability of Landsat images for long-term snow cover mapping and modeling studies (Rango, 1985), making it unsuitable for studies in the Southwestern United States, where the snow cover is sparse and ephemeral, and where the snow can melt away in less than 16 days.Nonetheless, Landsat images have been used to evaluate MODIS snow cover products at river-basin-scale areas (Justice et al., 1998) and as ground truth in accuracy studies of MODIS in conjunction with ground measurements (Hall and Riggs, 2007).
The AVHRR sensor, first launched in 1979 aboard the NOAA polar orbiting satellite, provides snow products (Carroll et al., 2001) with higher temporal frequency than Landsat for the continental US and Canada.The AVHRR sensor provides images a spatial resolution of 1 km and a revisit time of 12 h (one daytime and one nighttime pass).While its lower spatial resolution makes it difficult to use for snow mapping in small basins (Schmugge et al., 2002), studies on a 572.9 km 2 basin in the Pyrenees mountains of Spain found that the correlation between AVHRR and MODIS snow maps to be on the order of 0.8-0.9, even in smaller subbasins with areas of ∼ 8.3 km 2 (Gomez-Landesa and Rango, 2000).Similar comparisons in the Columbia and Missouri River basins showed that, on average, MODIS SCA images classified fewer pixels as cloud than the AVHRR images used in the National Operational Hydrologic Remote Sensing Center (NOHRSC) products, indicating it to be a better resource for snow mapping in forested watersheds (Maurer et al., 2003).
The (MODIS) sensors, launched in 1999 and 2002 onboard the Terra and Aqua Earth Observing System satellites, respectively, are designed to complement the Landsat-7 satellite in observing and monitoring Earth system changes (Justice et al., 1998).Daily, and 8-day composite, snow products (Hall and Riggs, 2007) are available globally at no cost and in a variety of resolutions and projections via the National Snow and Ice Data Center webpage (http:// nsidc.org/data/modis/index.html).Due to its higher spectral and spatial resolution it has become popular as an alternative to AVHRR-based snow maps.However, the accuracy of MODIS snow maps has been found to vary with land cover type (Justice et al., 1998).The most frequent errors are caused by difficulties in discriminating between snow and cloud and in the mapping of very thin snow (Hall and Riggs, 2007).Nevertheless, despite propagation of errors from the daily products, errors of commission (mapping pixels as snow where there is no snow) have been found to be very low in the 8-day composite products (Hall and Riggs, 2007).
In general, studies have found the accuracy of MODIS snow products to be persistently good (Bitner et al., 2002;Klein and Barnett, 2003;Simic et al., 2004;Tekeli et al., 2005;Zhou et al., 2005;Hall and Riggs, 2007).Klein and Barnett (2003) compared MODIS snow-cover products with NOHRSC maps, found that "comparisons . . .over the snow season show good overall agreement with accuracies of 94 % and 76 % for MODIS and NOHRSC, respectively".Zhou et al. (2005) conducted an evaluation of MODIS SCA products over the upper Rio Grande basin and reported statistically significant correlations when assessed against streamflow and SNOTEL (SNOw TELemetry) measurements.Hall and Riggs (2007) later found that MODIS fails to map snow at depths less than 4 cm.However, cloud obscuration was reported to be a major cause of reduced accuracy in all studies.

The problem of cloud obscuration
Overall, MODIS SCA products have been found to be useful for reducing uncertainty in knowledge of the extent and amount of snow.However, cloud cover has been found to significantly impact the usefulness of these data and to cause problems for assimilation into hydrological models.For example, McGuire et al. (2005) and Andreadis and Lettenmaier (2006) reported only being able to use SCA images for days where cloud obscuration was less than 20 % of the grid cell, while Rodell and Houser (2004) used 6 % as the threshold for minimum visibility.Consequently, several investigations into techniques for removing cloud cover from SCA images have been reported (e.g., Lichtenegger et al., 1981;Seidel et al., 1983;Molotch et al., 2004;Parajka and Blöschl, 2008;Gafurov and Bárdossy, 2009).
Of course, the cloud obscuration problem is not intrinsic to MODIS SCA products due to the fact that sensors for the visible portion of the electromagnetic spectrum cannot see through clouds and also due to the similarities in the spectral signatures of snow and clouds.The visible and thermal bands can be used to discriminate between snow and clouds (see reviews by Lucas and Harrison, 1990;Klein et al., 1998Klein et al., , 2000;;Riggs and Hall, 2002;Schmugge et al., 2002) and passive microwave can be used to infer snow extent, depth, water equivalent, and state (Schmugge et al., 2002).
Various studies have attempted to reduce the degree of cloud obscuration in SCA image products.Lichtenegger et al. (1981) and Seidel et al. (1983) used elevation, slope, exposure, and brightness information from a digital terrain model (DTM) to extrapolate snow cover from cloud-free to cloud-covered areas in digital Landsat multispectral scanner data, assuming that for each elevation zone, the regions with equivalent exposure and slope angle carry the same amount of snow.Molotch et al. (2004) filtered NOHRSC SCA maps (based on AVHRR/GOES data) using gridded positive accumulated degree days (ADD) and AVHRR-derived binary SCA to obtain a threshold for defining snow cover in the Salt and Verde rivers in Arizona.They reported temperature data to be helpful for estimating snow extent beneath clouds and to thereby improve spatial and temporal continuity of SCA and SWE products.Parajka and Blösch (2006) reported that although MODIS SCA had an average accuracy of 95 % when compared against snow data at 754 climate stations in Austria, clouds covered fully 63 % of the region (even worse in winter when interest in the snow product is higher).They tested a technique for cloud removal that combined MODIS Terra and Aqua data, used majority classification from the eight nearest pixels, and applied a 1-to 7-day temporal window.The approach was found to be remarkably effective at cloud reduction, although accuracy was found to decrease in an almost linear fashion with the application of each filtering technique (Parajka and Blöschl, 2008).Gafurov and Bárdossy (2009) applied the same techniques (with modifications) plus three additional filters -snow transition elevation, spatial combination of four neighboring pixels, and a time series of each pixel over an entire year -to the Kokcha Basin in Afghanistan.In a synthetic data evaluation they achieved complete removal of cloud cover and an overall accuracy of 91.49 %.

Study region and data used
The MODIS Terra sensor provides SCA map products having the best spatio-temporal resolution and accuracy.However, cloud obscuration reduces the value of these images and complicates the images use for data assimilation into models.In this study we develop and test a method for cloud removal from MODIS SCA images over the Salt River basin in central Arizona.

Study area
The upper Salt River basin (Fig. 1), having a drainage area of 11 152.5 km 2 , is a major source of surface water for the Phoenix metropolitan area.On an average annual basis, the precipitation varies spatially from ∼ 400 to 1200 mm yr −1 , with largest precipitation totals at high elevations in the watershed (PRISM Climate Group, 2006a, b), and runoff is 71 mm yr −1 .Annual minimum temperature varies spatially from −13.9 to 3.9 • C, and maximum temperature varies from 17.2 to 39.4 • C (PRISM Climate Group, 2006a, b).The elevation ranges between 674 m and 3472 m a.s.l., and land cover types are primarily ponderosa pine (65 %), chaparral (26 %), pinyon pine-juniper (10 %), and desert grassland (Rinne, 1975).Winter precipitation is of paramount importance since snowmelt can account for up to 85 % of the usable water (Hawkins, 2006), and so snow accumulation and ablation are very important to water resources management.
The local water and power utility, the Salt River Project (SRP), currently relies on sporadic helicopter flights to verify snow cover extent.In addition, daily snowpack data from four SNOTEL locations at higher elevations in the eastern part of the watershed (where snow water storage is greatest; Fig. 1) are used for making streamflow forecasts, but this can lead to erroneous underestimates of streamflow since the point measurements are not representative of the areal pattern of snow accumulation.Consequently, the utility is interested in using remotely sensed SCA imagery to improve the accuracy of their forecasts.

MODIS data
The MODIS SCA products used in this research (MOD10A1 and MYD10A1) are snow cover maps, created from spectral images obtained by the Terra and Aqua satellites.To generate these maps, NASA applies several algorithms based on the normalized difference snow and normalized difference vegetation indices (NDSI and NDVI, respectively), as well as a cloud mask to distinguish snow and cloud pixels (Hall and Riggs, 2007).The cloud mask algorithm uses 14 of the 36 MODIS bands in 18 cloud spectral tests that start with an initial guess of whether snow or cloud is being viewed so that appropriate processing paths and tests can be applied based on surface type, geographic location, and ancillary data.The result -"cloudy", "clear" or "probably clear" -is reported for each pixel (Ackerman et al., 1998).
The products are provided on a sinusoidal grid projection divided into 10 • × 10 • tiles (approximately 1200 × 1200 km at the equator) (Wolfe et al., 1998).The "h08v05" tile includes Arizona, Utah, Southern California, and parts of Northern Mexico.The data are provided in Hierarchical Data Format (HDF), as binary SCA, fractional SCA (FSCA), and snow albedo for each day, through the National Snow and Ice Data Center (NSIDC) webpage (http://nsidc.org/data/modis/index.html).We used SCA data from 1 October 2004 to 31 May 2005 (8 months).
ArcGIS9 ArcMap Version 9.3 was used to compute projection changes and for processing geographic information system data (e.g., digital elevation model).All imageprocessing computations were carried out using MATLAB 7.5.0.
Because the eventual goal is to assimilate the SCA images into a distributed hydrological model, Fig. 2 presents an example in which the SCA data were upscaled from its native resolution of ∼ 500 m (51 375 pixels) to grids of ∼ 12.5 km (116 grids).The upscaled data exhibited irregular behavior, with many grids going from 100 % FSCA to 0 % FSCA from one day to the next, not properly reflecting the gradual onset of snowmelt, as evident in the snow water equivalent (SWE) time series for a SNOTEL station located inside the grid (Fig. 2).Such erratic behavior illustrates the large extent to which clouds cover the watershed during the study period.
Figure 3 illustrates the number of days with snow and cloud over the study period.Figure 3a shows a pixel map of the number of "days of snow" across the basin between 1 October and 31 May, constructed from the raw SCA data, clearly indicating that snow is being accumulated in the right places (higher elevations on the eastern side).To complement this, Fig. 3b shows a pixel map of number of days with clouds over the same period; note that this number ranges from 70-120 days (29-49 %) of the total 24 days in the study period.Overall, clouds cover 39 % of the SCA image pixels, mostly during periods of active snowfall.Further, it can sometimes be several days before the cloud cover clears away, during which time some of the snow at lower elevations can melt away and therefore not be accounted for.
To simplify our analysis, the MODIS SCA images were first pre-processed to classify each pixel as belonging to only one of the following five categories -cloud (new code 1/original codes 11 and 50), error (new code 2/original codes 0, 254 and 255), snow (new code 3/original code 200), land (new code 4/original codes 25, 37, 39 and 100), and no decision (new code 5/original code 1).In addition, three new categories were added -missing day (new code 6), corrected snow (new code 7), and corrected land (new code 8).

SNOTEL data
SNOTEL sites are automated snow telemetry stations that record real-time SWE data as well as snow depth, precipitation, and min.and max.air temperature in their standard configuration.SNOTEL sites are spread throughout mountainous areas of the western US (more than 750 stations) and managed by the Natural Resources Conservation Service (Serreze et al., 1999).Twenty-one of these stations are in the state of Arizona, but only four fall within the boundary of the study area (see Fig. 1).

Methods applied for cloud removal
Based on the review of previously reported research, we investigated three methods for their effectiveness in cloud removal from MODIS imagery over our study area: 1. Combining Terra (MOD10A1) and Aqua (MYD10A1) imagery SCA products.
3. Nearest neighbor spatial interpolation informed by elevation and aspect.The aspect of a pixel is the exposure of the terrain represented by that pixel, in other words, the direction in which the slope of the terrain faces (e.g., north, south, northeast, etc.) when it is not flat.
In addition we developed and tested a fourth method: 4. Probability of snow estimation, via logistic regression informed by elevation and aspect.Each of these methods was first tested separately to assess the effectiveness of its performance.Subsequently, we implemented and tested a hybrid approach based on sequential application of the aforementioned methods.

Terra/Aqua combination
Combination of MOD10A1 and MYD10A1 SCA images takes advantage of the short (3 h) time interval between the two observations and the transience of clouds during that time.The principal assumption is that snow conditions remain essentially constant (no snow falls or ablates) during that period.If a pixel on the one image is classified as cloud (or error) but is reliably classified as land or snow cover on the other, then the problem pixel is reclassified appropriately.This method was adopted based on the experience of Parajka and Blösch (2008).The algorithm took approximately 12 min to execute on an Intel Pentium D 2.80 GHz computer.

Time interpolation
The time-interpolation method employs three temporal windows applied sequentially. .This approach makes the assumption that snow cover remains essentially constant during the spans of the different temporal windows.This assumption may be considered reasonable because the probability of snowmelt tends to be lower during cloudy days.This method was adopted based on the experience of Gafurov and Bárdossy (2009), who modified the temporal windowing approach applied by Parajka and Blösch (2008).The algorithm took approximately 3.2 min to execute on an Intel Pentium D 2.80 GHz computer.

Nearest neighbor spatial interpolation
The nearest neighbor spatial interpolation method uses information from the eight pixels surrounding (edge and diagonal) a cloudy pixel.If any neighbor pixel has snow cover, is at a lower elevation, and has the same aspect, then the cloudy pixel is classified as snow.Once the entire image has been processed, a similar logic is applied to any remaining cloudy pixels -if any neighbor pixel is land covered, at a higher elevation, and has the same aspect, the cloudy pixel is classified as land.
A 30 m resolution digital elevation model (DEM) of the area was obtained from the United States Geological Survey (USGS) and processed using the aspect tool of the spatial analyst extension in ArcGIS 9.3.Each pixel was assigned an aspect value from 0 • (due north) to 360 • (again due north) measured counterclockwise, or a value of −1 if the slope is flat.The pixels were classified accordingly as north, south, east, west, northeast, southeast, northwest or southwest facing.A resampling routine was then run to upscale the "elevation" and "aspect" to the 463.32 m resolution of the SCA images to allow comparison between the three layers of information.The resampling rules "bilinear interpolation" and "majority of pixels" were applied, taking into consideration the nature of the data.An elevation raster is a continuous surface and therefore it is most appropriately resampled using a bilinear interpolation, where the elevation values of the four nearest 30 m input cell centers (used to create the new 463.32 m cell) are averaged.This average is weighted, taking into consideration the proximity of the four nearest 30 m input cell centers to the new 463.32 m cell center.For the aspect raster, the "majority of pixels" rule was used.This rule determines the value of the new cell by choosing the most frequent values of the 30 m input cell centers that are within the new 463.32 m cell.
This method was modified from Gafurov and Bárdossy (2009), who used only elevation information.We included aspect information because topographic controls like elevation, aspect, and slope (not used here) can significantly influence energy exchange and melt across a watershed by modifying the exchange of direct-beam and diffuse shortwave radiation and longwave radiation (Dewalle and Rango, 2008).The algorithm took approximately 8 h to execute when run by itself on an Intel Pentium D 2.80 GHz computer.

Locally weighted logistic regression
The locally weighted logistic regression (LWLR) method uses relationships between the spatial and topographic attributes of pixels surrounding a cloudy pixel to estimate the "probability of snow occurrence" (PSO, Eq. 1).This method is adapted from Clark and Slater (2006), who used precipitation observations at sparsely located meteorological stations, as well as spatial maps of elevation, latitude and longitude, to estimate daily precipitation totals across complex terrain in western Colorado.Here we used elevation and aspect as the explanatory variables.
To estimate PSO at a cloudy pixel, the LWLR method weights the information from neighboring pixels inversely with distance, and fits the data to a logistic curve.The following equations are used to calculate the PSO of a cloudy pixel (PSOicloud): and where Z icloud is a vector of elevation and aspect information for the cloudy pixels indexed as icloud, X is vector of elevation and aspect information for the non-cloudy pixels, β is a vector of parameters, Y is a 0-1 vector indicating snow occurrence or not on the non-cloudy pixels, W is a diagonal matrix of weights to be assigned to each non-cloudy pixel, V is a diagonal matrix of the variance associated with the estimate of snow occurrence at each noncloudy pixel, π ipix indicates the PSO at each noncloudy pixel, d ipix is the distance from a non-cloudy pixel to the cloudy pixel, and MAXD is a coefficient specifying the window size used around the cloudy pixel (see Clark andSlater, 2006, andLoader, 1999, for details regarding implementation).
In regression analysis one makes use of one or more known independent explanatory variables X to predict the value of an unknown dependent variable Y .In logistic regression the model fitted to the data is a logistic function or logistic curve (Eq. 1) that represents the probability of occurrence of a "categorical" [0, 1] variable Y given X, and the regression coefficients are estimated using maximum likelihood estimation via an iterative process (Eqs.2, 3, and 4).Note that although X can take any value on [−∞, ∞], the value of Y is confined to the range [0, 1]; in our case Y represents the probability of snow occurrence.The regression is computed using information from the pixels that are not obscured by clouds, and then used to estimate probability of snow for the cloudy pixels in the image.In addition, the regression is computed and applied in a locally weighted manner, meaning that pixels closer to the cloudy pixel have more weight than pixels those that are further.The weights are calculated via Eq.( 3).
For this study we tested different values for the window size MAXD (from 5 to 45 pixels) for their ability to provide statistically robust results, reduce significant numbers of cloudy pixels, and require reasonable computational time.Reliability was evaluated using the Brier score (BS) verification statistic (see Wilks, 2005) and Clark and Slater (2006) computed from the joint distribution of the LWLR forecast probabilities and the observed snow/land pixels in the image.Although the best results (not shown) were obtained using MAXD = 45, we instead chose MAXD = 30 due to the relatively high performance achieved using only 1/3rd of the computer time (22 h as opposed to 67 h on a 2.66 GHz dualcore Intel Xeon computer).
Note that LWLR does not, by itself, automatically reclassify cloudy pixels, but only provides an estimate of the probability that the pixel is actually snow or land.Therefore, a minimum probability threshold must be selected to convert the probability to a binary outcome.To select this threshold, we conducted a sensitivity analysis by varying the threshold from 0 to 1 in steps of 0.025 and chose the threshold value that minimized the sum of the conditional probabilities of commission and omission errors computed over all non-cloudy pixels inside the window (for details see López-Burgos, 2010).The threshold was separately selected for each cloudy pixel to which LWLR was applied.Figure 4 shows an example of the transition from original to corrected MOD10A1 image (26 November 2004), along with maps of the estimated PSO and the values of the thresholds selected.

Results for each cloud removal method applied independently and in sequence
Each of the four methods was first tested separately to assess the effectiveness of its performance.Subsequently, a sequential approach was tested, in which the methods were applied in sequence.Two summary statistics are used to indicate the change achieved by implementation of each method:

Hydrol
where x i represents codes 1 through 6, and t represents the time period for which the statistic was calculated (e.g., year or month).
We first discuss the results of applying each method separately.The change in percentage cover achieved over the entire study time period by each method is shown in Table 1 for each category.Note that 39 % of the image pixels were initially classified as cloudy.
In brief: 1.The Terra/Aqua combination method had a percent change of cloud cover, error, and no-decision pixels of −23 %, −97 %, and −3.7 %, respectively.This translates to a cloud cover change from 39 % to 30 %.These results are similar to the cloud removal achieved with the same method by Parajka and Blöschl (2008), Xie et al. (2009), and Gao et al. (2010).Furthermore, the 2 days missing in the Terra images time series were not missing in the Aqua images and so a more complete time series was achieved.In addition, the combination of MOD10A1 and MYD10A1 increased the snow-covered pixels and the no-snow pixels by 34.5 % and 14.6 %, respectively (+1.55 and +8.13 snow and land cover).On a monthly basis, this approach works better during the fall and spring months versus the winter months with a maximum cloudy pixels reduction in May (−36 %) and a minimum cloud reduction in January (−16 %) (Fig. 5).
The monthly distribution of differences in cloud removal are similar to those of Parajka and Blöschl (2008), Xie et al. (2009), Gafurov and Bárdossy (2009), and Gao et al. (2010), while differing on the amount of cloud cover removed on each site.This might be due to differences in climate, topography, study area (km 2 ), and cloud dynamics between the watersheds during the different time periods used.
2. The time interpolation method had a percent change of cloud cover, error, and no-decision pixels of −43 %, −95 % and −46 %, respectively, with a cloud cover change from 39% to 22% while increasing the snowcovered pixels and no-snow pixels by 33 % and 28 %, respectively (+1.51 and +7.25 snow and land cover).These results are similar to the 1-day temporal filter used by Parajka and Blöschl (2008) and significantly less than their 7-day temporal, though the results are not strictly comparable since they applied the temporal filters to the already improved images by way of Terra/Aqua combination.The same applies to the cloud removal results obtained by Gao et al. (2010), which also used improved images as the inputs for their temporal filters.Xie et al. (2009) achieved a cloud cover reduction from 39.5 % to 6.0 %, but the amount of days used for each image and the weighted mean function used are not clear.This method removed more bad pixels than Terra/Aqua, but it does not give information on missing days.It also performed better during the fall and spring months than in the winter months, with a maximum cloud reduction in May (−84 %) and a minimum in February (−17.6 %).However, as shown on Fig. 5, cloud reduction in November was not as good as in October and December due to higher percentage of cloud cover in November.As is the case with the Terra/Aqua combination, this method performed less strongly during winter months due to cloud persistence beyond the temporal window used to improve the images.
3. The nearest neighbor spatial interpolation method is the least effective of the four methods.The percent changes on cloud cover, error, and no-decision pixels were only −5 %, −9 %, and −0.03 %, respectively, with a cloud cover change from 39 % to 37 %.These results are similar to that of Parajka and Blöschl (2008).Their spatial interpolation method reduced the cloud cover from 52 % to 46 %, although they applied this algorithm to the images after having improved them with the Terra/Aqua combination.Again, this method does not give information on missing days.Furthermore, this method does not show the seasonal difference in cloud reduction power of the three previous methods; its effect is almost constant throughout the time period.
4. The locally weighted logistic regression (LWLR) method achieved the highest effectiveness at the cloudcovered and no-decision pixels (−62 % and −60 %, respectively), but was significantly less effective than either Terra/Aqua combination or time interpolation at reducing the number of "error" pixels (−44 %).This method also had the highest increase in snow-covered pixels (96 %) and no-snow pixels (35 %) and shows a similar seasonal pattern of cloud removal with the maximum cloud reduction in April (−90 %) and a minimum in February again (−50 %).The change in cloud cover was from 39 % to 15 %.Notwithstanding, it does not take care of the missing days problem either.
Figure 5 presents the change in pixel classification for each month in the study period.Figure 5a shows how the percentage of pixels classified as cloudy in the original MODIS SCA images varies for each month, from a low of 24 % in April to a high of 69 % in February.Figure 5b-d shows the % change in clouds, snow, and land pixels achieved by each method for each month.Nearest neighbor spatial interpolation is consistently poor at cloud removal, while logistic regression is the most effective.For removing no-decision pixels, the Terra/Aqua and time interpolation algorithms perform best.All methods increased the amount of snow-covered pixels more than they increased the amount of land pixels -a reasonable result since the areas with the most consistent cloud cover throughout the study period are the areas of snow accumulation, mainly at higher elevations, and the period with more cloud cover is during the winter months when snow accumulates (Fig. 3).
The results above indicate that each method has different strengths and weaknesses.To synergistically exploit the strengths of all four methods, we next applied the methods in sequence (in the order Terra/Aqua combination, time interpolation, spatial interpolation, and locally weighted logistic regression), at each step retaining the results from the previous one.Table 1 shows that the sequential approach achieved a very high degree of cloud cover, error, and no-decision pixel removal (−94 %, −99 % and −64 %, respectively), while increasing the number of snow and land pixels by ∼ 154 % and ∼ 54 %, respectively -a considerable increase in the amount of SCA indicated by the images.Figure 6 shows the progressive improvement obtained by sequential application of the methods in terms of percent cover (i.e., how the cloud/snow/no-snow covers change with each step) and how this improvement distributes across months.
Overall, using the sequence the cloud cover was reduced from 39 % to 30 % using the Terra/Aqua combination, then to 14 % with the time interpolation, to 13 % with the nearest neighbor spatial interpolation and down to 2 % with the LWLR.On a monthly basis the cloud cover was reduced to less than 10 % for all months with the most cloudy month (February) having a final cloud cover of 9 % and October, April, and May having a final cloud cover of less than 1 %. Figure 6c and d look like a typical snow curve and an upsidedown snow curve.The increase in snow-covered pixels follows the expected monthly snow distribution except that more snow was added to the month of January than in February.This is likely due to the difficulty of removing clouds in February as it had a higher percentage of cloud cover.The increase in no-snow pixels also shows the expected monthly distribution with the fall and spring months showing a higher increase.Moreover, only 5 % of the days were left with a cloud cover greater than 10 % (not shown).These days are distributed between December and March with the majority being in February, the cloudiest month, and are grouped in doubles or triples (i.e., 2 or 3 consecutive days with heavy cloud cover).Therefore, it makes sense that the algorithms were not able to reduce the cloud cover to less than 10 %. Figure 7 shows the significant change obtained by the sequential method for the MODIS SCA image of 18 February 2005.The changes achieved by the hybrid approach are clearly substantial and in marked contrast to the changes achieved by each method applied independently.

Evaluation of accuracy of the results
Finally, we assess the accuracy of the results by comparing the SCA images with data from four available SNOTEL sites located in the mountainous zones of the Salt River basin (Fig. 1).Although the SNOTEL sites are effectively pointscale measurements, their pixel locations were localized using ArcGIS 9.3 by converting the station's point shapefile to a raster with the same extent, resolution, and projection of the images.This raster was then converted to an ASCII file for processing.Several evaluation statistics were calculated for (a) the original MOD10A1 and MYD10A1 images, (b) each cloud removal method applied separately, and (c) the sequential cloud removal approach.For each SNOTEL site pixel, the hit, false alarm, miss, and correct rejection rates were computed for each day of the available time series.A simple approach to evaluation would be to consider each SNOTEL site pixel to be snow covered if the corresponding SNOTEL station has recorded measurable SWE (SWE > 0).This assumption would imply that each station measurement, corresponding to an area of approximately 9 m 2 (snow pillow standard size is 3 × 3 m), is representative of conditions across the entire 500 × 500 m (250 000 m 2 ) image pixel.However, this assumption can be poor for several reasons -an important one being that the MODIS sensor fails to map snow when snow depths are less than 4 cm (Hall and Riggs, 2007).
Based on an average snow density of 0.3621 g cm −3 (average snow density for one of our SNOTEL stations during WY 2005), 4 cm of depth corresponds to approximately 1.4478 cm of SWE.Therefore, to establish a more accurate basis for evaluation, a sensitivity analysis was performed to find a threshold value of recorded SWE above which the pixel could be considered to be snow covered.To do this the SNOTEL station raster was used as the observed ground truth and the MODIS SCA images as the modeled forecast.The conditional probabilities of hits (observed = snow/forecast = snow), false alarms (observed = no snow/forecast = snow), misses (observed = snow, forecast = no snow) and correct rejections (observed = no snow/forecast = no snow) were computed for 20 different threshold SWE values (from 0 to the maximum value of SWE = 49.53 cm recorded at the stations over the observation time period), and an optimal threshold (SWE = 2.61 cm) was selected that minimized the sum of the conditional probabilities of misses and false alarms, while maximizing the sum of the conditional probabilities of hits and correct rejections.
The following evaluation statistics recommended by Wilks (1995) were then computed: and where a = number of hits, b = number of false alarms, c = number of misses, d = number of correct rejections, and n = a+b+c+d.The proportion correct (PC) and threat score (TS) are accuracy statistics, while B is a measure of bias, the false alarm ration (FAR) is a measure of reliability and the hit rate (H ) is a measure of discrimination.The proportion correct is a good measure of accuracy if the event (snow) and nonevent (no snow) occurred with equal frequency (i.e., 50/50).A completely accurate estimator will achieve a PC = 1 (b = c = 0), while a completely inaccurate estimator will have a PC = 0 (a = d = 0).This accuracy measure would be most useful during the winter months.During the fall and spring months, snow cover is less frequent than no-snow cover, thus the threat score is a better accuracy statistic for those months since it is good for situations in which the event to be forecasted (snow occurrence in this case) occurs less frequently than nonoccurrence.An accurate estimator will achieve a TS = 1 (b = c = 0), while an inaccurate estimator will have a TS = 0 (a = 0).An unbiased estimator will achieve B = 1 (b = c = 0); B > 1 indicates that snow is estimated more often than observed, while B < 1 indicates that snow is estimated less often than observed.A good estimator will also achieve an FAR close to 0 (no false alarms) while FAR close to 1 indicates very poor performance (no hits).Finally, a good estimator will achieve a hit rate close to 1 (no misses).
The results are summarized in Table 2.We see the following and Gao et al., 2010) This might be due to differences in evaluation statistics used, the amount of stations with ground truth observations used to calculate the accuracies and their distribution throughout the watersheds, and/or the fact that the other studies included the entire water year including summer where mapping accuracies are high because there is no snow and therefore omission and commission errors are practically non-existent.
3. The time interpolation method achieves the best overall accuracy, and provides consistently better evaluation statistics -better even than those of the original images.This is likely due to the seasonal persistence of snow at the SNOTEL sites (e.g., see Fig. 2), and the result may not be representative of mapping accuracy near the snow line.We will return to this point later in the paper.
4. The other three methods have varying results, with LWLR achieving the worst PC, TS, B, and H statistics.The lower accuracy results for the LWLR might be caused by the size of the window used for the analysis; a smaller window has better accuracy results but removes less clouds (López-Burgos, 2010).The most appropriate tradeoff between these two qualities can be chosen by the user based on the future applications of the final images.The lower accuracy of the LWLR may also be due to the thresholds used to decide if a pixel has snow or not.One way to improve this would be for the user to give more weight to the minimization of the sum of conditional probabilities of commission errors than omission errors (refer to Sect.4.4) since this would give more conservative results for snow cover.It is better to plan for less SWE in the watershed (underestimation) and find that there is more usable SWE than to plan for a higher amount of SWE (overestimation) and find out there was actually less.Development of a method that improves this step remains a topic for future work.
5. The sequential approach is second best in terms of accuracy (PC = 0.85, TS = 0.74), matches the Terra/Aqua combination in terms of bias (0.89), is middling in terms of false alarms (0.10) and second in terms of hit rate (H = 0.81).Table 3 shows how the evaluation statistics changed with each step during the sequence application.
As shown on the table, the first step (Terra/Aqua combination) worsens all the evaluation statistics except for the hit rate, which stayed the same.Following this step the time interpolation was applied and the evaluation statistics improved to values even better than the original Terra image, except for the FAR, but still showed an improvement on this statistic over the Terra/Aqua combination.The spatial interpolation did not cause any significant changes and the LWLR worsened the evaluation statistics of the sequence, but the final outcome is equal or better than the original Terra image.The final FAR statistic is higher than the original image and so is the bias.This is most likely caused by an increase in the commission errors for the FAR and also a decrease in omission errors for the bias.A way of improving or balancing this result was mentioned above.
6. Overall, all methods showed higher omission errors than commission errors.These results are similar to those of Parajka and Blöschl (2008), Xie et al. (2009), andGao et al. (2010).In addition, the final result of the sequence has similar or higher accuracy than the original MODIS images; this result is also similar to the aforementioned papers with respect to the methods they used and final products.

Differences in performance between the cloud removal algorithms
The cloud removal algorithms examined in this study gave differing results with regards to the number of cloudy pixels removed.In terms of number of cloudy pixels removed, LWLR was best, followed by time interpolation, and Terra/Aqua combination in that order.Spatial interpolation provided very little cloud removal performance.When removing clouds, all methods added more snow pixels than land pixels -a reasonable result since clouds tend to concentrate at higher elevations where more snow accumulates and the time period used only includes months when there is snow presence.The methods also share the same seasonal patterns of cloud removal and distribution of added snowcovered and no snow-covered pixels.Overall, the sequential approach achieved a 94 % reduction of cloud-covered pixels (from an overall cloud cover of 38.7683 % to 2.4084 %), which can be considered to be very successful.The results would, of course, be even better if we did not include images from days that were completely covered by clouds and for which the algorithms had little effect.Furthermore, the sequential approach reduced the cloud cover to less than 10 % for all months, with a maximum of 9 % in February followed by 4 % in January and a minimum of 0.15 % in April.It is likely that LWLR is able to remove more cloudy pixels than any other method because it uses information from a window of pixels of area ∼ 775 km 2 around each cloudy pixel (∼ 7 % of the total area of the watershed).This gives the method an advantage over the other algorithms, which only draw upon information at the same pixel or in the 8 pixels neighborhood having an area of 0.86 km 2 .Because the Terra/Aqua combination and time interpolation methods depend on the dispersion of clouds during the time step, cloud cover that persists for several days can confound the methods.In contrast, even if cloud cover persists on the same area for several consecutive days, the LWLR method can use information from non-cloudy pixels at lower elevations or high elevations elsewhere in the watershed.
On the other hand, LWLR is unable to remove cloud cover on days when clouds cover all or much of the watershed.In such situations, the Terra/Aqua combination and time interpolation methods can be effective, provided the clouds do not persist for too many hours or days.Further, the Terra/Aqua combination helps to complete the time series -a positive and useful aspect that none of the other algorithms possess.Taken together, it therefore makes logical sense that a sequential combination of these algorithms provides synergistic effects.Overall, spatial interpolation provides only minor improvements, and its abilities are partially replicated by the LWLR approach.Elimination of this method from the sequence will reduce processing time and have very little overall impact on the results.

Differences in evaluation results between the cloud removal algorithms
Time interpolation provides the best overall accuracy with consistently better evaluation statistics (better than that of the original images).LWLR removes the largest amount of cloud cover but has worse evaluation performance due inherently to the extrapolation of information from surrounding pixels.This may be aggravated by the large window size used.Terra/Aqua combination also reduces accuracy even though the observations are from the same day and the instruments should give similar results.One reason for this may be that the MODIS instrument aboard Aqua has degraded data quality in band 6 (70 % of the band 6 detectors have been identified as nonfunctional), pertaining to the shortwave infrared portion of the electromagnetic spectrum (∼ 1.6 µm) where snow surfaces have low reflectance, and used in the computation of NDSI used for SCA computations.As a consequence, band 6 was substituted with band 7 in an Aquaspecific algorithm to map snow.The 2005 accuracy of the Aqua MYD10A1 images has not yet been assessed (Hall and Riggs, 2007).Overall, sequential application of the methods achieved accuracy similar to that of original MOD10A1 image, with improvements in certain areas.Further, the threat score increased (which means that more pixels were correctly classified as snow), the bias came closer to one (meaning that in general there were relatively more hits than misses and false alarms and underestimation of snow events was reduced), and the hit rate moved closer to one (meaning that more pixels were correctly classified as snow).

Data assimilation applications
While satellite SCA estimates can be assimilated into a hydrologic model to update model states and improve streamflow simulations and forecasts (e.g., see McGuire et al., 2005;Andreadis and Lettenmaier, 2006;Clark et al., 2006), cloud obscuration reduces the value of SCA information by reducing the number of updates and by complicating procedures to interpolate/aggregate SCA estimates from the satellite pixels to the model grid cells.As such, cloud removal procedures provide additional information for assimilation into hydrologic models and simplify data assimilation procedures.
In a data assimilation context, the value of SCA information depends on the seasonality of snow and the availability of SCA images.For example, models and satellites typically have a high level of agreement in the portrayal of SCA throughout much of winter in seasonally snow-covered areas -both are close to 100 % -suggesting that assimilation will be most effective in locations where snow is ephemeral (see Clark et al., 2006, for more discussion of this issue).The key problem therefore is to map the snow line.
Based on this discussion, the use of SNOTEL sites in seasonally snow-covered areas is perhaps a suboptimal data source to evaluate the value of cloud removal methods for data assimilation applications.The persistence of snow at these SNOTEL sites means that temporal infilling methods should perform much better than spatial methods (as evident in the results presented in Table 2) -these site-specific results are not necessarily representative of the value of different cloud removal procedures for SCA assimilation in spatially distributed models.We speculate that spatial methods such as logistic regression may provide superior estimates of the snow line as well as provide an estimate of the observation error needed in more formal data assimilation schemes such as the Kalman filter (Andreadis and Lettenmaier, 2006;Clark et al., 2006) -further work is necessary to directly evaluate how different methods for cloud removal can improve SCA assimilation in spatially distributed models.

Conclusions
The sequential cloud removal approach has the potential to be applied successfully at other locations.Whereas the Terra/Aqua combination, time interpolation, and spatial interpolation methods have been previously tested in watersheds having different topographical and climate characteristics (Parajka and Blösch, 2008;Gafurov and Bárdossy, 2009;Xie et al., 2009;and Gao et al., 2010), the LWLR probability of snow method has not previously been used for cloud removal and SCA estimation and is the only method found so far that uses probabilities in the rules used to remove clouds.It should therefore be subjected to more extensive testing (e.g., how to choose the better thresholds to get binary results).It may also be interesting to examine whether the additional use of slope information and topographic shading as explanatory variables would prove helpful.
In addition, in order to make better comparisons with results of other papers, it is necessary to address the mapping accuracies of the methods on a monthly basis to find out the tradeoffs between cloud removal and mapping accuracies during transitional and stable periods.This step was not included in the original version of this study (López-Burgos, 2010) and therefore could not be added to this paper.However, the seasonal differences in mapping accuracies of the Terra/Aqua combination and time interpolation are expected to be similar to the ones found by Parajka and Blösch (2008) and Gao et al. (2010).In general, mapping accuracies for these methods are higher during the snow-covered periods and lower for the snow accumulation and snow ablation periods.This is due to the sudden change of snow cover from one day to the other during the transition periods.These two other studies also found that increasing the temporal window used in the time interpolation/filter reduces the accuracy of the images.This would be even more important during the transitional periods.
The most closely related study that we could use to give an estimate of the seasonal mapping accuracy of the LWLR is that of Parajka et al. (2010), who used a regional snow-line method for estimating the snow cover from MODIS during cloud cover.This method would be similar because it also uses snow cover information and elevation from cloud-free pixels to decide if there is snow or not on the cloudy pixel.However, this method does not take into consideration the aspect of the pixels, and as the topography and snow regime in Arizona and Austria are different, the mapping accuracies would not be strictly comparable (Parajka et al., 2010).
Our future work will include the assimilation of corrected SCA images into a distributed hydrologic model to reduce the uncertainty of streamflow forecasts.McGuire et al. (2005) and Andreadis and Lettenmaier (2006) have demonstrated that although assimilation of MOD10A1 images into the variable infiltration capacity model can provide favorable results, cloud cover reduces the assimilation frequency.Since removal of cloud cover can result in considerably larger estimates of snow, application of the algorithm developed in this work should help to simplify the assimilation process while improving the model estimates of various hydrological states and fluxes.

Fig. 2 .
Fig. 2. Fractional snow-covered area time series for one grid, overlayed upon a Snow Water Equivalent time series for a SNOTEL station located inside the grid.This plot illustrates the extent to which clouds cover the watershed during the study period with up to 69 % cloud cover in February.

Fig. 3 .
Fig. 3. (a) Pixel map of number of "days of snow" across the basin between 1 October and 31 May, indicating that snow is accumulated at higher elevations on the eastern side; (b) pixel map of number of days with clouds over the same period.

Fig. 4 .
Fig. 4. Example of the transition from original to corrected MOD10A1 image (26 November 2004), along with maps of the estimated PSO and the values of the thresholds selected.

Fig. 5 .
Fig. 5. (a) Shows how the percentage of pixels classified as cloudy in the original MODIS SCA images varies for each month during the study period; subplots (b), (c), and (d) show the corresponding percentage change in clouds, snow, and no snow achieved by each method for each month.

Fig. 6 .
Fig. 6.Progressive improvement obtained by sequential application of the methods and how this improvement distributes across months.

Fig. 7 .
Fig. 7. Example of the sequence results for 18 February 2005; (a) original image, (b) after T/A combination, (c) after time interpolation, (d) after spatial interpolation, and (e) after logistic regression (final result).
1.Both Terra and Aqua images have an overall high fraction of correctly classified snow-covered and snow-free pixels (PC = 0.85, 0.82) but tend to underestimate the occurrence of snow (B = 0.85, 0.91); therefore, commission errors are low.The latter is shown by low FAR values (0.08, 0.15).The H values (0.78, 0.77) are consistent with the B values; both Terra and Aqua images slightly underestimate the occurrence of snow.Overall, the Terra images have better accuracy than the Aqua images and both images show higher omission errors than commission.These results are consistent with those of Parajka and Blöschl (2008), Xie et al. (2009), and Gao et al. (2010) although they used other measures of accuracy.

Table 1 .
"%Cover" change of the MOD10A1 images, achieved by each method for each category for the entire study period.

Table 2 .
Evaluation Statistics for the original images and individual results.The data in bold is the method with the best accuracy results.

Table 3 .
Evaluation statistics for each step of the sequence.