Snow Property Inversion From Remote Sensing (SPIReS): A Generalized Multispectral Unmixing Approach With Examples From MODIS and Landsat 8 OLI

Spectral mixture analysis has a history in mapping snow, especially where mixed pixels prevail. Using multiple spectral bands rather than band ratios or band indices, retrievals of snow properties that affect its albedo lead to more accurate estimates than widely used age-based models of albedo evolution. Nevertheless, there is substantial room for improvement. We present the Snow Property Inversion from Remote Sensing (SPIReS) approach, offering the following improvements: 1) Solutions for grain size and concentrations of light absorbing particles are computed simultaneously; 2) Only snow and snow-free endmembers are employed; 3) Cloud-masking and smoothing are integrated; 4) Similar spectra are grouped together and interpolants are used to reduce computation time. The source codes are available in an open repository. Computation is fast enough that users can process imagery on demand. Validation of retrievals from Landsat 8 operational land imager (OLI) and moderate-resolution imaging spectroradiometer (MODIS) against WorldView-2/3 and the Airborne Snow Observatory shows accurate detection of snow and estimates of fractional snow cover. Validation of albedo shows low errors using terrain-corrected in situ measurements. We conclude by discussing the applicability of this approach to any airborne or spaceborne multispectral sensor and options to further improve retrievals.

Snow Property Inversion From Remote Sensing and Applications [2] ranked snow accumulation, melt, and sublimation among the "Most Important" objectives.
Remote sensing of snow from aerial and spaceborne sensors for runoff forecasts has a long history going back to the 1950s [3]. A longstanding approach applied to multispectral instruments, the Normalized Difference Snow Index (NDSI) that Dozier [4] introduced, is R λ is the apparent planetary reflectance, i.e., the directionalhemispherical reflectance [5] if angular distribution of the reflected radiation is assumed isotropic. VIS identifies a band in a visible wavelength, and SWIR identifies a band in the shortwave infrared. In the original formula for the Landsat 5 Thematic Mapper, the VIS band was TM2 (0.525-0.625 μm) and the SWIR band was TM5 (1.55-1.75 μm); similar spectral bands are used when applied to other sensors. The NDSI was the first spectral index designed specifically to identify snow from space and is still widely used to classify snow cover, for example, the moderate-resolution imaging spectroradiometer (MODIS) algorithm in NASA's Earth Observing System [6]. Snow-covered pixels at resolutions greater than a few meters often contain rock, soil, or vegetation (e.g., Fig. 1), so techniques to account for mixed pixels must be used [7]. For MODIS, a regression-based approach was developed to convert NDSI to fractional (i.e., subpixel) snow-covered area ( f sca ) [8], but this technique showed high root mean square error (RMSE) values [9] so the MODIS data system no longer provides a fractional snow product. A nonlinear (sigmoid) regression of f SCA derived from the high resolution (0.5-2 m) Pléiades constellation versus NDSI calculated from the Sentinel-2 satellites showed RMSE values of 0.25-0.38 in the Pyrenees, with the least certain results in regions of rugged topography [10].
The NDSI approach is attractive because of its simplicity, but using only two bands does not take advantage of all available spectral information and provides no information on other relevant snow properties such as grain size or the presence of light absorbing particles (LAPs).
Spectral mixture analysis [11] has been successfully used to map snow cover with multispectral sensors such as the Landsat Thematic Mapper [12] or MODIS [13] and spectroscopic This  airborne instruments such as airborne visible and infrared imaging spectrometer (AVIRIS) [14]. An independent study identified spectral analysis as providing more accurate results than regression against NDSI [15].
Analyses of the spectra also provide information about snow grain size [16], [17] and LAP [18]. Retrievals of albedo from spectral unmixing have been shown to be more accurate for estimating broadband snow albedo when compared with commonly used models based on aging [19].

II. SPECTRAL UNMIXING APPROACHES
Spectral mixture analysis was created to identify chemical components in mixtures [20] and has since been applied to many studies including subpixel vegetation mapping [11] and geologic classification of rocks [21]. The reflectance R λ of a heterogeneous instantaneous field-of-view at wavelength λ comprising N endmembers k each with reflectance r k,λ and fractional coverage f k is The residual error at wavelength λ is ε λ . Equation (2) can be written as a system of linear equations, e.g., for M bands where M ≥ N ⎡ By minimizing ε λ simultaneously across all bands, for example, through sum of least squares, and careful endmember selection, modeled spectral mixtures are computed. Additional usual constraints are that 0 ≤ f k ≤ 1 and f k = 1.

A. Decision Trees
Rosenthal and Dozier [12] applied spectral unmixing to the Landsat 5 Thematic Mapper using empirical endmembers for: coarse grain snow, fine grain snow, shaded snow, rock, and vegetation. From the modeled f sca , they used classification and regression trees to mask clouds and map the f sca as a faster alternative than iteratively solving (3) for every pixel in a scene.

B. MODSCAG and MODDRFS
In MODSCAG, the MODIS Snow-Covered Area and Grain size algorithm [13], the modeled snow reflectance, i.e., r k,λ for the snow endmember in (2), is computed from Mie scattering [22] and the discrete ordinates radiative transfer (DIS-ORT) model [23]. From these known reflectance values, three endmember fractions ( f k ) are estimated by minimizing ε 2 : snow-covered area ( f sca ), rock/soil ( f rock ) or vegetation ( f veg ), and photometric shade ( f shade ). The observed reflectance R λ comprises the MOD09GA surface reflectance product [24]. The rock/soil and vegetation endmembers are from field and laboratory measurements.
MODSCAG assumes the snow is clean and that the reflectance r snow,λ of the snow endmember varies as a function of the grain radius r g , illumination, and viewing geometry φ is the solar azimuth, μ the cosine of the illumination angle, φ v the sensor view azimuth, μ v the sensor zenith cosine, and ζ atm atmospheric absorption and scattering properties. Generally, directional effects of the reflected radiation can be ignored [25], so set φ − φ v = 0 and μ v = 1. Most snow contains dust or soot in sufficient quantities to degrade its visible albedo [26]- [29]. To account for LAP, MODDRFS-MODIS Dust and Radiative Forcing in Snow [18]-accounts for degradation of albedo by LAP at wavelengths up to 0.876 μm. The MODDRFS product is created by comparing observed dirty snow spectral reflectance and modeled clean snow spectral reflectance, using pixels with nearly 100% snow cover, to estimate broadband albedo based on the integration of differences across visible bands. Rather than using the grain size from MODSCAG, MODDRFS computes a Normalized Difference Grain Size Index (NDGSI) From this grain size, a clean snow reflectance is computed using (4). Unlike in MODSCAG, directional effects for reflected radiation are not ignored in MODDRFS. Using scalar values that are not provided, r (clean) snow,λ and the observed reflectance R λ are converted to spectral albedos α (clean) snow,λ and α (dirty) snow,λ . Finally, the difference between the observed dirty snow reflectance from MODIS and the modeled clean snow, integrated across the visible spectrum, is expressed as vis These quantities can be used to calculate broadband snow albedo as α = α (clean) snow −c vis , where c = 0.63 [19], the fraction of solar energy in the spectral range of (7).

C. MODImLAB
Instead of using a surface reflectance product, the solution from the MODIS Image Laboratory [30] at the University of Otago starts with top-of-atmosphere MODIS Terra Level-1B swath products and applies atmospheric and topographic corrections that account for shadows and multiple reflections at the surface. Using pixels with positive NDSI only, (3) is solved iteratively using four snow endmembers (ice, medium granular, coarse granular, and "transformed," apparently dirty). Two vegetation endmembers (pasture, rain forest + brush) were also used. Because the study focused on glacierized areas of New Zealand, glacier debris and rock endmembers were also used for a total of eight endmembers.
After snow and other endmember fractions are computed, a final cleaning step is applied. Pixels with Euclidean norms of ε λ that exceed the Euclidean norm of R λ are set to zero f sca . A 3 × 3 majority filter is applied to a mask of the snow cover. Then image dilation and intersection with the f sca values remove isolated spurious snow-covered pixels.

III. SPIRES APPROACH
Snow Property Inversion from Remote Sensing (SPIReS) builds on these previous spectral unmixing approaches in snow, optimizes processing for large numbers of pixels, and solves for LAP and grain size (thereby albedo) simultaneously. SPIReS estimates surface properties and fractional cover of the snowpack, adaptable to any multispectral sensor that includes visible through shortwave-infrared bands. Therefore, solutions are obtained through a direct, physically based approach that removes the need for empirical relationships.

A. Surface Reflectance
The standard surface reflectance product MOD09GA Collection 6 for MODIS [24] and the on-demand surface reflectance product Collection 1 for Landsat 8 operational land imager (OLI) [31] are used. The Landsat Analysis Ready Data Surface Reflectance product [32] can also be used. While there are limitations of these surface reflectance products in complex snow-covered terrain, such as not accounting for multiple reflections [30] and problems with aerosol optical thickness [33], we find that the utility of these mature products and their global availability outweigh their limitations.
Snow albedo, a hemispherical measure at the receiving instrument, and the directional estimate of surface reflectance at the receiving instrument, is assumed equivalent. Since neither MODIS nor the Landsat 8 OLI sensors measure the entire hemisphere, this assumption is not correct, but the standard atmospherically corrected reflectance products do not provide the reflectance distribution function. Moreover, for sensors of ∼30-m spatial resolution (e.g., for Landsat 8 OLI), digital elevation models are not of sufficient quality globally to correctly calculate the φ and μ terms in (4). Whereas these terms are accurate on flat terrain, uncertainty in slope and aspect leads to larger errors in mountainous terrain. At even finer spatial resolution such as measured by lidar or photogrammetry [34], the slope and aspect of snow-covered terrain are different from the snow-free surface below.
For sensors of coarser (250-1000 m) resolution, the pixelaveraged values of the φ and μ terms in (4) are better estimated, but the terrain slope and aspect vary across the pixel. Calculating the albedo rather than the hemisphericaldirectional reflectance factor [HDRF,5] introduces less error given the geometric uncertainties in lighting geometry and uncertainties in ζ atm in (4). Sirguey et al. [30] make the same assumption, reasoning that the computational cost of implementing a model that calculates the HDRF is not worth the benefit.

B. Cloud and Other Masks
Because both clouds and snow are bright and cold, cloudsnow discrimination remains an unsolved problem in remote sensing [35]. Because of the differences in spatial and temporal resolution, we apply different cloud masking procedures to MODIS and Landsat 8 OLI imagery. Other masks used include the water mask MOD44W for MODIS [36] and the National Land Cover Database [37] for Landsat 8 OLI, so that water-covered pixels are not processed. We use the National Hydrography Data set [38] to add additional water bodies missed by the standard products. Some dry lakes that falsely showed persistent snow cover were manually masked.
Contrary to approaches that assume only positive NDSI for pixels containing snow [30], pixels with small but detectable snow fractions can have negative NDSI values [39]. Thus, all pixels containing N DSI ≥ −0.5 are processed, while those with N DSI < −0.5 are skipped.
1) Landsat 8 OLI: Whereas some cloud cover issues in MODIS can be addressed using temporal interpolation (Section III-N), long periods between repeat acquisitions (16 days) for Landsat 8 OLI make temporal interpolation less useful. Therefore, users must supply a cloud mask for Landsat 8 OLI, or they can select cloud-free imagery.
2) Initial Cloud Mask for MODIS: For MODIS, an initial cloud mask C (initial) MODIS is applied using the MOD09GA supplied mask and an additional criterion MOD09GA is the MOD09GA cloud mask, obtained from the MODIS quality flags for that pixel's reflectance value; ∧ is the logical AND operator; R 6 is the reflectance from MODIS band 6 (1.628-1.652 μm); and t SWIR is an SWIR threshold, set to 0.2. Because no cloud mask accurately discriminates clouds from snow in MODIS based on reflectance alone [35], this approach was developed using trial-and-error experimenting with different values.

C. Modeled Snow Endmembers
The modeled reflectance of the snow endmember assumes: 1) Scattering properties of individual snow grains can be estimated from Mie theory with the properties of a nonspherical scatterer mimicked by the equivalent radius of the sphere with the same specific surface area [40]. 2) The radiative transfer equation applies, whereby the center-to-center spacing between the snow grains is large enough, compared to the wavelength, that near-field effects can be ignored, as Warren [41] discusses extensively.
3) The LAP can be treated as independent scatterers, recognizing that absorbers within a snow grain will cause more absorption than absorbers between grains [42]. Limitations are: 4) For the same mass concentration, LAP of a smaller size causes more absorption; therefore; the objective is to mimic the degradation of absorbers on snow albedo rather than to estimate the actual concentration, size, or even species of the absorbers. 5) When absorbers are within the instantaneous field-of-view, ascertaining whether they are in the snow, next to the snow, or beneath a shallow snowpack is difficult with multispectral imagers such as MODIS or Landsat 8 OLI.
With these constraints, the scattering properties of individual snow grains and light-absorbing particles are computed using Mie theory. Snow grains are large compared to the wavelengths of shortwave radiation (i.e., the Mie parameter, the ratio of the circumference to the wavelength, is large), so the complex angular momentum (CAM) approximation [43] provides efficiently computed results that smooth out some wiggles in Mie calculations that otherwise require integration over a size distribution of snow grains. For small Mie parameters, less than ∼20, the CAM approximations loses accuracy, so for the LAP, which are smaller, Wiscombe's [22] method is used. In either case, the needed wavelength-dependent variables comprise the scattering Q sca and extinction Q ext efficiencies, the single scattering albedo = Q sca /Q ext , and the Mie asymmetry parameter g. The calculations use the complex refractive indices of ice [44] and the contaminating absorber, dust [45] (for dust from the Colorado Plateau), or black carbon [46].
To calculate the scattering properties of the snow-LAP mixture, the mass fractions are treated as dimensionless numbers that sum to 1, and the volume fractions V j are the mass fractions divided by the densities. For each component, the geometric cross section is G j = 3V j /r j , and the weighted averages are To compute the spectral albedo from a snowpack, whereby the reflected radiance is integrated over the upwelling hemisphere, the two-stream approximations are appropriate [47]. Especially for forward-scattering, optically thick media like a snowpack, the delta-Eddington variation provides results that closely match more computationally intensive methods [48].
The scattering properties are modified as follows: For a deep snowpack, "semiinfinite" in the radiative transfer sense, the equations for direct and diffuse albedo, where μ is the cosine of the local illumination angle, are [49] Similar analytic equations are available for shallow snow [49]. The optical depth τ 0 of snow is W is the snow water equivalent (kg m −2 ), r g the grain radius (m), and ρ ice the density of ice (kg m −3 ).
To model the snow reflectance over band passes, the spectral reflectance is convolved with the spectrum of the solar direct and diffuse irradiance. For a band b extending from λ 1 to λ 2 with total spectral irradiance S λ and diffuse fraction q λ (both varying with solar zenith angle cos −1 μ 0 ), the reflectance at local illumination cosine μ is For the spectral distributions of S λ and q λ , we apply the SMARTS model [50] with a mid-latitude winter atmosphere at 3000-m elevation and a rural aerosol specification. For narrower band passes, such as those in the MOD09GA or Landsat 8 OLI surface reflectance products, a simpler approach can be taken using the spectral reflectance of the central wavelength for each band.

D. Mixture Model
A linear mixture model (2) is used. The RMSE ε λ is then minimized using a nonlinear optimization subject to the constraints that 0 ≤ f k ≤ 1 and k f k = 1 to estimate endmember fractions f k and other variables. The attained solution uses the Sequential Quadratic Programming Method [51]. The following four variables are solved: f sca , f shade , r g , and δ, where δ is the dust concentration in ppmw. Limits for f sca and f shade are set at 0 to 1, while limits for the grain radius r g (30 and 1200 μm) are set slightly below and slightly above realistic values such that solutions at the upper or lower limits can be identified and removed at a later step. Limits for the dust δ concentrations are set to 0 and 1000 ppmw.
For both MODIS and Landsat 8 OLI, seven bands from the surface reflectance products are used; thus, there are nine equations for each pixel, seven from (3) plus the two constraints on the fractions. This creates an overdetermined system with five more equations than unknown variables. In the presence of nonnegligible residuals, usually the case with scientific measurements, increasing the difference between the number equations and unknowns effectively minimizes residuals [52].
Linear mixing assumes single interactions between incoming light and a target, therefore, invalid under more complex interactions involving multiple reflections and/or transmission, for example in tree canopies [11]. For mapping snow, Rosenthal and Dozier [12] suggest that anisotropic reflectance and topographic illumination might comprise additional potential sources of nonlinear mixing, but Painter et al. [13] asserted that nonlinear effects in canopies can be treated in linear models using canopy-level endmembers up to canopy densities where no direct sunlight reaches the snow beneath. Quintano et al. [53] wrote that nonlinear mixing models have not been widely used in remote sensing because they are difficult to implement because of missing information and mathematical complexity, whereas linear mixing models have been successfully used for a wide range of applications, including those with nonlinear effects. To the best of our knowledge, all previous subpixel snow mapping studies [12], [14], [30], [54] use the linear assumption.
For the number of endmembers used, Occam's razor advises parsimony in endmember selection with 4 or fewer suggested [55]. Thus, two initial endmembers are used, with reflectances R (snow-covered) and snow-free R 0 (snow-free). These endmembers differ from traditional endmembers in that they are not pure substances, but rather mixtures. However, since an overarching goal of our approach is to map snow, not other substances, we find the approach justified. Furthermore, the use of a snow-free background measurement from the sensor as an endmember should account for some nonlinear effects, such as canopy transmission [13].

E. Shade Model
A third endmember, f shade , is used for shaded areas and shadows [21]. Based on the "ideal black surface" of Smith et al. [56], a value of zero reflectance is used in all bands for the shade spectra. Reflectances from shadows have different spectra based on what is being shaded, and some evidence shows that using multiple shade endmembers increases accuracy [57]. Thus, multiple shade endmembers could be a future improvement.

F. Adjustment for Perennial Snow and Ice
In glacierized areas or areas with perennial snow and ice, R 0 cannot be snow-free. If unaccounted for, the f sca estimates will be too low in these areas. Thus, f ice , the fraction of ice or permanent snow in a pixel, is estimated using GLIMS, Global Land Ice Measurements from Space [58]. Glacier outlines are converted to a binary mask with lower spatial resolution than the target sensor (e.g., 100 m for a 463-m MODIS pixel or 10 m for a 30-m Landsat 8 OLI pixel), then converted to a fractional value through upscaling via Gaussian pyramid reduction [59], creating a fractional ice estimate. For areas where the GLIMS database entries are created from Landsat imagery, this approach does not resolve at 10-m resolution, but for many areas of the Western U.S., glacier outlines are available at 10 m [60].

G. Canopy Cover
To account for snow hidden by a forest canopy, it is assumed that snow detected within canopy gaps persists below the canopy cover as well. To estimate the canopy fraction, the viewable gap fraction (VGF) is used, representing the fraction of the surface between the canopy gaps that contributes to the reflectance of the pixel.
The VGF varies with canopy structure, topography, and satellite view angle. To adjust the VGF based on the MODIS viewing geometry, we use the Geometrical Optical model [61]. For the canopy crown shape parameter b/r , the ratio of the average vertical crown radius (b) to the average horizontal crown radius (r ), we use a value of 2.7, based on mean values for Lodgepole Pine measured during the NASA Cold Land Processes Experiment [62]. Static fractional canopy cover estimates f (static) cc are from the MODIS Global Vegetation Yearly product [63]. As the MODIS off-nadir view angle increases, the canopy cover viewed by the sensor increases and the VGF decreases [64]; therefore, for MODIS, a dynamic f cc estimate is computed that changes daily as the viewing geometry changes (Section III-N).
The observed snow-covered area f (obs) sca is adjusted to f sca in a single step with the estimates of shade, perennial snow and ice, and canopy cover. Fig. 2 shows how the values affect the estimate of the fractional snow cover Fig. 3. Degenerate reflectance for patchy snow versus optically thin snow. The second endmember is the Manzanita shrub (Arctostaphylos) from the ECOSTRESS library [83]. In the snow/shrub mix, the f sca = 0.6, but the snow is optically thick, while in the shrub under snow, f sca = 1 but has only 1 mm of water equivalent. Other relevant parameters are: r g = 44 μm and μ = 2/3. MODIS band passes are shown in gray.
For Landsat 8 OLI, the Global Forest Cover Change product [65] is used, and a more aggressive correction is employed. If snow is detected in a pixel with canopy, i.e., f (obs) sca as a minimum detection threshold (Section III-K), then it is set to f sca = 1. This correction yielded nearly unbiased results for Landsat 8 OLI (Section V), but the same correction caused large positive biases with MODIS.

H. Interpolants
Mie scattering calculations and the iterative inversion process (3) are computationally intensive. To speed computation, linear interpolants are used to estimate the Mie scattering parameters for snow, dust, and soot at a range of radii of the scatterer.
Interpolants were created for each MODIS or Landsat OLI band used in the analyses from (13), for a range of grain sizes, illumination angles, and dust or soot concentrations. Some assumptions are necessary to effectively use the interpolants. 1) An optically thick snowpack is assumed. This is often not the case, for example, near the snow line. The spectral albedo of thin snow can be modeled, but the reflectance spectra for patchy thick snow and continuous thin snow are too similar to distinguish in the inversion (Fig. 3). 2) It is assumed that albedo reduction from LAP occurs from dust with a radius of 3 μm with measured complex refractive indices from the Colorado Plateau [45] for the Mie scattering calculations. This approach sometimes produces unrealistic LAP concentrations, and usually both dust and black carbon are present throughout the season at levels high enough to affect the snow albedo [28], [66]. Like shallow snow, the effect of LAP mixtures on the snow can be modeled, but adding an additional unknown to an inversion with a multispectral sensor is unlikely to improve LAP estimates. Moreover, considering just one absorber mimics effects of multiple absorbers on snow albedo, and the objective is to estimate albedo rather than size and concentrations of absorbers. For the multispectral instruments tested (MODIS and Landsat 8 OLI), almost identical albedo estimates are produced when varying the absorber (dust or soot) and particle size, owing to degenerate spectra for many combinations of dust and black carbon in snow (Fig. 4).

I. Clustering
To further reduce computation time, pixels with similar characteristics are grouped. Reflectance values in each scene are rounded to the nearest hundredth and arranged into matrices R (measured reflectance) and R 0 (background or snow-free reflectance), each of size N × M. The M columns cover the bands and μ (the cosine of the illumination angle for R). There are N rows, one for each pixel. A unique row function F u selects the pixels that are unique within a scalar tolerance t u The F u function compares each row to all the other rows and only keeps those that have one or more columns that differ by ±t u . The result matrix U contains many fewer rows, with each row having one or more columns with values that are unique above tolerance t u , which was set at 0.05. The function also returns vector i u , which contains indices for all rows to the quasi-unique rows.
Each row in U is then solved for f sca , f shade , r g , and δ. The remaining rows in the image that were part of R are then filled in using i u . This approach is about 50× faster than iterating over every pixel. The broadband albedo can then be calculated using r g , δ, and μ [19].

J. Limitations When Solving for LAPs and Grain Size
With relatively coarse spectral resolution, multispectral sensors such as MODIS and Landsat 8 OLI cannot distinguish between dark objects in the visible wavelengths that emerge as the snowpack recedes (such as vegetation) and increasing LAP concentrations at the surface. Both phenomena typically occur in the spring as local sources are uncovered [18], [67]. Likewise, analysis of MODIS time series results from SPIReS of mixed pixels sometimes shows a grain size reduction with decreasing f sca , a physically implausible result. The explanation is that in mixed pixels with receding snow, reflectance values in the SWIR bands increase as soil or vegetation is exposed. Thus, dirty snow can sometimes be confused with fine grained snow which has elevated SWIR reflectance compared to coarse-grained snow. Mixed pixels with fine and coarse snow have different spectra, but the MODIS and Landsat 8 OLI sensors lack the spectral accuracy and resolution to distinguish between them, resulting in the implausible receding fine-grained snow. Moreover, MODIS often mistakes clouds for fine-grained snow, which further causes a potential negative grain size bias. A related problem is that grain size estimates are too small when compared with in situ estimates [19] when the shade endmember is used. This artifact results from darkening in the SWIR bands, which grain growth causes, being mistakenly attributed to the shade endmember during unmixing.
To address these problems, 1) Relatively pure pixels ( f shade < 0.01 and f (obs) sca ≥ 0.90) are used to interpolate r g and δ to pixels that are more heavily mixed; 2) LAP concentrations are set to zero for smaller grain sizes (r g < 400 μm) since LAP degradation of snow albedo increases with grain size [68] and occurs more frequently at larger grain sizes [69]; 3) During melt out, the MODIS estimated r g is fixed at a measured peak value; 4) Likewise, the value of δ is fixed at that time, given the difficulty of accurately measuring LAP concentration and r g during rapid ablation.

K. Detection Limits
If no lower threshold is applied, any spectral unmixing approach will map small portions of snow everywhere. Under ideal conditions, SPIReS can detect snow down to 0.01 for MODIS and 0.07 snow cover for Landsat 8 OLI (Section V). However, because of clouds, skewed viewing geometry, and sensor noise, we apply a conservative f (min) sca = 0.10 to the MODIS retrievals and set all pixels below this threshold to f sca = 0. The Landsat 8 OLI f (min) sca was set to 0.07. Grain sizes and dust concentrations are then set to null for values of zero f sca .

L. Elevation Filter
Without an elevation filter, persistently cloudy areas are often mapped as snow, such as the California coast or California's Central Valley during Tule fog. Thus, a minimum elevation filter is used. For example, for the Sierra Nevada, all pixels below 1000-m elevation are set to f sca = 0.

M. Persistence Filter
The daily temporal resolution of MODIS is leveraged to reduce errors through persistence filters and time-space smoothing.
Even after filtering for detection limits and elevation, invariably some pixels that contain clouds are misidentified as snow.
Since clouds tend to persist for shorter periods than snow, a persistence filter is applied to eliminate false positives (FPs) for snow. The persistence filter F p creates a 3-D logical matrix (mask) P with dimensions x, y, and time. Starting with an initial mask P 0 , and operating along the time dimension, the filter modifies the mask using a centered moving window. Pixels that fail to show snow for at least t p total days are set to false across the window. Initial P 0 is set with f sca ≥ f (min) sca ∧ r > r g,min as snow (true) and false otherwise. The minimum grain radius r g,min is set to 40 μm. The threshold t p is set at 13 days of snow cover and the window size is set at 45 days in the forward and reverse time (90 total days). Pixels that fail the persistence filter are assigned f sca = 0. For regions like the Eastern Himalaya that experience frequent snow cover of short duration during the summer monsoon [70], these thresholds likely require adjustment.

N. Smoothing Over Time
Temporal smoothing uses weighted splines for each pixel, with weights based on viewing geometry and other sources. Weights, normalized from 0 to 1, are computed from the following attributes: cloud, cloud shadow, land versus water, aerosol specification at 550 nm, presence of cirrus, atmospheric correction, and adjacency. Viewing geometry is specified by a pixel weight w p = ( p p ⊥ ) −1 that penalizes off-nadir views. The weights depend on the normalized alongtrack p and cross-track p ⊥ dimensions of the pixel [71].
These pixel weights are then used to smooth f sca , r g , and δ in the time dimension. Because MODIS is a scanning instrument, considerable spatial distortion already occurs in the projected MOD09GA values, and since r g and δ are already interpolated spatially (Section III-J), we minimize spatial smearing by smoothing only in the time dimension. f sca is smoothed first with splines. Grain sizes r g and δ outside realistic ranges are set to null values to be interpolated. For r g , the acceptable range was set at 40-1190 μm. Although snow grains can be smaller, grains below 40 μm could be ice clouds instead. The upper 1190-μm limit was selected to be just below the maximum value of 1200 μm, which identifies a poor numerical solution. Similarly, for δ, the range was set at 0-950 ppmw, just below the 1000-ppmw maximum to identify poor solutions. After interpolation, any remaining values above/below the limits were set to the maximum/minimum values.

IV. VALIDATION
Validation focuses on: 1) viewable snow-covered area f (obs) sca ; 2) snow cover f sca (including under canopy and shadowed areas), i.e., from (14), and (3) broadband snow albedo. f (obs) sca for Landsat 8 OLI and MODIS is validated with nine days of fine-spatial resolution multispectral optical data from WorldView-2 and -3. This approach directly compares f (obs) sca from SPIReS with f (obs) sca from WorldView. To validate snow cover even beneath trees, we use snow depths measured from the Airborne Snow Observatory (ASO) at 3-m spacing. f sca from SPIReS using Landsat 8 OLI is compared with 5 ASO flights, the low number owing to the lack of Landsat 8 imagery on the days of ASO acquisition. f sca from SPIReS using MODIS is compared with 71 ASO flights. ASO only collects data in the spring.
Albedo is validated with terrain-corrected measurements from CUES, the CRREL-UCSB Energy Site [72] on Mammoth Mountain, California (CRREL is the U.S. Army Cold Regions Research and Engineering Laboratory, UCSB is the University of California, Santa Barbara).
To validate snow detection, f sca was converted to a binary mask, true when f sca > 0. To validate viewable snow detection, f (obs) sca was converted to a binary mask, with lower values set to zero using a range of f (min) sca . From the mask, identifications of TP (true positive), TN (true negative), FP, and FN (false negative) support calculations of Precision, Recall, and the F statistic All measures vary from 0 to 1. High Precision indicates few FPs (a pixel mistakenly classified as snow-covered). High Recall indicates few FNs (a pixel mistakenly classified as not snow-covered). The F statistic balances Precision and Recall. Optimal f (min) sca values were determined by maximizing the F statistic.
To account for geolocational uncertainty, SPIReS output is coarsened spatially by a factor of 4: 120 × 120 m 2 for Landsat 8 OLI and 2 × 2 km 2 for MODIS. Validation data from WorldView are available at about 0.5-m resolution; ASO data are available at 3-m resolution.

A. Validation of Per-Pixel Viewable Snow Cover
To validate f (obs) sca , panchromatically sharpened cloud-free WorldView-2 and -3 data within 2 days of a Landsat 8 OLI acquisition or near-nadir MODIS acquisition were used. We manually checked imagery to verify no snow fell between the two acquisitions. The spatial resolution of these validation data ranges from 0.34 to 0.55 m, depending on the view angle of WorldView. Table I shows the summary statistics for nine combinations of WorldView, Landsat 8, and MODIS scenes. Validation images from December to June were selected to account for variability in illumination conditions, snow cover, and snow albedo. The imagery spans diverse locations across California's Sierra Nevada that represent heavily forested western slopes, higher elevations, and drier eastern slopes. The WorldView images range from well illuminated alpine scenes above the tree line in June to heavily shadowed scenes below the tree line in December. The snow-covered WorldView pixels are assumed to be pure endmembers of 100% snow, which are then coarsened to Landsat or MODIS spatial resolutions. Neither WorldView, Landsat, nor MODIS can see through thick tree canopies, so the validation compares snow that an optical sensor identifies.
The volume of the WorldView data (e.g., 16 million 0.5-m pixels comprise a 2-km pixel) and the size of the WorldView scenes (on average 4.3 billion pixels) mean validation data cannot efficiently be generated via manual classification of pixels by an analyst. Two prior studies [7], [73] demonstrated that pan-sharpened RGB WorldView data product can be used to validate f sca , but both were limited to mostly unvegetated, sunlit, snow-covered slopes in alpine regions.
Here, multispectral (including near-infrared bands) WorldView data are used, extending the use of this imagery for snow cover validation in more complex scenes with shade and forest cover. These multispectral data are delivered as an 8-band product at ∼2-m resolution alongside a panchromatic band at ∼0.5-m resolution. The 2-m data were pan-sharpened using the Gram-Schmidt mode 2 algorithm with Generalized Laplacian Pyramid decomposition (GS2-GLP) available in the MATLAB pan sharpening toolbox v1.3 [74]. WorldView-2 and -3 bands 2-7 overlap the bandpass of the 0.45-0.80-μm panchromatic band, so only these five bands were sharpened. The WorldView-3 multispectral bands include a red-edge band from 0.706 to 0.746 μm and a near-infrared band from 0.772 to 0.890 μm, aiding in classification of snow cover.
A k-means classification was used with a final manual adjustment to automate most pixel classifications and retain the ability to fine-tune pixel output classes. This approach balanced validation accuracy with data volume to quickly generate a large library of diverse validation scenes.
A small fraction of each WorldView scene was randomly subsampled and run through a cluster evaluation algorithm that performed 5 replicates of k-means with 5-30 clusters. For each replicate and each set of clusters, the Caliński-Harabasz evaluation criterion [75] was calculated. This criterion, the ratio between the within cluster variance and the between cluster variance, is largest at the number of clusters where within-cluster variance is minimized and between-cluster variance is maximized. The maximum of this score during cluster evaluation was used to determine the optimal number of clusters for generating the validation data. Once the optimal number of clusters was determined, k-means was run on a data subset to determine the final centroids for each cluster. These centroids were then used to classify the entire image into individual clusters. To account for the range of illumination in each scene and to accurately cluster all the viewable snow pixels separately from other land cover types, multiple classes for snow cover were used within each scene. To assign these individual snow clusters to the final snow mask, each cluster was classified as snow-covered or snow-free using random subsamples of 5000 × 5000 pixels, with masks of the individual clusters to guide the decisions. This reduced the manual classification burden from billions of decisions for individual pixels within a scene to a maximum of 30 classification decisions per image. The final validation snow mask generated was the combination of all clusters that mapped snow in the image.

B. Validation of Snow Cover Even Under Trees
The ASO maps snow depth and models snow density to calculate snow water equivalent (SWE) using a comounted lidar and a spectrometer [76]. Estimated snow depth and SWE for specific watersheds across the Western U.S. are publicly available at the National Snow and Ice Data Center. We mapped snow cover using the snow depth product at 3-m resolution [77]. The lidar can penetrate forest canopies, and returns are unaffected by terrain illumination, thus providing the best available map of snow. The ASO 3-m snow depths were treated as snow cover if depth exceeded 0.1 m; then, they were coarsened to the analysis resolution and compared to f sca . ASO limits its acquisitions to mostly high-altitude basins in the spring when water managers plan releases from reservoirs, but snow exists within the images across a range of canopy covers. TN values ( f sca,ASO = 0) were excluded from the RMSE because some scenes contain large snow-free areas and most snow mapping algorithms work well at identifying them [9].

C. Validation of Snow Albedo
Broadband snow albedo measurements were estimated from radiometer measurements at CUES over the 2017-2019 water years. Since the snow surface seen by the down looking radiometer is rarely flat or level in windy areas, the incoming radiometer measurements were corrected to the slope and aspect of the snow surface [78] using hourly scans from a terrestrial laser scanner that operates automatically. We applied other filters to eliminate shadows cast by trees and to remove times with diffuse lighting. Details are provided in [79] with the following changes: 1) only broadband albedos at the local MODIS overpass time (∼10:30 AM) were measured; 2) an adjustable self-leveling boom was used to keep the down looking radiometer about 1 m above the snow surface such that trees and other objects were out of the field of view; 3) measurements when the snow surface was heavily sun cupped (usually beginning in June) were discarded. Previous work [19] shows that these features cause unrealistically dark snow and nullify the plane assumption for the snow surface correction. Also, to account for geolocational uncertainty and spatial variability of the snowpack, the closest match from a 9-pixel neighborhood around CUES was used. Albedo from Landsat 8 OLI was not validated since none of the coincident WorldView or ASO scenes covered CUES.  The WorldView snow masks are used to determine the lower detection limit of f (obs) sca , classifying snow only from the optical properties of an individual pixel and the detection performance with this lower threshold. The f (obs) sca values have no other thresholds or adjustments applied: no minimum elevation, canopy cover, or ice fraction adjustments. To determine the lower detection limit, the F statistic is calculated with a lower detection limit threshold set on both WorldView and MODIS or Landsat 8 OLI data at the sensor resolution. Pixels with snow-covered area below the threshold in either data set were considered free of snow, and any at or above the threshold were considered snow-covered. Fig. 7 shows the F statistic across a range of possible limits for the lower detection limit of SPIReS from 0.01 to 0.20 snow-covered area for Landsat 8 OLI and MODIS. Table I summarizes the binary statistics from the WorldView scene validation for both Landsat 8 and MODIS at their corresponding lower detection limits determined in Fig. 7.   The F statistic is maximized at 0.858 for a detection threshold of 0.01 snow-covered area for the MODIS sensor at 463-m resolution. The F statistic is maximized at 0.900 for a detection threshold of 0.07 snow-covered area for the Landsat 8 OLI sensor at 30-m resolution.
SPIReS f sca values were validated against the ASO-derived f sca with no lower detection limit applied to the ASO data.   Table II summarize the validation for 71 MODIS images that corresponded to an ASO image, and Fig. 8(b) and Table III show validation statistics for the five Landsat 8 OLI scenes that corresponded to an ASO image. Fig. 9 shows the F statistic across binned values of snow-covered area and canopy cover. For both MODIS and Landsat, the F statistic improves slightly with f sca . The F statistic decreases as canopy cover increases because the canopy obfuscates the snow. The RMSE and bias do not show trends related to either snow-covered area or canopy density.
Albedo was validated at CUES for water years 2017-2019 across all 186 days with valid measurements. The mean bias is 0.67% and the mean RMSE across all days and all years is 4.6% (Fig. 10). The remotely sensed albedo tracks with measurements during accumulation and throughout the melt season, though it tends to underestimate the higher values earlier in the season and overestimate the lower values later in the season. Some of these differences owe from the coarser pixel scale of the SPIReS estimate versus the plot scale of the CUES data. VI. DISCUSSION AND CONCLUSION For mapping snow, SPIReS performs favorably using a physics-based approach that generates the "right answer for right reason" [80]. In other approaches, the quality of the validation data and error metrics used vary, making intercomparing difficult. Seemingly subtle differences can greatly affect results, such as including areas with no snow in RMSE calculation. Here, extensive fine-resolution snow cover validation data have been employed, derived from multiple platforms. As with all approaches using passive optical measurements, the ability to detect snow degrades as canopy cover increases. In forests and in other areas of low illumination, adjustments are made that assume similar snow cover to areas that can be observed.
Likewise, the broadband albedo errors compare closely with those from MODSCAG/MODDRFS in an extensive validation study that uses terrain-corrected albedos that are more accurate than those from empirical age-based formulae [19].
New in SPIRES are: 1) the ability to measure concentration of LAPs; 2) use of only snow and snow-free endmembers; 3) integrated cloud masking and smoothing; and 4) interpolants and clustering to reduce computation time. Since the source code is available in a version-controlled open repository, the improvements allow users to process imagery on demand.
We will expand the use of SPIReS to other multispectral and hyperspectral sensors. Future improvements include treatment of shadows, sophisticated methods for selecting snow-free endmembers, and sensor fusion.

DATA STATEMENT
The source code for v1.0 release of SPIRES is at https://github.com/edwardbair/SPIRES/releases/tag/v1.0. The spacetime cubes of fractional snow cover, grain size, and dust concentration covering the Sierra Nevada for water years 2001-2019 (October through September) are in the IEEE Dataport [81], along with an expanded Excel file that shows statistics for all scenes in Tables I and II. WorldView-2/3 scenes and corresponding MODIS and Landsat images are in Zenodo [82]. Snow depths at 3-m spatial resolution from the ASO images are available from the National Snow and Ice Data Center [77]. Landsat 8 OLI and MODIS data are available from several sources, for example, the U.S. Geological Survey Earth Explorer (https://earthexplorer.usgs.gov/).