Mapping snow depth and volume at the alpine watershed scale from aerial imagery using Structure from Motion

Time series mapping of snow volume in the mountains at global scales and at resolutions needed for water resource management is an unsolved challenge to date. Snow depth mapping by differencing surface elevations from airborne lidar is a mature measurement approach filling the observation gap operationally in a few regions, primarily in mountain headwaters in the Western United States. The same concept for snow depth retrieval from stereo- or multi-view photogrammetry has been demonstrated, but these previous studies had limited ability to determine the uncertainties of photogrammetric snow depth at the basin scale. For example, assessments used non-coincident or discrete points for reference, masked out vegetation, or compared a subset of the fully snow-covered study domain. Here, using a unique data set with simultaneously collected airborne data, we compare snow depth mapped from multi-view Structure from Motion photogrammetry to that mapped by lidar at multiple resolutions over an entire mountain basin (300 km2). After excluding reconstruction errors (negative depths), SfM had lower snow-covered area (∼27%) and snow volume (∼16%) compared to lidar. The reconstruction errors were primarily in areas with vegetation, shallow snow (< 1 m), and steep slopes (> 60°C). Across the overlapping snow extent, snow depths compared well to lidar with similar mean values (< 0.03 m difference) and snow volume (± 5%) for output resolutions of 3 m and 50 m, and with a normalized median absolute deviation of 0.19 m. Our results indicate that photogrammetry from aerial images can be applied in the mountains but would perform best for deeper snowpacks above tree line.


Introduction
Snow depth and snow water equivalent are essential monitored quantities used in many water resource applications. In mountain environments, snow depth is traditionally measured continuously at instrumented sites (e.g., Snow Telemetry in the United States) or periodically along transects due to the complexity of the terrain. These longterm records are valuable but tend to be at locations that are accessible, lower in elevation, and hold snow for longer than the surrounding terrain (Molotch and Bales, 2005). This limited spatial coverage leaves a poor understanding of snow depth distributions particularly at high elevations or remote and inaccessible locations. To address the gap in spatial coverage and map snow depths across large areas, remotely sensed surface elevation models have been used to retrieve depths by differencing the surface elevation when snow is present ("snow on") and absent ("snow free") (Deems et al., 2013).
The surface elevation differencing principle to retrieve snow depth has been demonstrated from several remote sensing platforms, spanning a range of spatial resolutions, areal coverage, repeat intervals, and varying accuracy. Examples of past used platforms include the Ice, Cloud, and Land Elevation Satellite (ICESat), a laser altimeter that could map snow depths along swaths with sub-meter accuracy (Treichler & Kääb, 2017). The data, however, were of limited utility for mountain regions due to the low temporal resolution and large ground footprint (70 m). Current satellite stereo-photogrammetry derived digital surface models (DSMs), such as those from World view or Pléiades, have the potential for higher spatial (< 1 m) and temporal resolutions (Marti et al., 2016;McGrath et al., 2019;Shaw et al., 2020a;Deschamps-Berger et al., 2020). Limitations, though, include terrain-dependent accuracy with best values of 0.2 m-0.4 m over flat surfaces (e.g., bedrock, ice) with shallow slopes (< 10°C; Shean et al., 2016). Additionally, satellites require a clear view of the target area at overpass time with no obstructions, such as clouds (Shaw et al., 2020b).
With no current space-borne platform providing a consistent combination of high temporal and spatial resolution required for accurate spatial snow depth mapping at watershed scales in the mountains, airborne campaigns have been established to address these needs. For example, the Airborne Snow Observatories, Inc., (ASO), a combined Light Detection and Ranging (lidar) and imaging spectrometer platform, delivers time-series of snow depth maps at 3 m resolution with centimeter accuracy. ASO is currently operating in selected watersheds primarily in the California Sierra Nevada and Colorado Rocky Mountains (Painter et al., 2016). Expanding operations, similar to ASO, to more areas is limited by equipment costs, airplane logistics, and operational expenses. Similar to satellite remote sensing, airborne platforms also have limitations due to weather, but a strength is the option for on-demand surveys once conditions improve.
Although smaller in spatial extent, it has been shown that DSMs from multi-view Structure from Motion (SfM) photogrammetry can map snow depth up to alpine catchments size (Bühler et al., 2016;Harder et al., 2016;Schirmer & Pomeroy, 2020;McGrath et al., 2022). The reported results using a Remotely Piloted Aircraft System (RPAS) mapped snow depth at sub-decimeter resolution with centimeter accuracy. A modern photogrammetry technique, SfM uses multiple overlapping images to re-construct the area of interest as a three-dimensional surface model. Further descriptions of SfM and its underlying algorithms can be found in the existing literature (Westoby et al., 2012;Bühler et al., 2016;Harder et al., 2016). Other applications of SfM include time-lapse photogrammetry using permanently installed cameras (Filhol et al., 2019). RPAS-based SfM studies have gained popularity due to the affordability of high-resolution consumer-grade cameras, ease of SfM application, and ability to create accurate data sets on demand (Gaffey & Bhardwaj, 2020). Compared to airplane surveys, RPAS have smaller areal coverage due to battery life, cannot operate safely outside the line of sight, and face strong regulations in populated or protected areas. These limitations lead to significantly smaller footprints and less consistent coverage relative to large-scale platforms.
Previous studies that retrieved snow depth using output DSMs from SfM and images taken from an airplane (Bühler et al., 2015;Nolan et al., 2015;Eberhard et al., 2021) demonstrated the utility of the method over larger areas, but gaps remain for understanding and efficiently applying SfM at the alpine watershed scale. For example, Nolan et al. (2015) had small study areas (10 km 2 and 35 km 2 ) over relatively flat tundra terrain in Alaska. Bühler et al. (2015) and Eberhard et al. (2021) had larger areas with representative alpine terrain, but the validation data used for accuracy assessment was only available for a small subset of the area. Additionally, both studies did not record the reference validation data on the same day as the images used with SfM. This precludes the ability for a methodological evaluation under identical snowpack conditions, as snow metamorphism is a continuous process, and the depth is unlikely to remain constant. From an image processing perspective, Bühler et al. (2015) excluded all areas with visible vegetation in the snow scenes, inhibiting an assessment of surface types that are challenging to reconstruct with photogrammetric methods. The DSMs were also created in chunks that created an unaccounted-for error, merging multiple DSMs into one output product. Additionally, the images were recorded with a multispectral line scanner, making the results challenging to compare to other SfM studies using RGB camera images. In contrast, Eberhard et al. (2021) had imagery from an RGB camera and reconstructed the study in one step. However, the image processing required manual ground control points (GCPs) placement for the scene to be geo-referenced accurately. Using GCPs makes it challenging to follow this approach in vast, remote, snow-covered areas. Their approach also aligned the Frontiers in Earth Science frontiersin.org 02 snow depth via cubic-convolution resampling, which left no understanding of the individual scenes' geo-location accuracy (snow-free and snow-on).
This work assesses the ability of SfM to retrieve snow depth distributions over an alpine watershed basin with high-resolution RGB imagery captured by an airplane and evaluates the accuracy relative to a simultaneously collected lidar-based retrieval. The comparison is spatially complete, across different output product resolutions, with identical snowpack conditions, sensor viewing geometry, environmental influences (weather, sun angle), and the same data processing for snow-free and snow-on scenes. The presented workflow is open source and could be readily adopted for any airborne collected geo-referenced image data set. Meyer & Skiles (2019) showed that SfM can accurately reconstruct snow surface elevations from airplane imagery over bright, freshly fallen snow, with a relative accuracy of 0.17 m at 1 m resolution for a 3.2 km 2 area, relative to lidar. This work extends the SfM assessment from snow surface elevations to snow depth, while also scaling up to a larger alpine watershed basin (300 km 2 ) with more variety in topography and vegetation. The broader application of SfM will expand our understanding of the strengths and weaknesses of photogrammetric-based techniques and provide an additional option for large-scale aerial snow observations.

Study area
The East River Watershed (hereafter ERW) basin (330,050 E,43,14,650 N;EPSG:32,613), located northeast of Crested Butte, Resort. The East River is one of two primary tributaries of the Gunnison River, which itself discharges into the Colorado River. The ERW is estimated 300 km 2 in size, has an average elevation of 3,266 m, and a vertical relief of 1,420 m (Hubbard et al., 2018). The vegetation varies across the elevation ranges and includes brush and grass land, aspen and mixed conifer, and alpine meadows. The East River was designated as a Scientific Focus Area in 2016, supported by the US-DOE Biological and Environmental Research Subsurface Biogeochemistry Program. The ASO flights, and subsequent data processing, were funded by the state of Colorado to map snow distribution patterns and support water supply forecast improvements.

Images
The ASO images were taken on the 24 May 2018 for the snow-on and 12 September 2018 for the snow-free scene. Flight patterns were almost identical on both days, covering the area with a 50% overlap lawn-mower pattern. Flight altitude varied slightly between the two flights, where the May flight was 6,400 m above sea level, and the September flight was at 6,100 m. There was also a difference in flight line orientation between the two dates, with the May flights in a North-South direction and the September flights in a Northwest-Southeast direction. The orientation for the flight lines during the snow season was selected based on lighting conditions for the ASO imaging spectrometer and flight efficiency, and no direct considerations for the camera.
The camera used by ASO is mounted inside the lidar instrument, creating identical view perspectives between the lidar scanner and the camera to the ground surface. Each image has dimensions of 10,328 × 7,760 pixels with a 16-bit color depth and size of 5.2 microns for an individual pixel. Underlying hardware consists of a medium format Phase One iXU 180-R CCD sensor camera with a Rodenstock 50 mm HR Digaron-W wide-angle view lens. The recording interval for the camera was twelve seconds for the snow-free flight, which resulted in 287 images, and six seconds for the snow-on flight resulting in 582 images. The flight operators selected different acquisition intervals to visually identify conditions beneath the plane during the flight, without considering SfM processing, as ASO does not use the images in its operational workflow. The average ground sample distance (GSD) was 0.31 m/pixel for the snow-on and 0.28 m/pixel for the snow-free images. An overview for both collections is shown in Table 1.

Comparison, reference, and classification data
The quality assessment for SfM used the snow depth (SD; Painter, 2018a;Painter, 2018b), snow water equivalent (SWE; Painter, 2018c), and digital terrain model (DTM, Painter and Bormann, 2020) products by ASO. All of them are published and distributed through the National Snow and Ice Data Center Distributed Active Archive Center and are referred to as ASO SWE, ASO SD, or ASO DTM from here on out. ASO uses the same DSM difference principle to calculate SD, where snow-on values are subtracted from the snow-free. The ASO products used in this study were the SD maps at 3 m and 50 m resolution, SWE at 50 m resolution, and the DTM at 3 m resolution. Other data sets directly acquired from ASO were a land surface classification raster at 3 m resolution and the lidar point cloud from the snow-on acquisition. The latter was only used for the co-registration step, which is described in section 4.2.
The land surface classification map is an output data set from the ASO imaging spectrometer spectral reflectance processing pipeline, categorizing each pixel into the basic land surface types; snow, rock, vegetation, and water. The spatial grid and resolution matched the 3 m ASO SD map and required no further postprocessing. We note that the imaging spectrometer mostly misclassified the water class; it was shading within vegetation, which we confirmed with a subset of all water classified pixels. Therefore, we treated both categories (water and vegetated areas) as one category. Additionally, the spectrometer data is not processed to map fractional land cover, and pixels classified as rock or vegetation indicate that this land surface type dominated the spectral signature within that pixel. Consequently, pixels classified as rock or vegetation can still contain a measured snow depth even if another surface type is dominant. For example, the lidar measured snow between the canopy, but that same snow is only a small fraction compared to the vegetation within one pixel for the spectrometer, and hence the land surface classification would be vegetation. More technical details on the ASO platform and the output product's processing steps can be found in Painter et al. (2016).

Image processing
The camera images from the ASO flights were processed using Agisoft Metashape (version 1.6.2) along with associated geo-location and orientation data from the airplane global navigation satellite system and inertial measurement unit. Metashape was used for feature matching, image alignment, dense point cloud creation, and geo-referencing the point cloud. We refer readers for more technical details on data preparation and settings for Metashape to the workflow in Meyer & Skiles (2019).

Co-registration
After the SfM snow-free and snow-on point clouds were generated, a co-registration step ensured the closest alignment of the surface models to a reference data set. This process minimizes relative geo-location error and improves accuracy for DSM difference products. The co-registration was performed using the Ames Stereo Pipeline (ASP; version 2.6.3), which internally uses the iterative closest point algorithm to determine the difference between a reference and movable point cloud (Shean et al., 2016;Beyer et al., 2018). The reference in this study was the ASO snow-on point cloud and reduced to control surfaces to compute the differences. These surfaces, such as exposed bedrock or roads, have consistent elevation across time and were identified from the ASO classification map. Further refinement to the control surfaces included removing any areas with snow in the ASO SD product and any slopes steeper than 50°C (Shaw et al., 2020a). The bounding box for the reference was extended beyond the ERW basin FIGURE 2 ASO imaging spectrometer classification of the study area. Most snow surfaces were in the higher elevations in the northern locations (A). Control surface (red) distribution over the extended ERW basin boundaries used for co-registration (B).

Frontiers in Earth Science
frontiersin.org 05 boundaries to increase the available area for co-registration ( Figure 2B). Co-registering the SfM point clouds to the ASO lidar point cloud also had the advantage of implicitly aligning the SfM with the ASO SD product.
The co-registered SfM point clouds were converted to a gridded raster product (GeoTIFF) with 1 m resolution using the Point Data Abstraction Library (PDAL; Butler et al., 2021), which provides the inverse distance weighting (IDW) algorithm for interpolation. The IDW algorithm is suitable with a point density of multiple points per square meter (Guo et al., 2010), and both SfM point clouds had sufficient density for its application (23.2 points/m 2 for the snow-on and 31.5 points/m 2 for the snow-free). With the conversion to a raster product, PDAL clipped the outputs to identical bounding boxes and transformed them to matching projection (WGS 84/UTM zone 13N; EPSG 32613). The SfM snow-on and snow-free raster products were downscaled (from 1 m to 3 m and 1 m-50 m) to match the ASO SD map resolutions using the bilinear interpolation option from the Geospatial Data Abstraction Library (GDAL). The final step was calculating the SfM SD by taking the pixel-wise difference in surface elevation between the snow-on and snow-free DSMs.

Comparison
The SD values from SfM were compared to the ASO SD product by treating ASO as the reference since SD mapping with lidar is the more established method. The SfM SDs were compared to ASO SDs using the entire study domain's mean, median, and standard deviation. This comparison included the snow-covered area (SCA), snow volume, and SWE to show how the SD differences propagate to downstream products.
The SfM SCA was calculated as a percentage of pixels containing snow with a positive depth value and relative to the total number of pixels with depths in the ASO SD product. To enable a SWE comparison, the mean basin snow density was estimated from the ASO SWE and ASO SD product at the 50 m resolution. This basin-average density was then used with the SfM SD to calculate the SfM SWE. The SWE, most relevant for water resource forecasting, was then aggregated as a basin sum of all snow pixels and reported as a volume (m 3 w.e.; Fierz et al., 2009). We want to note that this is not the process ASO is using for the official NSIDC SWE product, and details of which are described in Painter et al. (2016).
The terrain factor analysis used the 3 m resolution and first binned the depths by elevation bands to assess similarities in the vertical relief. The elevation data was extracted from the ASO DTM, and the binning width was set to 10 m. Next, a relative pixel-by-pixel difference, subtracting the ASO SD from SfM SD, compared the mean, median, standard deviation, and normalized median absolute deviation (NMAD; Höhle and Höhle, 2009). Median and NMAD for control surface differences, using coinciding areas between the SfM snow-free and snow-on DSM, determined the relative error of the two models and a measure of uncertainty for the calculated SfM SD values. Finally, the snow volume was compared basin-wide and grouped by different surface classifications (snow, rock, vegetation; Figure 2A).
Areas with negative or no SfM SD values, were classified as SfM reconstruction errors, co-registration errors, or both. Pixels with no SD were considered gaps and not further analyzed, whereas the negative SD were inspected by three terrain characteristics: elevation, slope, and aspect. This analysis was done for the full ERW basin and after classifying the negative SD by surface classification type. Aspect and slope were calculated from ASO DTM to create independence from the modeled values by SfM.

Co-registration
Control surfaces, used for co-registration of the SfM snowfree and snow-on scene to the lidar reference point cloud, encompassed 14% of the ERW basin boundaries when gridded at the 3 m resolution ( Figure 2B). The snow-free point cloud was shifted 0.02 m to the North, −0.20 m to the East and −0.41 m in the vertical direction, while the snow-on was 0.01 m to the North, −0.02 m to the East and 0.01 m in the vertical. After applying the translation to the respective SfM point clouds, the control surfaces in the SfM DSMs had a mean difference of 0.01 m with a standard deviation of 0.51 m, and median of 0.03 m. The NMAD of 0.19 m indicated that there were still some areas in the control surfaces that were not constant between the scenes, despite all the refinements to constrain those. With median and mean close to zero, however, the co-registration can be considered successful for the two scenes.

Structure from Motion snow products
Overall, where SfM and ASO products had pixels with measured SD, there was a good agreement ( Figure 3). Notably, deeper snow (> 1 m) and across higher elevations (> 3,500 m) showed similar pixel values for depth. Most pixel differences with negative SfM SD and positive ASO SD values were in areas with vegetation or shallow snow depth (< 1 m). The following sections explain each compared output product (SD, Snow Volume, SWE, and SCA) in more detail and how the two output resolutions (3 m and 50 m) impacted the magnitude in differences.  Table 2 and visual comparison for a sub-area is shown in Figure 4.
As a whole, SfM showed an underestimation of SD to ASO when compared on a pixel-by-pixel basis with SD pixels > 0 m at the 3 m resolution ( Figure 5A). The depth distribution was mainly in the 0 m-5 m range ( Figure 5B), and higher values were more dispersed and highly localized. Other studies that have used ASO SD maps also observed extreme outliers and considered these as spurious snow depth values (5 m in McGrath et al., 2019;6 m in Brandt et al., 2020). We did not remove any high outliers from both data sources for this study and included them in all comparisons.  The distribution of SD values across 10 m elevation bands showed higher accumulation in the upper elevations, but not necessarily at the highest in both data sets (Figure 6). SfM and ASO had increasing depths between~3,200 m and~3,800 m and decreasing values above this range. However, SfM had a higher spread in the elevations between 3,200 m and 3,500 m ( Figure 6A), with 74% of this elevation band classified as vegetation. Overall, the agreement between the two distributions showed that SfM is suited to map snow depth patterns across a range of elevations in complex terrain.
Categorized by the surface classification map and SD > 0 m for both SfM and ASO at the 3 m resolution had most SfM pixels as snow (81%), followed by rock (12%), and the vegetation (7%). For ASO, the pixels were distributed differently across land surface types; 69% in snow, 15% in rock, and 16% in vegetation. The

FIGURE 4
Comparison of snow depth between 3 m, and 50 m resolution. SfM did not accurately capture snow depth within vegetated areas but showed good agreement for open areas at higher output resolution. The higher 3 m resolution also retrieved more areas with snow depth around vegetation and had closer values of SCA compared to 50 m.

FIGURE 5
SfM snow depth values plotted against ASO snow depth values at the 3 m resolution in the coinciding area (A). The dashed line shows a hypothetical one-to-one relationship. SfM tended to underestimate the snow depth compared to ASO at the same pixel locations. The snow depth histograms showed a strong agreement between the ASO (blue) and SfM (beige), in the overlapping areas with measured depth (B).

Frontiers in Earth Science
frontiersin.org 08 higher ASO value in the vegetation class denotes the strength of lidar SD mapping and is further discussed in section 6.3. This comparison was not done at the 50 m resolution as the classification map was unavailable. A simple interpolation and down-scaling of the 3 m product was not considered a reliable classification base map.

Snow water equivalent and snow volume
The estimated amount of SWE was calculated by applying a constant snow density of 394 kg/m 3 , the basin-average value calculated from the 50 m ASO SWE and ASO SD product. This is a simple treatment of snow density, but in the spring snow density varies less than snow depth. With this single value applied to the 3 m ASO SD product, the calculation of basin SWE matched the 50 m ASO SWE product to within 3%. Based on this density estimate, the 3 m SfM product had a basin total SWE of 8.13 m 3 × 10 6 m 3 resulting in a difference of 1.53 m 3 × 10 6 m 3 relative to ASO. At the 50 m resolution, SfM SWE was 8.36 m 3 × 10 6 m 3 with a difference of 1.30 m 3 × 10 6 m 3 relative to ASO.

Structure from Motion snow product errors
The SCA from SfM was 27% less area than the ASO SCA at the 3 m resolution, of which 0.2% was a data gap (no pixel value) in SfM. The missed ASO SD had a mean of 0.47 m, a median of 0.39 m, and a standard deviation of 0.41 m. A pixel-by-pixel comparison showed the majority of the missing ASO area with shallow depths (< 1 m) in ASO ( Figure 7A). The largest negative SfM values (−4 m and −28 m) were primarily found in areas classified as vegetation ( Figure 7B), while smaller negative values (up to −1 m) were predominantly classified as snow.
The regression analysis to explain negative SfM pixel values (SD < 0 m) by the different terrain characteristics showed no strong relationships. Out of the investigated influence factors of aspect, elevation, and slope, the most visible trend was detected when values filtered to only snow (no vegetation or rock) were binned by slope angle and the median depth calculated. The median did not exceed −1 m up until around 60°C, then started to decrease sharply ( Figure 8). Terrain with missing SD and above this slope angle was rare in the study domain though, and most locations stayed below this angle. The observation of a decline in accuracy for photogrammetric reconstructions for steeper slopes is common (Bühler et al., 2012;Müller et al., 2014;Shean et al., 2016;Shaw et al., 2020b), and this The primary focus of ASO is the delivery of lidar-based snow depth and SWE maps, and the camera images are currently not used in data product processing (Painter et al., 2016). With the lidar and imaging spectrometer as the primary data streams, there is little consideration for the image overlap, illumination conditions for the camera sensor, or minimum GSD for further use with photogrammetric reconstruction. Given this image acquisition set-up and the presented results, we believe that although the results here are promising, there is room for improvement if flight campaigns were planned to produce snow depth maps with SfM. For instance, consistent image overlap can improve the quality of SfM output products (Bühler et al., 2016;Harder et al., 2016;Meyer and Skiles, 2019). The potential for airplane-based SfM DSM generation for snow-depth mapping has been demonstrated on a smaller scale by Nolan, et al. (2015), with a geolocation accuracy for the SfM DSMs of ± 0.3 m. Our NMAD of 0.19 m using DSM control surfaces over a larger target area denotes

FIGURE 8
The median snow depth difference for negative SfM values in open areas binned by slope angles stayed below -1 m (orange dashed line) until 60°C (grey dotted-dashed lines), before increasing sharply with most of the terrain below this slope angle (grey histogram). This trend was the only strong relationship detected when analyzing the negative SfM values for possible influences due to terrain.

Frontiers in Earth Science
frontiersin.org the scalability of this technical set-up, and is comparable to Eberhard, et al. (2020), which reported a 0.17 m between a reference and reconstructed snow surface. Another indicator was the point density of the two SfM point clouds; there was an average of 23.2 points/m 2 for the snow-on acquisition and 31.5 points/m 2 for the snow-free acquisition, which signifies well-reconstructed surfaces by SfM. This high point density would allow for higherresolution gridded output products than was analyzed here (to match the lidar resolution) and is a suggested path for further research.

Comparison to other platforms
The SfM NMAD from airplane imagery in this study (0.19 m) shows a higher accuracy compared to satellite-based stereo photogrammetric studies, where the NMAD ranges from 0.36 m (Shaw et al., 2020a) over 0.45 m (Marti et al., 2016) andup to 0.69 m (Deschamps-Berger et al., 2020). Reasons for the higher accuracy are primarily technical, as satellite stereo pairs have, for instance, a lower GSD in a single image, making it more difficult to capture great terrain detail for reconstruction. Additionally, DSM generation from satellite images has a different technical design where stereo photogrammetry uses up to three images (tri-stereo) (Shaw et al., 2020b;Deschamps-Berger et al., 2020;Bhushan et al., 2021). SfM can use any number of images and varies based on the image overlap for an area. On the smaller area scale, using RPAS, the accuracy is higher (cm-scale) than what we have achieved here (Bühler et al., 2016;Harder et al., 2016;Avanzi et al., 2018). The improved accuracy can be attributed to the lower flight altitude, resulting in higher image overlap and GSD. Some of the remaining challenges for RPASs include maintaining accuracy while covering larger areas, currently limited by battery life, and higher sensitivity to weather conditions in alpine areas (Bühler et al., 2016).

Application strength and challenges
Given the limitations of current technologies to remotely measure snow depth, we see SfM with images from airplanes filling an important gap between RPAS-SfM and satellite stereophotogrammetry. Across the two output resolutions and the accomplished NMAD here, we believe that SfM compared well against an active measurement instrument like lidar on larger scales. For one, the SfM product had few pixel gaps (0.2%), indicating that a high camera image sampling rate and sufficient flight line overlap can provide a reliable source for reconstruction. The mean and median snow depth differences for both resolutions were ± 0.05 m for the coinciding SCA by both technologies. These statistics were achieved for mostly open spaces, elevations above the tree line, flatter terrain, and a snowpack of more than 1 m. The higher error with shallow snow depth or vegetation was expected and can stem from multiple sources. For one, any elevation difference between the snow-on and snow-off DTM below the NMAD is less likely to be retrieved successfully. Another source is the SfM snowfree scene with shallow vegetation covering the ground surface with grasses or shallow bushes. These ground cover types are poorly reconstructed by SfM (Avanzi et al., 2018;Fernandes et al., 2018), where lidar can map more accurately (Harder et al., 2020). Snow deposition in the winter additionally compacts this vegetation and is a physical process that cannot be accounted for with the differencing principle (Nolan et al., 2015). Steep terrain is another aspect where accuracy for the results degrades, particularly on angles above 60°C. Here, we argue that the accumulation in those areas is low and exemplified by the snow pixel histogram in Figure 8.  (Deschamps-Berger et al., 2020). Additionally, the difference was small for the volume calculation; SfM had a 3% better match to ASO at the 50 m resolution relative to the 3 m resolution (84% at 3 m versus 87% at 50 m).
We argue that the chosen snow depth resolution also comes down to the use case. Inputs for hydrologic models, for instance, currently do not require or have not been assessed against output resolutions of less than 25 m (Behrangi et al., 2018;Hedrick et al., 2020;Pflug and Lundquist, 2020). The presented results for SWE and the small depth and volume difference between the 3 m and 50 m showed that there is no imminent need for a higher-resolution product. Thus, we see a strong use case for SfM-based retrievals for lower resolution applications, for example, validating or bias correcting snow depth in distributed models, as our results showed good agreement across the metrics of SWE, depth, and volume.

Availability of required data
We acknowledge that ASO has a unique setup, with images, lidar, and land surface classifier data recorded simultaneously.

Frontiers in Earth Science
frontiersin.org 11 The combination provided an opportunity to compare SfM to an established snow depth mapping technology. However, we note that simultaneous data recording is not needed to perform this methodology outside the ASO operation domains. For land surface classification, the snow-on and snow-free products can be directly classified using the photogrammetric DSMs and the image RGB information (Shaw et al., 2020a) or nearinfrared spectrum (Deschamps-Berger et al., 2020), where available. Producing the classification with this approach was beyond the scope of this work and warrants an accuracy assessment by itself before continuing to use in downstream products. Using the existing classification, we reduced a potential source of error for this assessment. Once classified, the DSM co-registration can use any externally sourced point or gridded referenced data set. Solutions that have used control surfaces from third-party sources already exist (Shean et al., 2016). For areas with little change to control surfaces (exposed rock surfaces or roadways), the reference DSM can further be from different recording years and does not have to be from the same year of the images (Midgley and Tonkin, 2017). In the end, the presented processing steps can be applied to any airborne collected and geo-referenced images. A simultaneously recorded lidar-based point cloud reference with image spectrometer classification is not required.

Absence of ground control points
Ground control points (GCP) are commonly used for RPAS-based studies to geo-reference the results, strongly influencing accurate geolocation (James et al., 2017). Our process explicitly did not use GCPs, which reduces manual data recording and data processing intervention and increases automation potential. We believe that image geolocation and perspective information (omega, kappa, and phi) combined with co-registration is a reliable substitute for GCPs, while not compromising output quality. Co-registration is common for photogrammetric snow depth products from satellite images (Shean et al., 2016;McGrath et al., 2019;Shaw et al., 2020b;Deschamps-Berger et al., 2020) and is equally applicable to airplane imagery. The SfM software performed very well with the image metadata; the snow-on model was close to the lidar (North: 0.01 m, East: −0.02 m, Up: 0.00 m), while the snow-free model had a higher shift, predominantly to the East (−0.20 m) and Up (0.41 m) and slightly in the North (0.02 m) direction. We hypothesize that exposed vegetation and ground cover in the snow-free scene degraded the geolocation accuracy. With both alignment adjustments low in magnitude, it is feasible to align the two models to each other and compute snow depth and volume in relative geolocation space. Generally, alpine areas with seasonal snow cover have a higher potential for identifiable control surfaces in both scenes for multi-view image processing (e.g., exposed rock along ridgelines), which can be utilized for co-registration. Alternative surface alignment approaches for environmental conditions with few or no overlapping control surfaces, such as ice sheets, have been developed (Howat et al., 2019;Shean et al., 2019).
From an operational perspective, methods that rely on the presence or collection of ground control data limit the ability to scale the application to larger areas or more frequent acquisitions. This limitation is from a logistical point of view, where more resources (people, equipment) would be needed to complete one survey. Additionally, mountain watersheds are often remote and challenging to access, and travel through mountains in winter requires specialized skills and training, limiting the capability to collect ground data. Current operational environments, like ASO, show that GCPs are not essential for delivering consistent results from regular repeated surveys.

Expanding snow science
The ability to "fill-in" missing information between point-based snow depth measurements using snow depth maps from airborne campaigns has improved our understanding of large-scale snow processes. Improvements include enhanced model capabilities to predict snow precipitation (Behrangi et al., 2018), observe snowfall distributions (Brandt et al., 2020), improve snow energy balance models (Brauchli et al., 2017;Hedrick et al., 2020), and quantify snow depth variability (Zheng et al., 2019). Spatially complete and temporally consistent records have also improved estimates of SWE (Margulis et al., 2019) and streamflow (Shaw et al., 2020a) through assimilation. Although SfM does not yet deliver similar accuracies across all terrain characteristics and land cover classes, it can be used to supplement or build upon lidar data sets. For instance, a first survey could use the higher accuracy lidar for the entire watershed basin, with successive observations using the more cost-efficient SfM for open areas with low vegetation on a subset of the domain that provides sufficient spatial snow observations to infer patterns for the larger initial domain (Pflug and Lundquist, 2020). The results of this work support that SfM could be an option for operations like ASO for repeated observations after the initial lidar flight. From a technical setup perspective, it may be feasible to use images or video from space-borne platforms (Bhushan et al., 2021), adding the option of temporally consistent broad-scale coverage and reducing operational requirements. More cross-platform (airborne vs. space-borne) and technology (photogrammetry vs. lidar) comparisons like Eberhard et al.

Frontiers in Earth Science
frontiersin.org (2020) are needed to understand and evaluate our capabilities.

Recommendations and further research
We would like to see photogrammetry, including SfM, applied to larger areas and more frequent image acquisition to improve our understanding of this technology at the alpine watershed scale. This work was the first investigation for one flight during the snowmelt season, which still leaves the need for applying the presented methods to different times, such as accumulation or peak snow depth. Our results showed an almost spatially complete snow depth product with sub-meter accuracy above the tree line and for deeper snowpack at two output resolutions. The 0.31 cm/pixel GSD for the snow-on scene imagery is coarser than the 16 cm/pixel in Meyer and Skiles (2019). This difference between the GSD is attributed to the lower flight altitude (1,555 m above ground level; AGL) compared to the normal ASO level in this flight (~3,000 m AGL), as the camera image recording interval was identical (12 s). Despite the different GSD, the achieved 0.19 m NMAD at 3 m resolution here was only slightly higher compared to the 0.17 m at 1 m resolution in Meyer and Skiles (2019). Our recommendation based on these findings is that a lower AGL is not required. In contrast, we suggest keeping the higher AGL like the flight used in this study, as it results in broader areal coverage for the same flight time and reduces processing time due to the reduced overlap between images. A faster processing time could also be possible by reducing the number of images used during DSM generation and is an avenue for continued research. The next potential options for this data set also include comparing at higher resolutions and applying the method to more ASO flights as available.

Conclusion
This study's motivation was to investigate whether Structure from Motion should be considered an additional remote sensing data source for snow depth monitoring on the alpine watershed scale. It also emphasized automating the data processing, as much as possible, to be scalable with area size and ready for operational use. Compared to simultaneous lidar-based measurements, the results had almost identical statistics for mean and median for depth, and a similar estimation in volume for open areas. This assessment holds true at 3 m and 50 m resolution. As with previous studies, vegetated, steep, or shallow snowpack areas had high reconstruction errors with no measured snow depth. Consequently, these terrain and snow depth characteristics also caused the missed volume, SCA, and SWE by SfM compared to ASO.
The importance of monitoring the mountain snow water reservoir is well recognized, with seasonal snow depth and snow water equivalent both being identified as "targeted observables" in the most recent Earth Science Decadal Survey (National Academies of Sciences, Engineering, and Medicine, 2018). This United States National Academy of Sciences survey guides upcoming scientific missions and goals for earth observations from space. Targeted observables are priority observations that may not yet have mature measurement techniques but could within the next 10 + years. This recognition brings attention to emerging technologies and incubation funding to mature their approaches and application. Surface Topography and Vegetation (STV) is one such incubation effort, which focuses on high-resolution global topography mapping and topography change. Photogrammetry, along with lidar and radar, was specifically targeted as a measurement technology with potential for maturation (Donnellan et al., 2021). Although the focus is on stereo-photogrammetry, for which satellite capability is well established, the potential for spaceborne multi-view photogrammetry is also promising and equally suitable. Comparisons and data as in this work contribute to the STV effort with methods undergoing active development. Ultimately, the goal is for global satellite-based time-series mapping of snow volume in the mountains and at resolutions needed for water resource management.

Data availability statement
Camera images, lidar point cloud files, and imaging spectrometer ground classification were acquired through personal communication with the ASO team. The snow depth map, SWE map, and digital terrain model are publicly available through the links in the references from the National Snow and Ice Data Center Distributed Active Archive Center download portal. Active Archive Center download portal: The software components used to process data are publicly available on https:// github.com/UofU-Cryosphere/snow-aso and analytical code for the output data is available on: https://github.com/jomey/raster_ compare.

Author contributions
JM and SMS conceptualized the overall study, with helpful contributions from JD and DS on parts. JM performed the data processing and analysis. KB provided ASO data and support. SMS provided financial support for the study. JM wrote the first draft of the manuscript, which was then contributed to by all authors.

Conflict of interest
Co-authors JD and KB were members of the NASA ASO team (which produced the data used in this study). JD is a cofounder of Airborne Snow Observatories, Inc., and KB is currently employed by ASO, Inc., formed as a result of the ASO NASA technology transition effort.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.