Applications of Snow-Covered Areas from Unoccupied Aerial Systems (UAS) Visible Imagery: A Demonstration in Southeastern New Hampshire

1. Introduction

The presence of snow-covered areas (SCA) and their accurate delineation are relevant to hydrology [1,2,3,4,5,6], land surface and atmospheric modeling [7,8,9], and ecology [10,11,12]. SCA variations can also regulate the climate [13,14], and their measurement helps to monitor the effects of climate change [15,16].

Satellite remote sensing has provided SCA estimates since the 1960s [17]. Daily snow cover observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) and semi-monthly observations from Landsat instruments have been available for decades [18,19]. Observations from the Visible Imaging Radiometer Suite (VIIRS), the next-generation replacement to MODIS [20], and the Sentinel-2 Multispectral Instrument (MSI) have also been used for SCA delineation [21]. Commercial satellites offer additional near daily imagery [22] to map SCA [23]. More recently, snow cover products have been produced from optical imagery using Unoccupied Aerial Systems (UAS) or Vehicles (UAV), commonly referred to as ‘drones’ [24,25,26,27,28,29].

Satellite SCA observations often have a tradeoff between their spatial and temporal resolution. Sentinel-2 and Landsat have moderate spatial resolutions (10s of meters) with revisit intervals of approximately one week. Moderate spatial resolution can provide information on snow variability at field scales (< 500 m), though cloud obscuration paired with weekly revisit intervals leads to substantial information loss over time. In contrast, MODIS and VIIRS have daily or sub-daily revisit intervals with coarse spatial resolutions (100s of meters to kilometers). Commercial satellites offer near daily imagery at high to moderate spatial resolutions (~5 m) [22], but they are more costly than publicly available satellite imagery. Daily observations are ideal for observing melt evolution, though coarse spatial resolution limits their ability to provide relevant information on snow cover variability at field scales. While the aerial extent of UAS observations is smaller than satellites and crewed aircraft, they readily observe optical snow properties at high (i.e., fine) spatial resolutions (<1 m). UAS are also low cost, offer flexible flight planning, can carry a range of light-weight sensor payloads tailored to the needs of a particular study, and are not limited by cloud cover like satellites [30,31,32,33,34]. However, UAS can be constrained by the availability of trained operators, site accessibility, and weather conditions (e.g., cold temperatures, precipitation, high winds) which can limit or curtail UAS flights [35].

The value of UAS is well documented for mapping snow depth using UAS-based Structure-from-Motion (SfM) photogrammetry or Light Detection and Ranging (LiDAR) [35]. Although UAS platforms are routinely outfitted with electro-optical (EO) cameras that could provide a low-cost method of monitoring field-scale snow cover, surprisingly few studies have used UAS to estimate SCA. Because UAV payloads typically lack shortwave infrared observation data needed to calculate NDSI, alternative approaches to estimate SCA using only visible imagery are emerging [36]. Niedzielski et al. (2018) used UAS visible imagery for snow cover detection with cm-level resolution [26]. Belmonte et al. (2021) used snow cover derived from UAS visible imagery to explore snow persistence and its drivers in forested environments [29].

While coarse satellite records are sufficient for global-scale applications [37], there are cases where finer spatial resolutions are needed [38,39]. High to moderate spatial resolution observations have enabled studies on microclimates [40,41], small basin hydrology [42,43], agricultural management [44], and mobility planning [45]. For these applications, snow cover mapping requires approaches with resolutions that can capture heterogeneous snow cover in areas with mountainous terrain [46,47,48], forests [23,28], and transitional zones between fields and forests [49]. For example, forest edge effects and microtopography can influence snowpack energy balances and drive melt variations over distances of less than 1 m [49,50]. Accurate representation of snow patch distribution is particularly important for estimating snow depth, snowmelt rate, and runoff because snow/bare soil patches cause considerable local variations in surface temperatures [51,52,53,54,55,56,57,58,59,60,61].

UAS imagery has recently been used to downscale MODIS snow cover products in the Tibetan Plateau [24,25]. The spatial variability of snow has also been modeled extensively using regression-based and statistical models [62,63], distributed physics-based numerical models [64], and machine learning models. Since 2020, machine learning methods like random forests (RF) and clustering have been widely applied to model snow [65,66,67,68,69]. Machine learning models have been used to downscale MODIS SCA using Landsat and Sentinel-2 SCA observations by training predictive models over mountainous terrain [48,70]. These models have also identified drivers of fine-scale snowpack variations using UAS snow depth observations [71]. High-resolution UAS snow cover observations can be used as training datasets for machine learning models that downscale satellite observations to sub-meter scales. Furthermore, UAS can observe microtopography and vegetation at equivalent scales and, since these variables are the dominant controls on snowmelt and accumulation patterns, provide useful model inputs [50,72,73,74,75]. Despite the value of static features, such as vegetation type or slope, as time-stable controls of snow distributions across seasons [71,76], capturing a changing snowpack still requires some measure of dynamic conditions [77]. Using high-resolution UAS observations of land-surface features and routine satellite observations of snow cover to train machine learning models can effectively combine static and dynamic variables for monitoring snow cover evolution.

Combining satellite and UAS snow cover observations leverages a satellite’s routine temporal coverage and expansive geographic extent with a UAS platform’s high resolution and limited extent to understand snow process and assess the performance of snow models [78,79]. Here, our goal is to demonstrate the value of UAS observations for mapping SCA at high resolutions, identifying fine-scale controls on snow persistence, and modeling snow cover depletion patterns. First, we apply and evaluate Niedzielski et al.’s (2018) k-means clustering approach [26] to detect snow cover and patterns during ablation (i.e., melt) from repeat-pass UAS visible imagery at fine scales. Next, we compare UAS-derived fSCA to satellite remotely sensed snow cover from Landsat and Sentinel-2’s moderate-scale observations and MODIS and VIIRS’s coarse-scale observations as well as UAS LiDAR derived SCA. Finally, we train a RF model using the UAS snow cover observations, physical surface features mapped using UAS LiDAR, and Sentinel-2 NDSI observations and assess the RF model’s ability to capture fine-scale SCA and snow-free patches. The results are used to quantify the performance and establish the role of UAS in mapping snow-covered areas. This work uses a case study conducted in a relatively homogenous open area in an ephemeral snow region [37,80] in Southeastern New Hampshire, USA, to demonstrate our approach.

2. Materials and Methods

2.1. Site Description

The University of New Hampshire’s Thompson Farm Research Observatory, which served as the field site, is located in southeastern New Hampshire, United States (43.1089° N, 70.9485° W) with elevations ranging from 10 to 35 m (Figure 1). Analyses were performed using UAS observations from the open field portions of the site indicated by the hatched area in Figure 1. This portion of the field site is approximately 7.6 hectares (76,000 m² or ~0.1 km²). The site has limited topographic relief with gentle slopes where low-lying areas within the field and clayey soils result in areas of surface water collection and icing during the winter months [81]. Field vegetation varies between low-cut pasture grass and tall (~1 m) unkempt grasses. The surrounding forest canopy is 10 to 30 m tall. This forest contains mixed deciduous and coniferous trees primarily consisting of white pine (Pinus strobus), northern red oak (Quercus rubra), red maple (Acer rubrum), shagbark hickory (Carya ovata), and white oak (Quercus alba). Vegetation shading provides the primary control on solar exposure, with most of the field experiencing two to five hours per day of shade with up to eight hours per day at the forest-field boundary in the north portion of the field [71]. Pixel coverage areas from this study’s satellite remote sensing platforms, which range from approximately 400 m² (Sentinel-2) to 200,000 m² (MODIS), are indicated in Figure 1.

The Thompson Farm field site is located within the ephemeral snow class, which is characterized by relatively warm temperatures, shallow snow, and a total duration of snow cover lasting for an average of 70 days per year [37,80]. Based on a 20-year (2002–2021) winter climatology at the co-located U.S. Climate Reference Network site (USRCN, site ID: 1040), January is the coldest month with an average temperature of -4.1 °C followed by February (-3.1 °C) and December (-1.2 °C). These are the only months with average temperatures below 0 °C at the site. From December to February, monthly mean high temperatures range from 0.6 °C to 3.4 °C with lows ranging from -9.3 °C to -6.0 °C. The average monthly liquid equivalent precipitation is highest in December (122 mm) and lowest in January (78 mm). February, March, and April have approximately 100 mm of precipitation on average (98 to 108 mm). The site’s combination of sub-freezing temperatures and regular precipitation results in frequent snowfall. Conversely, because above-freezing temperatures are also common, there are intermittent melting and rain-on-snow events throughout the winter. Based on MODIS observations 2002 to 2021, snow cover occurs in November (average monthly fSCA is 0.02), December (0.36), January (0.70), February (0.73), March (0.48), and April (0.04).

For the study period, December 2020 through April 2021, precipitation, air and soil temperatures, and solar radiation were obtained from the USCRN site. Snow depths were obtained daily at approximately 10 a.m. local time using Wingscapes Timelapse Camera time-lapse images of a marked snow stake.

2.2. Data

2.2.1. UAS Observations

A series of Red-Green-Blue (RGB) visible imagery and LiDAR observations were collected by UAS over Thompson Farm during the winter of 2020–2021. From February 4 to April 2, 2021 (snow-off flight), there were 15 UAS RGB and 17 UAS LiDAR data collections. Eleven of those data collections had both RGB and LiDAR on the same day (Figure 2). The flights were scheduled throughout the winter to capture a range of snow conditions with near daily flights during the March snowmelt. The first UAS flights occurred on February 4, 2021, immediately following a small snowfall event. The final (15th) flight, on April 2, 2021, was flown shortly after all snow had disappeared from the study area. An additional snow-off LiDAR flight was collected before the study period on December 14 to observe terrain and vegetation structure.

The UAS RGB imagery was collected using a DJI Phantom 4 RTK (real-time kinematic) outfitted with a 20-megapixel EO camera with an 84° field of view (FOV) at an altitude of 100 m above ground level (AGL). This altitude was maintained using a reference bare earth Digital Terrain Model (DTM) spanning the site to achieve an approximate pixel size of 5 cm. Each survey consisted of approximately 600 images which were imported to Agisoft Metashape Pro to combine collected imagery into orthomosaics. The 3-band RGB imagery was the primary input to the snow cover detection algorithm.

LiDAR data was acquired using the Velodyne, Inc. VLP-16, flown at 8 m/s on a DJI Matrice 600 Pro at an altitude of 65 m AGL. The VLP-16 has vertical and horizontal FOVs of 30° and 360°, respectively. The sensor was mounted with the scan direction perpendicular to the flight direction and provided approximately 300,000 observations per second with an accuracy of ±3 cm at a range of 100 m. Inertial navigation system (INS) information of the sensor payload was collected with an integrated Trimble APX-15. For each point cloud, ground returns were classified using the progressive morphological filter algorithm [82], and snow depth was calculated by subtracting a snow-off DTM from a snow-on DTM as by Jacobs et al. (2021) [81]. The LiDAR snow depth estimates compared well to in-situ observations, with a mean absolute error (MAE) of 2.3 cm, root mean square error (RMSE) of 2.8 cm, and a bias of +0.2 cm [83]. While the LiDAR data were primarily used to characterize field topography, LiDAR snow depth observations within the field were also used to derive a 1 m UAS SCA product (hereafter referred to as UAS-LiDAR-SCA). UAS-LiDAR-SCA is a binary product where pixels having snow depths greater than zero were considered to be snow covered and zero or less were considered to be snow-off (i.e., bare earth).

2.2.2. Satellite Observations

The suite of satellite snow cover products used in our study, MODIS, VIIRS, Landsat, and Sentinel-2, have been extensively validated using other satellite products, aerial photography, field cameras, National Weather Service data, Snow Telemetry (SNOTEL) stations, in-situ observations, and operational model output, and at a range of locations (e.g., [17,18,20,21,84,85,86,87,88,89,90,91,92]). The satellite NDSI products all use a common approach of differencing visible (VIS) and shortwave infrared (SWIR) band reflectance values as follows:

$$NDSI={\displaystyle \frac{VIS-SWIR}{VIS+SWIR}}$$

(1)

All available satellite products from MODIS, VIIRS, Landsat, and Sentinel-2 were extracted for the field site around the UAS collection dates, extending more than a week prior (January 24, 2021) and a few days after (April 6, 2021) to ensure observations were included both before and after field collections (74 total days). Only pixels intersecting the field (Figure 1) were retained, and imagery from all products was filtered to retain only observations collected in clear sky conditions (see Figure 2). Satellite snow cover values were averaged within the field area. Satellite observational overpass dates are detailed in Figure 2.

• MODIS

Daily NDSI observations with a nominal 500 m spatial resolution from the MODIS Terra cloud-free products (MOD10A1F, [93]) were extracted for the field area. The cloud-free products use the latest non-cloud-obstructed daytime observations (10:30 local solar time) from band 4 (green, 545–565 nm) and band 6 (SWIR, 1628–1652 nm) via the NDSI to create a continuous daily record of snow cover [94] in the range 0 (snow-off) to 100 (completely snow covered). There were four MODIS pixels intersecting the field site, each covering mixed forest and open field areas. The NDSI for the field area was calculated using an area-weighted average of the MODIS NDSI pixel values. Of the 74 days considered during the study period, 26 (35%) were cloud-free.

2.

VIIRS

Daily cloud-free NDSI products (VNP10A1F, [95]) with a nominal 375 m spatial resolution from the VIIRS instrument onboard the Suomi-NPP satellite were extracted for the field area. The VIIRS NDSI is calculated using band I1 (red, 600–680 nm) and band I3 (SWIR, 1580–1640 nm) with daily overpasses around 13:30 local solar time. Five VIIRS pixels overlapped the field site and were used to produce an area-weighted average NDSI for the field. Of the 74 days considered during the study period 35 (47%) were cloud-free.

3.

Landsat

Two Landsat products, fSCA and NDSI, were included in the study. The Landsat Collection 2 Level 3 fSCA product [96] directly provides a within-pixel fractional snow cover estimate at a 30 m spatial resolution. A 30 m NDSI product was calculated from the U.S. Analysis Ready Data (ARD) products [97] by differencing the green and SWIR band reflectance values observed by the Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and the Landsat 8 Operational Land Imager (OLI) sensor. Landsat 7 ETM+ bands 2 (green, 520–600 nm) and 5 (SWIR1, 1550–1750 nm) or Landsat 8 OLI, bands 3 (green, 530–590 nm) and 6 (SWIR1, 1570–1650 nm) were used to calculate NDSI. NDSI values, which typically range from -1 to 1, were scaled by 100, and negative values set to 0.

Landsat pixels that overlapped the field site were extracted, then the average fSCA and NDSI for the field area were calculated using their area-weighted average. Both the Landsat 7 ETM+ and Landsat 8 OLI sensors have global repeat observation cycles of 16 days, and the combined products have an 8-day revisit frequency over the study area. Of the ten overpasses during the study period, six were cloud-free.

4.

Sentinel-2

Sentinel-2 pixels that overlapped the field site were extracted, then the average fSCA and NDSI for the field area was calculated using an area-weighted average of the pixel values. NDSI values were calculated from the 20 m spatial resolution Sentinel-2 (2A and 2B) Level 2 processed reflectance products (MSIL2A, [98]). Like Landsat, Sentinel-2 NDSI was calculated via equation 1 using Sentinel-2’s band 3 (green, 543–578 nm) and band 11 (SWIR, 1565–1655 nm) reflectance. Sentinel-2A and 2B operate on the same 10-day sun-synchronous orbital cycle, acquiring observations at approximately 10:30 a.m. local time with an offset of 180°. Sentinel-2 was designed to have a 5-day global repeat frequency, but due to the study area falling within a swath overlap area, imagery was collected more frequently (approximately every 2.5 days). Of the 29 Sentinel-2 observations during the study period, 15 (52%) were cloud-free.

Because Sentinel-2 does not provide an fSCA product, we estimated Sentinel-2 fSCA values from NDSI values using a hyperbolic tangent function following Gascoin et al.’s (2020) method [99] where:

$$fSCA=c(\tanh \left(a NDSI+b\right)+1)$$

(1)

equation 2 was fit using UAS-based fSCA values within all of the 20 m Sentinel-2 pixels and the Sentinel-2 NDSI observations for the six dates with coincident UAS and Sentinel-2 observations. fSCA values greater than 1 (100% snow-covered) were set to 1. Parameters a, b, and c were optimized by minimizing the MAE between the predicted fSCA and the training fSCA (i.e., the UAS fSCA observed within each Sentinel-2 pixel). The optimized parameters, a = 3.6, b = -2.0, and c = 0.56, had an MAE of 6.6% and a correlation coefficient (R²) of 0.91.

2.3. Analysis Methods

2.3.1. UAS-Based Snow-Covered Area Mapping

Snow cover maps with approximately 5 cm spatial resolution were created using k-means clustering [100] of the UAS RGB orthomosaics (15 orthomosaics in total, Figure 2). The k-means algorithm was run in MATLAB using built-in functionality based on the k-means++ random seeding algorithm [101]. After masking forest cover and open water, four classes were used to discern bare earth, ice, snow, and shadows. Areas having shadows were retained as a separate class because shadows made the underlying cover class difficult to discern. A manual review of the classification was performed to assign consistent classes across all dates. For example, when the entire field was snow-covered, the k-means algorithm identified multiple sub-classes with snow-covered areas. In these cases, the sub-classes were combined into a single snow class.

Snow cover maps produced from the UAS RGB orthomosaics using k-means (hereafter referred to as UAS-RGB-SCA) were validated using three 15 x 15 m test areas (~75,000 classified 5 cm pixels) in the field with varied shading, snow, and ice characteristics. Each test area was classified manually for six dates having snow cover at a different stage of melt (18 total segments). The test areas were carefully delineated by hand for the four k-means classes (bare earth, ice, snow, and shadow) and converted to raster format. Agreement between the manually delineated and k-means-based UAS-RGB-SCA classes was quantified using overall accuracy, defined as the number of agreeing classifications divided by total classifications, and F1-score. F1-score was calculated for each of the four classes as the harmonic mean of precision and recall as:

$$F1 score=2\times {\displaystyle \frac{Precision\times Recall}{Precision+Recall}}$$

(3)

where precision is the number of correctly classified positives (instances when the UAS-RGB-SCA and the manual classification agree on a given class) divided by all positive UAS-RGB-SCA classifications and recall is the number of correctly classified positives divided by all positive classifications in the manually delineated reference.

Field fSCA values from the UAS were calculated as the total area classified as snow, divided by the area of the field. Using this approach can lead to an underestimation of field fSCA when shadows are prevalent as the calculation excludes shadowed areas that contain snow. However, this influence was negligible for most dates (Table S1) Other limitations include the uncertainty created by the inclusion of the ice class, due to the similar optical properties of ice and snow.

2.3.2. Field-Aggregated Satellite and UAS Snow Cover Observation Comparisons

Field-averaged satellite and UAS observations of fSCA and NDSI were compared on the 15 UAS RGB flight dates. For each satellite product, values for dates without clear sky observations were estimated using a linear interpolation of the field-averaged values from the closest date prior to and after the UAS flight. Correlation analyses were performed for only clear sky observations collected on the same date and all dates including clear sky and the interpolated values.

2.3.3. Predictors of Field-Scale Snow Cover Persistence

Microtopography and vegetation influence snow accumulation and melting processes [72,102]. Ridges, valleys, slopes, or flat terrain have been consistently highlighted as a driver of snow depth spatial variability [47,50,66,67,103,104]. Forest edge effects also contribute to high variability [105,106] due to shading/solar exposure [29,47], canopy interception [50], and wind attenuation [107], which can be detected by UAS at fine scales [29,71,108].

We used the UAS LiDAR observations of bare earth (i.e., snow-off) from the fall before the snow season (December 14, 2020) to produce seven static spatial predictors for each 20 cm pixel (Figure 3A-G). An additional predictor, Sentinel-2 NDSI (Figure 3H), varied by date. Topographic features were calculated from DTMs and/or Digital Surface Models (DSMs). Using the GDAL software in QGIS we calculated four topographic predictors from the 20 cm resolution DTM as follows:

• Northness: The cosine of aspect, in which values range from 1 (north-facing) to -1 (south-facing)
• Roughness: The largest absolute inter-cell difference in elevation between a central pixel and its surrounding 8-neighborhood pixels in meters
• Topographic Position Index (TPI): A prominence index, calculated as the average difference between a central pixel’s elevation and its surrounding pixel elevations in meters for each 20 cm pixel using a 10 m moving window
• Slope: The average of the absolute value of all slopes from a central pixel relative to the surrounding pixels for each 20 cm pixel using a 20 m moving window

Analyses showed that 10 m (TPI) and 20 m (slope) moving windows were most predictive of snow cover persistence patterns at the field site (Figure S1). The remaining three static predictors were derived as follows:

• Grass Height: The elevation range (in meters) of UAS LiDAR returns within each 20 cm pixel
• Distance to Forest Edge: Euclidean distance in meters from all within-field 20 cm pixels to the forest edge
• Sun Hours: Using the UAS LiDAR DSM (1 m pixels), the TopoToolbox software package shadow function in MATLAB [109] was used in combination with the National Renewable Energy Laboratory’s Solar Position Algorithm [110] to produce a daily average of unshaded (sun) hours during the UAS flight period

2.3.4. Random Forest Snow-Covered Area Downscaling

RF modeling is a supervised machine learning technique that relies on an ensemble of decision trees to capture non-linear relationships across many input predictors [111]. When downscaling satellite observations of snow cover, RF models have proven to be more effective than neural networks and support vector machines [25]. From early applications in snow cover mapping using Landsat data with decision trees [18] to more complex RF models driven primarily by spatial variations in topography [75,102,112,113], tree-based classifiers accurately delineate SCA with low computational cost. Such models have been widely applied to enhance the spatial resolution (i.e., downscale) of snow cover datasets [25,75,102,113].

We trained a RF model using MATLAB’s TreeBagger function with the UAS-RGB-SCA observations and predictors shown in Figure 3. An ensemble of weak learners (i.e., decision trees) was produced by randomly sampling a subset of all predictors (i.e., bagging) at each node (i.e., decision rule or split point) and defining a split rule that best partitioned the dataset between UAS-RGB-SCA snow and bare earth classes. The eight predictor variables were normalized using the ranges shown in Figure 3 and resampled from their native resolutions (ranging from 20 cm to 20 m) using nearest neighbor interpolation to match the resolution of the reference UAS-RGB-SCA dataset (~5 cm). In training the model, six dates with variable snow cover conditions were used in which UAS observations coincided with clear sky observations from Sentinel-2. The initial training set contained nearly 10 million feature vectors. This was reduced to less than 3 million by drawing a random sample from each NDSI bin grouping (-0.55 to 1, with increments of 0.05), which ensured the model was trained based on a diverse but balanced distribution of NDSI values. Finally, the number of model input feature sets was reduced to approximately 300,000 by randomly sub-sampling 10% of the ~3 million feature sets. We confirmed that the distribution of predictors in the final feature sets captured the initial training feature set’s range of variability.

Three model parameters, the number of predictors to sample at each tree node (NumPredictors), minimum leaf size (MinLeafSize, or smallest number of features allowed in a single tree endpoint or leaf node), and the number of decision trees in the ensemble (NumTrees), were optimized considering the model accuracy computed on hold-out data (25% of the training set). NumPredictors was optimized as 3, which is approximately the square root of the number of input predictors (8). MinLeafSize was optimized to be 0.1% of the total dataset size. The difference between using 0.1% (~1000 features per split) and 0.001% (~10 features per split) improved model accuracy by only 2% and overfit the model. The model improvement due to increasing NumTrees plateaued after approximately 25, though 100 was used to ensure model robustness. Predictors were sampled using an in-bag fraction of 1, allowing all data from the selected predictors to be sampled for defining each split. Predictor importance (i.e., permuted delta error) was determined by randomly permuting all values of each predictor and determining the average loss of predictive skill (i.e., accuracy) across all individual trees within the random forest. The permuted delta error values were then normalized relative to the predictor with the highest error.

The RF model output was used to produce 5 cm resolution snow cover maps (RF SCA) over the field site for each of the six days having both Sentinel-2 and UAS-RGB-SCA observations. The RF model produces 100 decision trees with each tree providing a binary SCA value (1 = snow-on, 0 = snow-off). The final output from the RF model is the proportion of decision trees predicting snow for each 5 cm pixel (RF Proportional SCA). The RF Proportional SCA value is used in combination with the 20 m Sentinel-2 NDSI values to determine whether each of the 5 cm UAS pixels within a Sentinel-2 pixel is snow-covered or snow-free.

The modeling workflow, shown in Figure 4 (and detailed further in Figure S2), illustrates the steps used to produce downscaled snow cover estimates. First, the RF Proportional SCA value is determined for each 5 cm pixel. Then, the 20 m Sentinel-2 predicted fSCA value is determined from its NDSI value using the NDSI-fSCA conversion function (equation 2). Next, the RF Proportional SCA snow-on/snow-off threshold is set to a value between 0 and 1. The threshold is set such that the proportion of all (~160,000) 5 cm pixels within the 20 m Sentinel-2 pixel classified as snow-covered reproduces the Sentinel-2 predicted fSCA value for that pixel. Finally, each 5 cm pixel within the Sentinel-2 pixel is designated as snow-on if its RF Proportional SCA value exceeds the threshold and snow-off if its value is below the threshold. RF Proportional SCA thresholds had a large range (0.05 - 0.90), with the majority (73%) falling between 0.20 and 0.80. The model accuracy is computed first using a ~75,000 feature subset (25% holdout) of the data during the training process, then using a pixel-wise comparison between the applied model RF SCA value and the UAS-RGB-SCA reference.

2.3.5. Bare Patch Fractal Geometry

The fractal geometry of the emerging bare (snow-off) patches was analyzed for select RF SCA and UAS-RGB-SCA maps during the ablation period, February 20 – March 11, 2021. Snow and bare patches can be characterized as fractal objects where there is a power relationship between the patch area and patch perimeter across scales as well as between the patch area and the distribution of sizes [51]. If the patches behave as fractals, the patch perimeter (P) is a function of its area (A) as:

$$P={kA}^{{\displaystyle \frac{D_p}{2}}}$$

(4)

where k is a unitless scale-dependent multiplier and D_p is the fractal dimension. For natural objects, D_p falls between 1 (i.e., a straight line or single flat edge) and 2 (i.e., an infinitely convoluted and complex perimeter). Similarly, the fraction of the number of patches with a size equal to or greater than the area, F(A), is related to the patch area following Korcak’s law:

$$F\left(A\right)={cA}^{{\displaystyle \frac{{-D}_k}{2}}}$$

(5)

where c is a unitless scale-dependent multiplier, and D_k is the fractal dimension. Larger D_k values reflect a more uniform distribution of areas.

To calculate D_p and D_k for RF SCA and UAS-RGB-SCA, bare patches were identified and outlined by applying morphological operations to the bare earth class of each map. Edge maps of the bare earth class were then created through convolution with the Sobel kernel. Continuous bare patches were formed through morphological operators applied to the edge maps. These were used to dilate regions so nearly connected bare patches were merged. Finally, any interior holes within the dilated bare patches were removed, and an erosion operation was applied to reduce the patches to their original size. The slope of the resulting power relationship between the perimeter or areal frequency distribution and area of each outlined bare patch was calculated for each survey date.

3. Results

Observed in-situ meteorological data for the study period, December 2020 through April 2021, are shown in Figure 5. The largest single-day precipitation (37 mm) occurred on January 16, 2021, in a rain-dominant event (daily average air temperature of 5 °C). The peak snow depth of 21 cm occurred following a series of precipitation events between February 2nd (-7 °C) and 9th (-3 °C), 2021. The average daily incoming solar radiation ranged from clear sky values of ~400 W/m² in December to ~800 W/m² in April.

3.1. UAS-Based Snow-Covered Area

Three UAS SCA test areas were used to evaluate the accuracy of the k-means-based UAS RGB snow cover estimates (UAS-RGB-SCA). Figure 6 shows the UAS-RGB-SCA products for four example dates during the ablation period, February 20 – March 11, 2021. Side-by-side comparisons of 5 cm RGB orthomosaics, UAS-RGB-SCA, manually delineated classes, as well as the agreement between products, are shown for the three test areas on each date. The UAS-RGB-SCA product was comparable to the analyst-delineated SCA reference capturing spatial patterns. Overall, the UAS-RGB-SCA classification performance was strong, with the shaded north test region having a reduced agreement with the manual classification of 80% compared to the 86% and 88% in the central and southern test sites, respectively (Table 1). Days with the most homogenous surface cover conditions (e.g., February 20-28 and March 11, 2021) had the highest agreement. The agreement ranged from 100% on February 20, 2021, when the south test region was completely snow-covered to a minimum of 33% on March 7, 2021, in the north test area. The reduced performance in the north test area on March 7th is primarily due to the difficulty in distinguishing shadows from ice, snow, and bare earth classes using only visible imagery. March 7th had a large relative proportion covered by shadows and ice (63%) in the northern test area (Table 1). F1 scores calculated for snow detection showed similar results, with the highest scores (0.85 – 1) during periods with majority snow cover between February 20th and March 3rd, 2021, and the lowest (0.19 – 0.30) occurring on March 11th when snow cover was less than 5% in each test area.

Having established that k-means clustering produces good agreement with manually delineated SCA, it was used to produce UAS-RGB-SCA maps and classify their persistence for all fifteen dates during the study period (Figure 7) as well as their patchiness on four dates during snowmelt (Table 2).

Prior to ablation, there were already over 2,000 small patches whose area totaled less than 1% of the study area. The patches’ size and count increased during the melt period, with the greatest number of patches in early March. We found that the snow patches exhibited self-similarity over fine to moderate scales (i.e., 0.01 m² to 100 m²) with the area-perimeter fractal dimension (D_p) ranging from 1.24 to 1.40 and area-frequency (D_k) fractal dimension from 1.55 to 1.73. While most of the snow had melted by early March, persistent snow was observed in protected areas north of forest edges into early April (Figure 7C). Lower persistence was observed in areas having high solar exposure, including areas south of forest edges and south-facing slopes (e.g., the south-central field near the roadway). Areas with taller grasses also had lower snow persistence, though higher persistence was frequently observed in adjacent depressions (e.g., west of the roadway near the USCRN site). Similarly, the effect of grass being mowed in rows, which can be seen in the TPI metric in the northeastern field (Figure 3C), resulted in snow persistence varying in organized patterns over short distances. The shaded area at the north-facing forest-field boundary in the western part of the field was correctly classified as a shadow on multiple days but was also incorrectly classified as bare earth on others (Figure 7A). This reduced the calculated snow persistence at that location. Regions where water pooled and re-froze, forming persistent surface ice, were also discerned in numerous locations (Figure 7B).

3.2. Fractional Snow Cover from UAS and Satellite Observations

In this section, the fifteen UAS-RGB-SCA observations are compared to the satellite snow cover metrics and the UAS-LiDAR-SCA using the average fSCA values for the study site. Results indicate that all UAS and satellite products captured the general pattern of snow accumulation and ablation when averaged over the field area (Figure 8). Even products with a coarse spatial resolution (MODIS/VIIRS) containing mixed forest-field pixels effectively capture fSCA over the field. The observational counts during the snow-on period were limited for Landsat (3 total snow-on observations) as compared to Sentinel-2 (9), MODIS (10), VIIRS (15), and UAS-RGB-SCA (15). The snow accumulation period (February 1–7, 2021) had four UAS flights, but no cloud-free satellite observations. In contrast, the two-week melt period, characterized by relatively warm, clear sky days and cloud-free satellite observations, is well captured by the entire suite of products. While all satellite products indicated that the study area was snow-free between March 10th and 12th, 2021, UAS-RGB-SCA showed persistent snow in protected areas (fSCA < 10%) on March 10th and 11th.

For the field averaged values, all snow cover products and the reference UAS-RGB-SCA had a statistically significant linear relationship (α = 0.05) (Figure 9), indicating all products captured snow cover variability. The correlation between VIIRS, MODIS and Sentinel-2 NDSI and the observed UAS-RGB-SCA fSCA exceeded 0.90 on coincident dates. However, the range of the average NDSI values varies by product due to different pixel footprints, sensor properties, and orbital periods. The correlation between the derived Sentinel-2 fSCA values and the observed UAS-RGB-SCA fSCA also exceeded 0.90. Landsat fSCA and NDSI (Figure 9C and 9D) also had statistically significant correlations to the UAS-RGB-SCA reference but had lower correlation values of 0.73 and 0.84, respectively, despite having no coincident values.

The moderate agreement between fSCA from UAS-LiDAR and UAS-RGB-SCA (R² = ~0.80) indicates that UAS-LiDAR snow depths provide an imperfect estimate of fSCA even with low UAS-LiDAR snow depth errors (~ 2 cm) (Figure 9A). The LiDAR observations also provide insight into the relationship between snow depth and fSCA from the UAS-RGB-SCA (Figure 9B). At this field site, there was a clear positive linear relationship between the fSCA and field average snow depth (Figure 9B), especially for shallow depths (< 10 cm) and UAS-RGB-SCA fSCA below 0.80. As the average snow depth at the site decreases below 10 cm, snow cover becomes more discontinuous. This suggests that UAS fSCA and snow water equivalent (SWE) derived from UAS LiDAR snow depths could be used to parameterize snow depletion curves once bare patches emerge.

3.3. Random Forest Modeling of Snow-Covered Area Using Surface Features

Sentinel-2 was chosen to demonstrate the capacity of a moderate-scale SCA satellite product to capture snow cover within the field area. MODIS’s and VIIRS’s coarse resolutions are not adequate to capture the spatial distribution of snow presence within this small study site. Landsat’s 30 m and Sentinel-2’s 20 m moderate resolutions both capture the spatial distribution of within-field SCA and largely avoid the need to interpret mixed forest-field pixels. However, Sentinel-2 has a higher spatial resolution (20 m versus 30 m) and temporal repeat (3 days versus 7 days) than Landsat 7/8.

Thus, the random forest (RF) SCA pattern model was trained using the Sentinel-2 NDSI product and the static physical metrics. Sentinel-2’s NDSI was the most important predictor of snow cover in the RF model and the only dynamic metric included (Figure 10). The static physical metrics that had relatively high importance are the distance to forest edge, slope, TPI, and sun hours. Roughness, northness, and grass height provided limited predictive value in modeling fine-scale snow variations. The RF model was used to produce downscaled 20 m Sentinel-2 SCA values. The RF model SCA results for the 25% holdout data had an overall agreement with the reference UAS-RGB-SCA classifications of 88.9%, with 92.9% and 81.3% agreement for snow-off and snow classes, respectively.

Comparisons among the UAS-RGB-SCA, Sentinel-2, and the RF SCA values were made for the six days having clear sky UAS and Sentinel-2 observations. These six dates have a range of snow cover conditions, transitioning from nearly complete snow coverage in late February, to sporadic snow cover in March, to snow-off conditions by early April (Figure 11). Qualitatively, the three SCA products have similar snow cover variations in the field with the primary difference being the UAS and RF model’s high resolution and Sentinel-2’s moderate resolution. The RF Proportional SCA values readily identify regions with early melt on February 26, March 3, 8, and 10 along the northern edge of the field adjacent to coniferous trees. While Sentinel-2 indicates that these areas had reduced NDSI, the RF Proportional SCA output was able to distinguish which portions of the 20 m pixel had snow and to translate those results to a 5 cm SCA product. During the 11-day melt period, the snow reduction at these forest-field transitions is clear in the RF SCA product and shows a reasonable progression of expanding bare areas.

A quantitative assessment of Sentinel-2 and the RF SCA products compared to the UAS-RGB-SCA reference was conducted (Table 3). Sentinel-2 pixels having fSCA values that exceeded a threshold of 0.5 (50% snow-covered), which corresponded to an NDSI of >0.52, were classified as snow-on with all other locations classified as snow-off. The best performance occurred on sampling dates with nearly completely snow-covered or snow-free conditions. There was over 82% agreement between the UAS-RBG-SCA and Sentinel-2, except on March 8, when conditions were a mixture of snow, ice, and bare earth. The RF SCA agreement with UAS-RGB-SCA was comparable to Sentinel-2 during periods when the site was dominated by snow or bare earth conditions but exhibited an improvement of 4 to 8% when conditions were mixed (February 26 – March 8). The lowest agreement between RF SCA and the UAS-RGB-SCA was approximately 73%, occurring on March 8th when snow cover was highly discontinuous. In such mixed conditions, satellite-derived snow cover did not reliably capture surface conditions. However, for all approaches, when one cover class was dominant, the minority cover type was poorly captured (< 62%, Table 3). For example, on March 10th, most of Sentinel-2’s snow-off pixels matched the UAS-RGB-SCA bare earth class (97%), but Sentinel-2 SCA also misclassified 73% of the UAS-RGB-SCA snow or ice pixels as snow-off. RF SCA was able to reduce the misclassification to about 59%.

While the Sentinel-2 SCA and RF SCA had very similar agreements with UAS-RGB-SCA, Figure 11 shows a marked difference in patchiness. The Sentinel-2 product, which cannot resolve any snow cover variations at scales finer than its 20 m pixel resolution, had mean patch areas that exceeded 500 m² and at most ten bare patches in the field on any date. Fractal dimensions could not be calculated for Sentinel-2 due to the limited number of unique patch sizes. The RF SCA product’s 1,500 to 7,600 bare patches were much higher than Sentinel-2, but somewhat lower than the UAS-RGB-SCA’s 2,300 to 29,000 patch counts (Table 2). The RF SCA patches were between two and five times larger than the UAS-RBG-SCA’s mean patch areas on all dates. RF SCA had strong linear power relationships over five orders of magnitude (0.01 m² to 100 m²), indicating that the RF model replicated the fractal nature of the field site’s bare patches. The average fractal dimension from the RF SCA model’s D_p was nearly identical (1.30) to the UAS-RGB-SCA observed values, indicating that the model produces similar complexity bare patch perimeter to that which was observed. The UAS-RGB-SCA’s has a higher D_k value (1.64) than the RF SCA model (1.40), indicating that the model has a less uniform distribution of patch areas with fewer small patches resulting in fewer patches overall.

In summary, the use of fine-scale vegetation or terrain information from snow-off UAS-LiDAR observations in the random forest model effectively downscaled Sentinel-2 NDSI observations to sub-meter resolutions with modest changes to overall accuracy (-0.5 to +7.9%) and considerable enhancement in capturing the number and size of the bare ground patches. The RF SCA model was also able to qualitatively capture locational differences in melt timing. The first regions to melt out were the exposed forest edge at the northern site boundary (i.e., high solar exposure and high canopy snow interception) and the steeper, more south-facing slopes in the southwestern field (i.e., high solar exposure). The areas with low solar or wind exposure, such as the north-facing forest edges and the small surface depressions, were observed and modeled to melt out last. However, we acknowledge reduced model performance in these areas due to frequent shading, which led to some errant UAS-RGB-SCA bare earth classifications during periods of snow cover (see February 26 and March 8 in Figure 11). Even small snow-covered areas on the forest edge and in the depressions between rows of grass were effectively predicted during the late stages of melt (see March 10 in Figure 11).

4. Discussion

4.1. Performance of UAS-Based Snow-Covered Area Mapping

Few studies have reported performance for snow cover products from UAS optical imagery which may be due to challenges, subjectivity, and variability in manual classification depending on the operator and site conditions [26,27,28,29]. Our study’s 5 cm UAS-RGB-SCA product had an 86% overall accuracy compared to the manual classification. This was somewhat lower than Eker et al.’s (2019) overall accuracy which exceeded 95% under ideal conditions (i.e., flat, unvegetated site) [27]. In Northern Arizona, Belmonte et al. (2021) reported a 90.2% overall accuracy for snow cover from 15-cm resolution multispectral (green, red, red edge, and near-infrared or NIR) UAV orthomosaic images compared to manually digitized polygons [29]. Niedzielski et al. (2018) mapped SCA using both UAV RGB and RGB combined with NIR in the Izerskie Mountains in southwestern Poland [26]. While they did not report traditional agreement metrics, both approaches had about a 5% error in the areal means compared to manually classified snow cover, indicating that their snow mapping performance was not notably improved using NIR. Similar to Niedzielski et al. (2018) [26] and Malle et al. (2019) [28], we found that it was difficult to distinguish between shadows and bare earth, as well as ice and snow, due to the low optical contrast in RGB imagery. Overall, preliminary results from UAS snow cover studies suggest that, compared to satellite snow cover mapping, UAS snow cover mapping is less impacted by cloud cover, forested areas, complex terrain, and shallow snow conditions (< 10 cm depth). In cases with close proximity to tall vegetation and or steep terrain, cloudy conditions may be preferable to clear sky conditions because cloud cover provides diffused light, which reduces the presence of harsh shadows.

It is challenging to generalize performance from previous studies due to differences in SCA classification methodology and performance statistics. While satellite approaches typically use a Normalized Difference Snow Index method (NDSI, [114]) to map snow cover using green and SWIR bands, there is no single recommended strategy for using UAS RGB imagery to map snow cover. In this study, we followed Niedzielski et al.’s (2018) k-means clustering approach [26] to distinguish among snow-on, snow-free, shaded, and icy regions. A strength of k-means clustering is that it does not require training data. However, each image needs separate, manual interpretation of the resulting classes. Other promising methods for differentiating between snow and non-snow in UAS RGB imagery include blue-channel thresholding, random forest supervised classification, and supervised classification via image processing software (e.g., ENVI), all of which require training data [24,26,27,29,115].

4.2. The Role of UAS in Snow-Covered Area Mapping

Despite increasing interest in UAS applications, questions remain on how to best use these systems to advance cryosphere science [33,35]. This research contains three examples of how UASs can be incorporated with satellite observations to support snow cover mapping: (1) supporting validation and calibration of satellite snow cover products across spatial scales, (2) translating satellite indices (e.g., NDSI) to physical variables (e.g., fractional SCA), and (3) developing models to downscale satellite SCA observations using UAS-observed land and vegetation structure (via LiDAR) and detailed imagery (via optical/RGB sensors).

4.2.1. Evaluating Satellite Snow Cover Products

Traditionally, coarse-scale snow cover maps have been validated using high-moderate resolution imagery such as that from Landsat [89,116], SPOT [117], ASTER [94], and/or imagery collected from piloted aircraft [18]. To our knowledge, only three other studies have linked UAS and satellite remote sensing observations of SCA [24,25,115]. Liang et al. (2017) trained MODIS fSCA models using UAS-observed 5 cm resolution snow cover collected over 15 MODIS pixels [24]. Similarly, Liu et al. (2020) trained a random forest fSCA model using MODIS reflectance, elevation, and UAS snow cover [25]. Jacques-Hamilton et al. (2025) used UAS fSCA to evaluate MODIS, VIIRS, and Sentinel-2’s fSCA and snow melt timing and found that satellite approaches typically overestimate fSCA, with biases up to 60% for the coarser MODIS and VIIRS products [115]. In contrast to their results, we found that all averaged satellite products captured the evolution of fSCA when using clear sky observations (R² = 0.93 to 0.98) and even our linearly gap-filled satellite data had only modest reductions in correlations relative to using clear sky observations alone (-0.04 to -0.29). The gap-filled correlation would likely improve using more sophisticated SCA gap-filling approaches [30,31,34,118,119,120], especially when transitioning from snow-off to snow-on conditions.

Our finding that gap-filled satellite datasets induced error is not surprising because cloud cover is widely noted as the greatest source of uncertainty in the optical detection of snow-covered areas [32,86,91,121]. Due to the seasonal link between fSCA and SWE, regular cloud cover can also degrade the accuracy of SWE reconstructions [122]. UAS’s ability to make repeated, high-resolution snow cover maps during cloudy conditions has potential value for understanding uncertainties in satellite SCA estimates induced by cloud cover and quantifying temporal uncertainties in SWE estimates [123,124]. While the frequency of cloud cover differs widely by location [3], the seasonal variation can also be significant [125,126]. For the Southeastern New Hampshire, USA fields mapped in this study, satellite remote sensing observations of SCA were difficult to obtain during the cloudy accumulation season (68% cloud cover in mid-January to late-February) as compared to the snow depletion period (47%, late-February to early-April).

4.2.2. Translating NDSI to Fractional Snow Cover

The simple linear relationship between NDSI and fSCA is well-established and validated [127,128]. Historically, regression coefficients have been optimized for different sensors and used as the standard algorithms for MODIS fSCA products [89,129]. Recently, non-linear sigmoidal models have been developed to better capture the shape of the relationship between NDSI and fSCA [99] and quantify uncertainty in VIIRS snow cover estimates [121]. Our research demonstrated that UAS observations can be used to produce accurate, locally tuned functions for converting satellite NDSI observations to fSCA. Our optimized sigmoidal function using the ~5 cm UAS SCA observations within each Sentinel-2 pixel yielded strong performance (R² = 0.91, MAE = 6.6%) relative to the UAS-observed snow cover fraction. Our sigmoidal model also had better performance in high and low fSCA conditions (Figure 9 E) compared to linearly rescaling the range of observed NDSI values (Figure 8).

4.2.3. Downscaling Satellite Snow-Covered Area Products

A final application for UAS-RGB-SCA observations is training machine learning models to downscale satellite SCA observations. Previous studies have demonstrated the utility of decision trees and RF models for improving snow cover detection from satellite data at coarse and moderate spatial resolutions and often use higher resolution SCA (e.g., via Sentinel-2 or Landsat) with vegetation and terrain features as model inputs [18,47,75,102,130,131,132,133]. UAS optical imagery for snow cover modeling are orders of magnitude finer (0.05 m) than even the highest spatial resolution satellite products that are available publicly (10 m) and commercially (1 to 3 m). While UAS imagery has not been widely used to downscale snow cover, it has been applied to produce fSCA estimates for full MODIS pixels (500 m) using machine learning [24,25]. Our random forest model (RF SCA) produced realistic representations of snow cover at cm-level detail with an overall accuracy (90.8%) that is slightly higher than the binarized Sentinel-2 SCA (87.4%).

The most notable performance improvements were during periods with patchy cover (March 3, 8, and 10) when RF SCA accuracy values were higher by 5 to 8% and more realistic snow patches were produced as compared to Sentinel-2. Our RF SCA fractal dimensions (D_p = 1.30 and D_k = 1.40) were comparable to those found from the UAS imagery (D_p = 1.29 and D_k = 1.64) and similar to two other studies documenting the fractal behavior of snow patch area-perimeter (D_p = 1.46 to 1.68) and area-frequency relationships (D_k = 1.21 to 1.79) [51]. Their work used moderate-scale pixels over a larger study area to show self-similarity over scales from 1 m² to 10,000 m², which overlaps with our scales of 0.01 m² to 100 m². Despite the evidence of self-similarly, it is possible that there are multiple physical processes operating over these scales that, in combination, produce the observed fractal behavior [79]. The detailed, fine-scale quantification of the overall snow cover fraction, sub-grid scale snow cover fraction, and bare ground and snow patches are essential for improving energy balance modeling and parameterization when snow is patchy [54,133,134,135]. UAV and field camera optical imagery datasets from ablation periods are well poised to provide these characteristics and to train satellite downscaling models.

The RF model development demonstrated here also provides insights into the fine-scale physical controls on snow cover depletion. Our finding that Sentinel-2’s NDSI was the most important input for downscaling agrees with Rittger et al.’s (2021) conclusion [102] that the most important predictor variables were MODIS fSCA and day of year because they capture the temporal evolution of snow cover. Yet, static variables (i.e., distance to the forest edge, slope, TPI, and solar exposure) are needed to predict fine-scale SCA patterns in the RF model. Not surprisingly, these static variables agree with Clark et al.’s (2011) summary [72] of a legacy of studies documenting terrain and vegetation controls on snow spatial variability and the importance of scale. While all static variables were defined at high spatial resolutions, the most important variables considered the influence of adjacent pixels, as opposed to variables (e.g., surface roughness, aspect, and grass height) that only measured conditions for individual 20 cm pixels. This suggests that there may be a declining value to using UAS observations at scales finer than 20 cm.

5. Conclusions

Snow cover persistence is useful for detailing snow spatial variability over large areas, and has direct relevance to agriculture, hydrology, mobility, and land surface modeling. This work provides a case study analysis of snow cover estimates across observational scales ranging from UAS’s high resolution (5 cm) to MODIS’s coarse resolution (500 m) observations and demonstrates how UAS may be used to validate satellite snow cover products, convert satellite NDSI observations to fSCA, and downscale satellite SCA observations. SCA derived from UAS RGB imagery achieved more than 80% accuracy in almost all test conditions, including those with variable illumination. UAS RGB performance, even at fine scales (<1 m), was similar to satellite snow cover products despite the lack of a NIR band, and performance was improved relative to satellite products when clouds were present or under shallow snow conditions.

Fine-scale UAS products are readily aggregated to moderate or coarse scales, and the site average values provide confidence that satellite fSCA and NSDI captured the general snow cover evolution during ablation. However, because MODIS and VIIRS were limited to four and five pixels, respectively, within the study area, only Landsat and Sentinel-2 could show spatial patterns within the field (i.e., at field scales, <500 m). As compared to Landsat, Sentinel-2’s more frequent observations were needed to capture the brief, two-week ablation period for shallow snowpacks that are typical of ephemeral regions such as our Southeastern New Hampshire, USA site. UAS-RGB-SCA was able to discern that the 20 cm random forest-based satellite SCA downscaling model improved the overall accuracy by 5 to 8% and created more realistic bare ground patches as compared to the traditional Sentinel-2 NDSI thresholding produced SCA maps during ablation period.

This study’s findings point to emerging opportunities for UAS RGB imagery during winter. There is already a wealth of UAV and field camera RGB datasets with high-resolution snow information, but relatively few studies have explored or applied these datasets. This case study demonstrates some ways in which UAS optically-sensed snow cover maps can contribute to cold weather research and applications and recommends future research avenues that incorporate these datasets for providing high-resolution snow cover information.