Next Article in Journal
The Effect of Forest Mask Quality in the Wall-to-Wall Estimation of Growing Stock Volume
Next Article in Special Issue
Drone-Based Hyperspectral and Thermal Imagery for Quantifying Upland Rice Productivity and Water Use Efficiency after Biochar Application
Previous Article in Journal
Dynamics of Vibrio cholerae in a Typical Tropical Lake and Estuarine System: Potential of Remote Sensing for Risk Mapping
Previous Article in Special Issue
Hyperspectral and Thermal Sensing of Stomatal Conductance, Transpiration, and Photosynthesis for Soybean and Maize under Drought
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 13
Issue 5
10.3390/rs13051036
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

UAV-Based Estimate of Snow Cover Dynamics: Optimizing Semi-Arid Forest Structure for Snow Persistence

by
Adam Belmonte
1,
Temuulen Sankey
1,*,
Joel Biederman
2,
John Bradford
3,
Scott Goetz
1 and
Thomas Kolb
4
1
School of Informatics, Computing, and Cyber Systems, Northern Arizona University, 1295 S. Knoles Ave, Flagstaff, AZ 86011, USA
2
USDA Agricultural Research Service Southwest Watershed Research Center, Tucson, AZ 85719, USA
3
Geological Survey Southwestern Biological Center, Flagstaff, AZ 86001, USA
4
School of Forestry, Northern Arizona University, Flagstaff, AZ 86011, USA
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(5), 1036; https://doi.org/10.3390/rs13051036
Submission received: 20 January 2021 / Revised: 10 February 2021 / Accepted: 4 March 2021 / Published: 9 March 2021
(This article belongs to the Special Issue Ecohydrological Remote Sensing)
Browse Figures
Versions Notes

Abstract

:
Seasonal snow cover in the dry forests of the American West provides essential water resources to both human and natural systems. The structure of trees and their arrangement across the landscape are important drivers of snow cover distribution across these forests, varying widely in both space and time. We used unmanned aerial vehicle (UAV) multispectral imagery and Structure-from-Motion (SfM) models to quantify rapidly melting snow cover dynamics and examine the effects of forest structure shading on persistent snow cover in a recently thinned ponderosa pine forest. Using repeat UAV multispectral imagery (n = 11 dates) across the 76 ha forest, we first developed a rapid and effective method for identifying persistent snow cover with 90.2% overall accuracy. The SfM model correctly identified 98% (n = 1280) of the trees, when compared with terrestrial laser scanner validation data. Using the SfM-derived forest structure variables, we then found that canopy shading associated with the vertical and horizontal metrics was a significant driver of persistent snow cover patches (R2 = 0.70). The results indicate that UAV image-derived forest structure metrics can be used to accurately predict snow patch size and persistence. Our results provide insight into the importance of forest structure, specifically canopy shading, in the amount and distribution of persistent seasonal snow cover in a typical dry forest environment. An operational understanding of forest structure effects on snow cover will help drive forest management that can target snow cover dynamics in addition to forest health.

1. Introduction

Runoff from seasonal snowpack in semi-arid regions provides drinking water and agricultural irrigation for at least one-sixth of the world’s population [1]. Snowpack also provides important ecosystem services, including water for vegetation, aquatic ecosystems, and shallow groundwater recharge [2,3,4,5,6]. Throughout the western United States, higher annual mean temperatures and changes to winter air humidity are contributing to an earlier and faster spring snowmelt [7]. Regional climate projections indicate that these effects will become more pronounced [8,9,10] and inconsistent snowpack will contribute to higher rates of drought-induced tree stress and mortality [11].
In the ponderosa pine (Pinus ponderosa) forests of northern and central Arizona, roughly half of the annual precipitation falls as rain during the summer, and the other half as snow during the winter [12]. For mature ponderosa pine trees, winter precipitation via snowpack is the dominant source of water throughout the year, underscoring the importance of snowpack to forest health [13,14]. In turn, forest cover influences the accumulation and ablation of snow on the ground surface, primarily by controlling canopy interception and by partitioning the solar radiation available at the ground surface [15,16,17,18,19,20]. Additionally, the size, shape, spacing, and structure of forest patches (i.e., groups of trees) influence snow accumulation and ablation processes [21,22,23,24,25].
At the landscape-scale, snowpack distribution, depth, and snow water equivalent (SWE) are most accurately estimated using airborne lidar datasets. Similarly, recent advances in forest inventorying and monitoring have leveraged both terrestrial and airborne lidar datasets to produce georeferenced and scaled measurements of individual trees at the landscape scale (400+ ha) [26,27,28,29,30,31,32,33,34]. However, lidar is often expensive to acquire. Optical remote sensing of snow-covered area (SCA) is a proven method that can provide accurate estimates of spatial and temporal snowpack accumulation and distribution [35,36,37]. However, SCA estimates lack the snow depth and SWE dimensions. In addition, airborne and satellite-based optical images generally have coarse spatial resolutions, fixed time intervals, and cloud interference that limit their applications in the southwestern USA, where snow rapidly melts in <2 weeks following a storm [37]. It is important to identify alternative means to rapidly and cost effectively monitor landscape-scale snow extent and persistence for assessing the relationship between SCA and forest structure changes.
As a lower-cost alternative to airborne data, unmanned aerial vehicles (UAVs) are increasingly used to acquire high resolution images [38,39,40,41,42,43,44,45,46,47]. UAV-borne photogrammetric Structure-from-Motion (SfM) modelling has also been successfully used in fine (0–4 ha) and mid-scale (4–400 ha) forest structure assessments [39,40,48,49,50]. Compared to lidar-derived 3D point cloud models, UAV SfM often provide similar accuracies in estimating individual tree canopy heights, diameter, and volumes [41,49,50]. In addition to characterizing forest canopy cover, tree canopy height and diameter, tree density, patch and gap sizes and geometry [38,39,40,41], UAV datasets can also be used to measure snow cover and depth [51,52,53]. UAVs can provide very high resolution and near real-time SCA data [54] and can be examined with similar methodologies used in airborne and satellite optical remote sensing to derive SCA. UAVs offer a flexibility essential for capturing both ongoing forest structure changes as well as day-to-day variability of the highly dynamic snow cover in semi-arid Southwestern forests [39,41,55].
Here we quantify snowpack dynamics in a recently thinned and unthinned ponderosa pine forest to evaluate the potential benefits of the forest restoration treatment for snowpack. We use UAV-derived multispectral orthomosaic imagery to quantify snow cover dynamics and identify persistent snowpack, and three-dimensional Structure-from-Motion (SfM) models to examine the effects of forest structure, specifically via shading, on snowpack persistence. Our objectives were to:
  • Quantify snow cover following three different winter storms and identify persistent snowpack across the study site;
  • Examine forest structure shading effects on snowpack persistence;
  • Model and predict persistent snowpack using the most influential forest structure metrics.
To assess the importance of specific characteristics and their spatial distribution, forest structure metrics were separated into groups that emphasize vertical and horizontal characteristics. We hypothesized that horizontal forest structure metrics would be more influential for predicting persistent snowpack than vertical forest structure metrics due to the variation in ground shading exhibited by trees, no matter their size or shape. In addition, we hypothesized that crown base height and crown volume would be important for predicting the size of persistent snow patches.

2. Materials and Methods

2.1. Study Area

Our study area is located in the Coconino National Forest in northern Arizona (Figure 1). The study area includes 76 ha of forest, characterized by relatively flat topography ranging in elevation between 2200 and 2275 m, with slopes of 0–10% across most of the study area. The vegetation is characterized by ponderosa pine (Pinus ponderosa), which dominates both the region and the study area, and includes sporadic Gambel oak (Quercus gambelii) and Rocky Mountain juniper (Juniperus scopulorum). The climate is sub-humid and characterized by distinct seasonal trends including early summer drought, with an average of 560 mm of precipitation (Western Regional Climate Center, 2020) falling half as snow during winter months and half as rain during the mid-summer monsoon season. The site had been undisturbed since its last naturally occurring fire in 1876, except for selective historical firewood harvesting [56]. However, the site was subjected to a prescribed fire in 1976, in which 63% of the smaller surface fuels and 69% of the woody surface fuels (up to 8 cm in diameter) were consumed [56].
A mechanical thinning restoration treatment at the study area began during the fall of 2017 and was completed by the spring of 2018, yielding a 53-ha treated portion and 23-ha untreated portion for this study (Figure 1). Similar to other regional restoration treatments that promote diversity in tree group and interspace size, shape, and spacing, this treatment aimed to reinstate pre-settlement forest conditions and included a range of thinning goals that would promote healthy overstory vegetation and the regeneration of understory vegetation [56,57,58,59]. Specifically, the restoration treatment prescription at our study area emphasized irregular tree group delineation, expansion of interspace, retention of all non-ponderosa pine species (e.g., Gambel oak and juniper), and significant reductions in smaller ponderosa pine trees within groups and interspaces.

2.2. Data Description and Processing

2.2.1. Snow Covered Area and Persistent Snow Patches

Our examination of the relationship between persistent snow cover and forest structure began with quantifying snow covered area (SCA) across the entire study site throughout the melt-off period following three significant (>20 cm new snowfall depth) independent winter storm events during the winters of 2017–2018 and 2018–2019. For each storm, a set of 15-cm resolution multispectral UAV orthomosaic images were acquired from the first day following the storm (i.e., maximum SCA) until the snow cover had fully disappeared (i.e., minimum SCA). To capture temporal changes in SCA following each storm, the interval at which images were acquired depended on the weather conditions in the preceding days, but we attempted to keep image date intervals consistent among independent storm events (Table 1). The final SCA dataset consisted of three ‘snow-series’ (one snow-series per storm event), yielding a total of 11 different orthomosaic images.
Data comprising the snow-series orthomosaic images were collected using a Sensefly eBee RTK fixed-wing UAV platform (SenseFly, Lausanne, Switzerland) fitted with a Parrot Sequoia multispectral sensor (Parrot Drones SAS, Paris, France). Each flight mission was pre-programmed using Sensefly’s proprietary software, which controlled flight plans and customized all flight and data parameters. Flight specifications included high latitudinal and longitudinal overlaps (85% and 90%, respectively) in order to generate the SfM outputs, an average flight altitude of 120 m, and operation centered around solar noon on each image date. There was an average of ~350 images captured for each flight based on the pre-determined flight path, with wind conditions largely responsible for the total number of images collected. During each flight, four georeferenced images were recorded at each photo location in the green (530–570 nm), red (640–680 nm), red edge (730–740 nm), and near infrared (770–810 nm) spectral bands. These images were post-processed using Sensefly’s eMotion 3 software, which automatically excluded a few distorted images from each flight (SenseFly, Lausanne, Switzerland).
The images were then used to create the final datasets in Agisoft PhotoScan v1.4.0 photogrammetric processing software (Agisoft LLC, St. Petersburg, Russia) for each image date. All images from each flight were scanned for matching ‘tie-points’, oriented in three-dimensional space via bundle-adjustment, and then mosaicked together [39]. The general workflow in the software includes image alignment to create a sparse point cloud, incorporation of GCP locations, image alignment optimization, gradual filtering out of inaccurate and error-inducing points, and a full image realignment [39]. This workflow generated a final orthomosaicked image in 15 cm spatial resolution in all four spectral bands and dense 3-dimensional (3D) point cloud data photogrammetrically generated from the high resolution images using Structure-from-Motion (SfM) algorithms. We optimized the SfM algorithm parameters based on our previous results [39].
The snow-series images were used to first classify SCA across the study site in each image and subsequently delineate persistent snow cover patches. All images (n = 11) were classified into five classes: vegetation, sunlit bare ground, shaded bare ground, sunlit snow, and shaded snow, using a random forest supervised classification performed in Google Earth Engine [60,61]. The random forest classifier was parameterized based on previous remotely sensed image classification examples, with the number of trees set to 100 and the number of variables per split set to the square root of the number of variables [60,61]. The training dataset consisted of manually digitized polygons (n = 492 total) for each of the five classes. The image resolution (15 cm/pixel) allowed for the training dataset classes to be easily delineated. The final binary snow versus non-snow rasters were created by combining pixels classified as both sunlit and shaded snow into a single ‘snow’ class, while pixels in the remaining classes were combined into a single ‘non-snow’ class. The accuracy of this binary classification was assessed using a set of randomly generated snow/non-snow samples (n = 500) from the image with roughly half its pixels snow covered, allowing for an unbiased accuracy assessment. Image classification relied on a multi-band image stack created from the four original spectral bands (green, red, red edge, near infrared) and six additional bands derived from the original four: normalized difference vegetation index (NDVI) [62], soil adjusted vegetation index (SAVI) [63], Gray Level Co-Occurrence Matrix (GLCM) [64] variance, GLCM homogeneity, GLCM contrast, and GLCM entropy. The GLCM texture bands were included to increase the effectiveness of the classification algorithm in discriminating the different combinations of bare ground, snow, and their shaded equivalents.
Finally, persistent snow cover was delineated using a single snow cover image composite. This image composite was created by stacking each binary classification (n = 11) and counting when each pixel was snow-covered out of the 11 total dates in the image composite. A simple post-processing procedure was used to eliminate spurious pixels and reduce noise along edges of snow pixels. Specifically, we conducted a progressive moving-window majority filtering using windows of 3 × 3 and 5 × 5 cells. From this image composite, persistent snow patches were identified by selecting isolated groups of adjacent pixels with snow cover in 10 or 11 out of the total 11 images (Figure 2). Groups of pixels less than 200 m2 were numerous and observed to melt completely between image dates, thus they were eliminated from the analysis. The resulting patches were refined using a general boundary cleaning procedure to eliminate spurious edges, then vectorized and attributed with an identification number and an area estimate (m2).

2.2.2. Forest Structure Metrics

Forest structure metrics were quantified across the entire study site using the UAV SfM point cloud data and validated using both field measurements and a terrestrial laser scanner (TLS) point cloud dataset. The vertical forest structure metrics examined in our study are tree crown height (CH), crown diameter (CD), crown base height (CBH), the ratio of tree base height to crown height (CBH:Z), and crown volume (CV). The horizontal structure metrics in our study include trees per hectare (TPH), canopy cover (%) per area (CC), average of five nearest neighbor distances (KNN), Clark and Evans Index (CEI), a canopy height-based solar radiation footprint (SRF), and an average distance from the northern canopy edge (NCE).
The UAV SfM point cloud data were collected across the entire study site with the same platform and sensor used for snow-series data collection and following completion of the mechanical thinning restoration treatment [39]. The final orthomosaic has error estimates of 1.14 m and 1.80 in the X, Y, and Z dimensions, respectively [39]. In addition, the SfM point cloud was used to create a digital terrain model (DTM) using the Cloth Simulation Filter (CSF) tool in CloudCompare v2.9.1 [39]. The CSF tools was parameterized using ‘relief’ for Scenes, a Cloth Resolution of 1 m, and a Classification Threshold of 0.7 m. The accuracy of the DTM (R2 = 0.95 and a RMSE = 2.98 m) was assessed by comparing points extracted from the DTM to corresponding points from differentially corrected Trimble GeoXH GPS elevation values collected from the trees in the field-based validation dataset.
Individual tree segmentation was performed on both the SfM and TLS point clouds using the Li et al. [65] algorithm in the lidR package [66] in RStudio (R-Studio Team, 2015). The segmented point clouds were then used to estimate each tree’s location (X, Y in UTM 12N m), maximum crown height (CH in m), and average of its widest and narrowest crown diameters (CD in m) using the rLiDAR package in R-Studio [67]. Each tree’s crown base height (CBH in m) was estimated from the point cloud data using a novel multiple changepoint detection algorithm [68]. This algorithm used the area of two-dimensional convex hulls, which were calculated at every 50 cm along a tree’s height to isolate the height at which trunk points transitioned into canopy points. Finally, each tree’s CV (m3) was then estimated by calculating the volume of the three-dimensional convex hull made from each tree’s points above its CBH value.

2.2.3. Forest Structure Metrics Validation

We assessed the accuracy of the SfM point-cloud-derived tree metrics by comparing them to both the field-measured and terrestrial laser scanner (TLS) point cloud-derived tree metrics (Figure 3). The TLS point cloud data provided accurate high-resolution estimates of all field-measured forest structure metrics and spatially extended the field-measured metrics while also alleviating known accuracy issues associated with UAV SfM-derived crown base height (CBH) and crown volume (CV) estimates [40]. Since we used a newer TLS instrument that has not been previously evaluated, an assessment of the omission and commission error rates was performed between the field-measured and TLS-derived datasets as well as between the TLS- and SfM-derived datasets. This facilitated direct comparison between trees from the TLS and SfM point clouds.
The TLS point cloud data was acquired within 0.64 ha field plots (n = 16) covering a total area of 5.17 ha and spread along a gradient of forest density present across the study site. The TLS data were collected using a Leica Geosystems BLK360 Imaging Laser Scanner (Leica Geosystems AG, 2020), which has a range of up to 60 m radius. The BLK360 captures 360,000 points per second at 830 nm wavelength with 300 and 360 degrees of vertical and horizontal field of view, respectively, and with 6mm-10m accuracies. Each TLS validation plot included at least three scans, which were tied together and georeferenced using four distinct reference targets. Reference targets were 50 × 50 cm reflective panels bolted to adjustable tripods at 1 m above ground, and their positions were recorded using a Trimble GeoXH GPS unit. We determined the locations and configuration of the scanner and reference targets using guidelines developed in our previous study [34]. Using Cloud Compare software, the scans were co-registered using the target center points, then merged into a single point cloud and georeferenced using the GPS coordinates of the target; the overall X, Y, and Z positional accuracy of the 16 point clouds is 0.46 m. The average point spacing of all point clouds was 0.03 m with a range of 0.025–0.046 m, and each was subsampled to 0.01 m to reduce redundant points and ensure consistent point spacing across scans.
The TLS point cloud data were validated using field-based measurements taken from 137 trees within 30 × 30 m validation plots (n = 16), covering a total of 1.44 ha. The field-based validation plots were centered within the TLS validation plots, providing as much overlap between the datasets as possible. Using the Trimble GeoXH handheld GPS unit with an attached Trimble laser range finder module, each tree’s geographic position (X,Y,Z coordinates), CH, CD, and CBH were measured and differentially corrected in GPS Pathfinder Office.

2.3. Data Analysis

2.3.1. Forest Structure Predictor Variables

The effect of forest canopy shading on persistent snow patches was quantified using distinct tree shading influence areas. Previous research shows that a tree’s ground shading influence extends between 1 and 2 times its crown height during the winter season in forests of the southwestern U.S. [55,69,70]. To support and refine this, we used three different spatial extents of localized tree shading to assess the size of persistent snow patches. These tree shading influence areas (TSIAs) were established with respect to each individual persistent snow patch and included trees located within a distance of 1, 1.5, and 2 times their crown heights. To ensure that only the trees capable of influencing a persistent snow patch were included, trees were selected based on whether their shadows were within solar azimuth angle extents specific to study site location and snow-series image date. The resulting minimum bounding extents were termed TSIA1.0, TSIA1.5, and TSIA2.0 (Figure 4). Once trees were selected in each TSIA, their vertical forest structure metrics were averaged to produce a final dataset.
In addition to the vertical forest structure metric summaries, each TSIA was assigned a set of horizontal arrangement metrics to quantify the spatial distribution of trees. Forest density was estimated by calculating trees per hectare (TPH) and average canopy cover (CC) in percent, while spacing of trees was measured using the average nearest neighbor (KNN) distance in meters from each tree to its five nearest neighbors. The clustering of trees was measured using a unitless Clark and Evans Index (CEI) value, with values < 1 indicating ordering and values > 1 indicating clustering. The amount of ground shading was expressed using an average solar radiation footprint (SRF) value in w/m2 for each TSIA polygon. The SRF was created using a sitewide canopy height model and the solar radiation toolset in ArcMap 10.8 to calculate the cumulative amount of incoming solar radiation (w/m2) during each snowstorm event [71,72]. The ground surface pixel values were averaged across the three storms and then across each TSIA polygon. Finally, each pixel’s distance to the northern canopy edge (DNCE) was calculated and averaged across each TSIA using a binary canopy cover raster based on the methodology developed by Mazzotti et al. [73]. The DNCE values provided a fine-scale metric that describes the directional within-stand differences of forest canopy spatial arrangement, and are well suited for assessing the effect of forest canopy radiative transfer.

2.3.2. Model Framework

An initial exploratory data analysis revealed sources of significant multicollinearity between forest structure metrics (independent variables) as well as numerous complex and non-linear relationships with persistent snow patch area (dependent variable). A variable selection process was used to eliminate the highly correlated and less descriptive forest structure metrics using variable inflation factor (VIF) scoring. To find a suitable model, we parameterized generalized linear models, random forest, support vector machine, and Multivariate Adaptive Regression Spline (MARS) algorithms for regression and tested their predictive accuracies using the Caret package [74] in RStudio. We selected the MARS model framework since it provided both the most accurate predictive accuracy as well as offering variable importance scores, which helped us interpret results into meaningful forest management recommendations.
The MARS model performs nonparametric multivariate regression without underlying assumptions of the data and is useful for regression problems with high-dimensional datasets and multiple predictor variable interactions. Our modeling framework assessed the relationship between persistent snow patch size (m2) and forest structure metrics for each of the three different tree shadow influence areas (TSIAs). Initial attempts at model fitting and prediction presented inconsistent results, with model hyperparameters and predictive accuracies fluctuating based on different randomly selected training and testing datasets, which can be apparent when using machine learning methods on smaller datasets [75]. Due to our relatively small sample size (n = 99), we adopted a split-sample validation procedure that was repeated 500 times. This provided 500 different randomly selected testing and training data partitions, model construction, and predictive accuracy assessments. While the number of repetition is usually arbitrarily chosen, we found that 500 repetitions provided a large enough sample size to clearly discriminate between trends in variable importance. Each model iteration utilized an exhaustive grid search to select the optimal set of hyperparameters, allowing for the inclusion of all possible predictors as well as for two-way interactions between predictors. An optimal model was then selected from each iteration by minimizing the generalized cross-validation error estimate, from which variable importance scores ranging from 0 to 100 were calculated for each model term.
Model results were summarized based on the frequency of both the hyperparameter value and the variable importance scores. Separately, the optimal model framework for each iteration was used for prediction on the testing dataset. Instead of arriving at a single ideal model that definitively explains the relationship between persistence snow patch area and forest structure, this approach provides a robust conceptual understanding of the relationship.

3. Results

3.1. Snow Cover Classification and Persistent Snow Patches

The snow cover classification accuracy assessment was performed using the third image from the second storm due to its mostly even distribution of snow/non-snow pixels. The snow classification performed well, with an overall accuracy of 90.2% and a kappa coefficient of 0.80. Relatively low rates of omission error were observed for the non-snow and snow classes at 13% and 6%, respectively, while similar rates of commission error were observed for the non-snow and snow classes at 6% and 14%, respectively (Table 2).
With the average daily high temperatures similar across all storm events (6.1 °C, 5.8 °C, and 5.7 °C for storms 1–3, respectively), storm 1 had the lowest initial snowfall amount (28.9 cm) and the lowest reduction in site-wide SCA (−51%) over 8 days. In contrast, storms 2 and 3 had greater initial snowfall amounts, 73.4 cm and 93.9 cm, respectively, and site-wide reductions in SCA were consistently higher and for longer time periods: −59% over 15 days and −70% over 26 days, respectively (Table 1).
The proportion of SCA in the treated portion of the study site on Day 1 was consistent across all storms, ranging from 91.2% to 89.4% to 90% for storms 1, 2, and 3, respectively. The proportion of SCA across the untreated portion of the study site on Day 1 of each storm was consistently lower, with values ranging from 43.4% to 62.3% to 58.1% for storms 1, 2, and 3, respectively. The effect of treatment condition on the magnitude of SCA reduction is evident both within and across storm events. The treated portion of the study site consistently exhibited a wider range of reductions in SCA within each storm compared to the narrower range in reductions observed in the untreated portion. The average of each storm’s total SCA reduction in the treated portion was −76.5% compared to −38.6% in the untreated portion (Figure 5).
Distinct persistent snow patches were identified from the final composite of classified snow-series images from grouped pixels that were covered in snow for 10 or 11 of the total days (n = 11) (Figure 6A). A total of 99 individual snow patches were delineated across the entire study site, covering 8.4% (6.36 ha) of the total area and ranging in size from 203 to 2699 m2 (SD = 469.7 m2), with an average size of 646.9 m2 and a majority (82%) being less than 1000 m2. The patches delineated across the treated portion of the study site (n = 67) had an average size of 722 m2 and covered 9.2% (4.84 ha) of the treated area (Figure 6B). In contrast, the patches in the untreated portion of the study site (n = 32) had an average size of 474 m2 and covered 6.6% (1.52 ha) of that area (Figure 6B).

3.2. Forest Structure Metrics Summary

The TLS-derived point cloud correctly identified 88% (n = 120) of trees from within the field-measured validation plots. The final TLS dataset consisted of 1293 trees identified within the TLS plot data extents (n = 16 validation plot, totaling 5.17 ha). The SfM-derived point cloud correctly identified 98% (n = 1280) trees within the TLS point cloud data extent, and the final dataset included a total of 8847 trees across the 76-ha study site. The individual tree metrics compared between all datasets are identical except for CV, which was calculated from TLS and SfM point clouds (Figure 2).
The accuracy of the individual tree structure metrics used in the final model are a product of a multi-scale accuracy assessment, which compared the field-measured, TLS-derived, and SfM-derived estimates. The relationships between the field-measured and TLS-derived tree metrics were generally strong, while the relationships between the TLS- and SfM-derived tree metrics varied.
The final set of forest structure predictors used for modeling were selected to balance both accurate tree dimension estimation as well as to maximize the model’s predictive capacity. The vertical forest structure metrics were selected by first examining the relationships between SfM- and TLS-derived estimates, then by assessing their multicollinearity and variance inflation factor (VIF) scores. For example, the weak relationship between SfM- and TLS-derived CD (R2 = 0.14; RMSE = 1.62 m) indicated that it should not be used as a predictor variable. In addition, CD was highly correlated with CV (R2 = 0.72) and had a high VIF score (12.3), indicating that it would negatively impact model accuracy and interpretability. For the entire study site, the complete set of vertical and horizontal forest structure metrics included as model predictors are summarized by TSIA size (1, 1.5, and 2).

3.3. Relationship between Snow and Forest Structure

The most frequent number of model terms for TSIA 1.0 was 5 (50% of all model iterations), for TSIA 1.5 were 6 (35% of all model iterations), and for TSIA 2.0 were 7 (35% of all model iterations), while the most frequent interaction degree term was 1 for all TSIA groups. The average model error estimates calculated across all model iterations grouped by TSIA size were 0.14 for TSIA 1.0, 0.11 for TSIA 1.5, and 0.14 for TSIA 2.0. As the influence region size increased, more predictor variables were needed to stabilize the model and there was less consensus in the number of predictors selected.
For all TSIA sizes, the forest structure metrics with the highest variable importance scores were the mean canopy cover (CC), mean solar radiation footprint (SRF), and trees per hectare (TPH) values (Figure 7). While this remained relatively consistent overall, differences in variable importance scores existed across TSIA sizes as evidenced by the differences between variable importance scores from TSIA groups within CEI, mean SRF, and mean TPH metrics. Specifically, the greater importance were observed for CEI at the largest spatial scale and the lesser importance of mean SRF and mean TPH at the smallest spatial scale.
The models describing forest structure at TSIA 1.5 performed best overall with a mean prediction accuracy R2 = 0.70 (RMSE = 267 m2) (Figure 8), followed by those at TSIA 1.0 and TSIA 2.0 with mean prediction accuracies of R2 = 0.66 (RMSE = 286 m2) and R2 = 0.61 (RMSE = 307 m2), respectively. Of note is the increasing variability in model predictive accuracy on both ends of the persistent snow patch area, especially on the larger end where wide fanning indicates greater uncertainty.

4. Discussion

Seasonal snowpack provides vital water resources for both human- and ecosystem-oriented services in semi-arid environments. Quantifying the spatially heterogenous and often ephemeral nature of this snow at the mid- (10–1000 ha) and landscape-scales (400+ ha) requires accurate spatially extensive and temporally dense datasets [37]. The structure and arrangement of trees in dry forests are central components influencing snow cover and persistent snowpack. In this study, we demonstrate an accurate method to provide near real-time estimates of mid-scale snow covered area (SCA) and assess how forest structure influences persistent snow patch size in a thinned ponderosa pine forest on the southern edge of the North American continental snow distribution. We found that forest structure metrics emphasizing the spatial arrangement of trees and tree groups were more influential on persistent snow cover (Figure 7), and that these effects were most pronounced when considering trees within 1.5 times their height to persistent snow cover (Figure 8).
We first quantified snow-covered area (SCA) across the discontinuous forest of our 76-ha study site using UAV snow-series datasets. In forested environments, SCA is often underestimated in remote sensing data propagating from lower resolution satellite and airborne imagery, with trees masking the ground surface and tree shadows being misclassified as ‘non-snow’ [76,77,78,79]. The overall accuracy of our high-resolution binary snow classification was 90.2%, indicating strategic UAV data acquisition, associated high-resolution imagery, and a relatively straightforward classification process can be used to quantify SCA in forested areas [80]. Furthermore, by capturing inter-storm reductions in SCA (Figure 5), we demonstrate this approach can effectively delineate regions harboring the most persistent SCA and do so at a fine spatial resolution (Figure 6). This level of detail and accuracy provided the foundation for our assessment of forest structure impact on persistent snow cover patches.
Except for crown diameter (CD), we found good overall agreement between the UAV SfM estimates and our validation datasets. Horizontal forest structure metrics were more influential than vertical structure metrics in predicting persistent snow patch size (Figure 7). Specifically, tree canopy cover (CC), trees per ha (TPH), and solar radiation footprint (SRF) were the most explanatory variables. This result supports our hypotheses that snow cover and, therefore, subcanopy shortwave radiation, is moderated by the spatial arrangement of trees more than by the vertical structure of individual trees.
While the horizontal metrics of CC, TPH, and SRF were the most influential predictor variables, they were also the metrics with the greatest variability when grouped by persistent snow patch size. These are metrics of foliar cover and tree stem density that directly influence subcanopy shading, which subsequently impacts the relatively shallow, intermittent, and spatially heterogeneous snow cover found at our study site. In addition to the horizontal metrics, we expected the mean CEI and mean KNN distance values to be influential because they describe the location and spacing of trees. However, their lack of influence might indicate a need to adjust the scale and scope at which these metrics are calculated. For example, the spacing of groups of trees might be a more valuable metric than the spacing of individual trees when considering sub-canopy ground shading.
In contrast to the horizontal metrics, variables in the vertical metric group were only slightly or moderately useful predictors of persistent snow patch size (Figure 7). We anticipated that forest patches comprised of trees with noticeably different vertical structures would result in differences in the forest canopy shading and ultimately statistically significant impacts on persistent snow patch size. This is illustrated in Figure 9, where the amount of shading seemingly provided by trees with higher CBH:CH ratio and lower CV values (Panel A) is noticeably lower compared to the amount of shading provided by trees with overall larger crowns (Panel B). The post-treatment forest structure conditions in Figure 9A are more common than those in Figure 9B. A potential cause for this could be a combination of the relatively homogenous vertical structure metrics observed across the study site (Figure 1B) and the relatively low agreement (R2 = 0.50; RMSE = 3.23 m) between the TLS and field-measured CBH values.
Mean canopy height (CH) had strong agreement between both the field-measured and TLS-estimated values and the TLS- and SfM-estimated values, as well as low variability when grouped by persistent snow patch size. Accurate tree height estimates are not surprising given the ability for SfM and TLS to accurately estimate tree height regardless of forest density [34,39,41] and the relatively homogenous post-settlement forest structure conditions present across the study site. Given the low variable importance scores, mean CH appeared unimportant when modeling persistent snow patch area. However, we believe the importance of mean CH is evidenced in the superior predictive ability of the TSIA 1.5 category. More specifically, using tree CH to define the north–south width of forest canopy gaps (in the northern hemisphere), specifically to 1.5 times the average crown height of adjacent trees, provides the most consistent estimates of both snow cover and snow persistence. Early research shows that snow persisted in ’zones of retention’ corresponding to 1.5 to 2 times the height of adjacent tree stands in ponderosa pine forests [81]. While these early studies assumed that persistent snow retention zones were driven by forest structure, mainly tree CH, we use high-resolution UAV-based estimates of CH to confirm the 1.5 CH value as well as to identify horizontal metrics contributing to persistent snow. Our results indicate that mean tree CH and tree patch spacing at 1.5 times CH are important considerations in forest restoration, in addition to the horizontal metrics described above, for maximizing snow cover retention and water yield.
As thinning-based restoration is expanded throughout dry forests that harbor seasonal snowpack, it is crucial to understand how restoration-driven canopy reduction will impact forest water balance [82] and snowmelt water inputs into shallow groundwater reservoirs. Previous research shows that snow accumulation in discontinuous or disturbed forests can be greater in less dense forests and within large canopy openings [37,73]. The restoration thinning at our study site significantly reduced both canopy cover (from about 40% to 10%) and tree density (from about 212 TPH to 65 TPH), while increasing the number of forest patches (from 39 to 133) and decreasing the mean forest patch area (from 0.68 ha to 0.13 ha) [39]. While this study focused on developing remote sensing techniques to quantify snowpack dynamics, contrasting the treated versus the untreated portions of our study area suggest important treatment impacts that should be confirmed by more replicated comparisons. Specifically, our results show there was more overall persistent snow cover in the thinned portion of the study site compared to the untreated portion (10% and 7%, respectively). Our results generally support and provide further insight into restoration-based reductions in forest density to promote snow cover. For example, the largest persistent snow patches were located adjacent to forest patches with 31–33% canopy cover (CC), while the smallest were located adjacent to patches with 50–52% CC, a trend which is supported by the optimal snow persistence CC value of 24% provided by regional satellite-derived estimates [39].
Our findings suggest that tailoring future forest restoration treatments to promote persistent snow cover in southwestern U.S. dry forest ecosystems should continue to focus on horizontal forest structure metrics like forest density and canopy cover at the landscape-scale (>1000 ha). Our results also underscore and refine the importance of horizontal forest structure metrics like CC, TPH, and SRF, regardless of spatial scale. However, we propose these criteria be emphasized during the creation of forest patches at relatively fine patch scale (<4 ha). It is critical to regional restoration efforts to continue implementing the commonly accepted standards in southwestern dry forest restoration, like promoting variation in interspace size, tree group size and within-group tree spacing [80]. In addition, we propose landscape-scale restoration treatment goals should also be implemented at fine spatial scale, operating within and among individual forest patches. For example, assuming a restoration treatment reflects diversity in forest patch size, spacing, and density, an overarching goal may be to achieve 24–33% CC across a 1000-ha treatment area. While this level of CC could be measured across the entire treatment area, having individual forest patches at 24–33% CC and distributed with distances at 1.5 times the average tree height within the patch should also promote localized persistent snow cover.

5. Conclusions

This study considered the inclusion of detailed forest structure metrics in quantifying persistent snow cover in a dry southwestern U.S. forest. We found the size (m2) of persistent snow patches can be effectively predicted using targeted forest structure metrics. Specifically, the most effective predictor metrics included tree shadows that are 1.5 times the tree heights, as well as tree density and canopy cover within this shaded area. While our findings underscore the importance of forest canopy shading on persistent snow cover, they also indicate the relationship between persistent snow cover and fine-scale forest structure is likely more complex, rooted in different variables, or present at different spatial scales. Maximizing persistent snow cover in dry forest environments can be achieved by controlling subcanopy shading at the ground surface through an optimal set of fine-scale forest structure and spacing metrics. Future research and restoration efforts can achieve this by coupling our UAV-based methodology for quantifying persistent snow cover with more descriptive measurements of snow dynamics, such as snow depth and snow water equivalent.
Our results support the utility of thinning-based forest restoration in dry southwestern forests to promote snow cover retention and forest health. We show there is a wide range of persistent snow patch sizes across thinned forest, and that differences in fine-scale forest structure are important for maximizing snow persistence. Adjusting existing restoration thinning prescriptions to reflect landscape-scale goals in fine-scale forest patches will help further this objective while promoting broader ecosystem health and resiliency.

Author Contributions

Conceptualization: T.S. and A.B.; Formal analysis: A.B.; Funding acquisition: T.S.; Methodology: T.S. and A.B.; Resources: T.S.; Supervision: T.S.; Validation: A.B.; Writing—original draft: A.B.; Writing—review and editing: T.S., J.B. (John Bradford), J.B. (Joel Biderman), S.G., and T.K. All authors have read and agreed to the published version of the manuscript.

Funding

The research was funded by the Salt River Project (Grant Number: 8200007407) and the Nature Conservancy (Grant Number: AZFO-170712). The UAV equipment used in this study was funded by Northern Arizona University Office of Vice President for Research.

Acknowledgments

Any use of trade, firm, or product name is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflicts of Interest

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

References

  1. Barnett, T.P.; Adam, J.C.; Lettenmaier, D.P. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 2005, 438, 303–309. [Google Scholar] [CrossRef]
  2. Bales, R.C.; Hopmans, J.W.; O’Geen, A.T.; Meadows, M.; Hartsough, P.C.; Kirchner, P.; Hunsaker, C.T.; Beaudette, D. Soil Moisture Response to Snowmelt and Rainfall in a Sierra Nevada Mixed-Conifer Forest. Vadose Zone J. 2011, 10, 786–799. [Google Scholar] [CrossRef]
  3. French, H.; Binley, A.M. Snowmelt infiltration: Monitoring temporal and spatial variability using time-lapse electrical resistivity. J. Hydrol. 2004, 297, 174–186. [Google Scholar] [CrossRef]
  4. Newman, B.D.; Campbell, A.R.; Wilcox, B.P. Lateral subsurface flow pathways in a semiarid Ponderosa pine hillslope. Water Resour. Res. 1998, 34, 3485–3496. [Google Scholar] [CrossRef]
  5. Price, A.G.; Hendrie, L.K. Water motion in a deciduous forest during snowmelt. J. Hydrol. 1983, 64, 339–356. [Google Scholar] [CrossRef]
  6. Wilcox, B.P.; Newman, B.D.; Brandes, D.; Davenport, D.W.; Reid, K. Runoff from a semiarid Ponderosa pine hillslope in New Mexico. Water Resour. Res. 1997, 33, 2301–2314. [Google Scholar] [CrossRef]
  7. Harpold, A.A.; Brooks, P.D. Humidity determines snowpack ablation under a warming climate. Proc. Natl. Acad. Sci. USA 2018, 115, 1215–1220. [Google Scholar] [CrossRef] [Green Version]
  8. Barnett, T.P.; Pierce, D.W.; Hidalgo, H.G.; Bonfils, C.; Santer, B.D.; Das, T.; Bala, G.; Wood, A.W.; Nozawa, T.; Mirin, A.A.; et al. Human-Induced Changes in the Hydrology of the Western United States. Science 2008, 319, 1080–1083. [Google Scholar] [CrossRef] [Green Version]
  9. Hidalgo, H.G.; Das, T.; Dettinger, M.D.; Cayan, D.R.; Pierce, D.W.; Barnett, T.P.; Bala, G.; Mirin, A.A.; Wood, A.W.; Bonfils, C.; et al. Detection and Attribution of Streamflow Timing Changes to Climate Change in the Western United States. J. Clim. 2009, 22, 3838–3855. [Google Scholar] [CrossRef]
  10. Safeeq, M.; Shukla, S.; Arismendi, I.; Grant, G.E.; Lewis, S.L.; Nolin, A.W. Influence of winter season climate variability on snow-precipitation ratio in the western United States. Int. J. Clim. 2016, 36, 3175–3190. [Google Scholar] [CrossRef]
  11. Allen, C.D.; Breshears, D.D.; McDowell, N.G. On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 2015, 6, 1–55. [Google Scholar] [CrossRef]
  12. Hereford, R. Climate Variation at Flagstaff, Arizona’s Mountain Town Compiled from National Weather Service Data; NOAA Flagstaff Weather Station: Flagstaff, AZ, USA, 2014.
  13. Ehleringer, J.R.; Phillips, S.L.; Schuster, W.S.F.; Sandquist, D.R. Differential utilization of summer rains by desert plants. Oecologia 1991, 88, 430–434. [Google Scholar] [CrossRef]
  14. Kerhoulas, L.P.; Kolb, T.E.; Koch, G.W. Tree size, stand density, and the source of water used across seasons by ponderosa pine in northern Arizona. For. Ecol. Manag. 2013, 289, 425–433. [Google Scholar] [CrossRef]
  15. Essery, R.; Pomeroy, J.; Parviainen, J.; Storck, P. Sublimation of Snow from Coniferous Forests in a Climate Model. J. Clim. 2003, 16, 1855–1864. [Google Scholar] [CrossRef] [Green Version]
  16. Harpold, A.A.; Biederman, J.A.; Condon, K.; Merino, M.; Korgaonkar, Y.; Nan, T.; Sloat, L.L.; Ross, M.; Brooks, P.D. Changes in snow accumulation and ablation following the Las Conchas Forest Fire, New Mexico, USA. Ecohydrology 2014, 7, 440–452. [Google Scholar] [CrossRef]
  17. Molotch, N.P.; Blanken, P.D.; Williams, M.W.; Turnipseed, A.A.; Monson, R.K.; Margulis, S.A. Estimating sublimation of intercepted and sub-canopy snow using eddy covariance systems. Hydrol. Process. 2007, 21, 1567–1575. [Google Scholar] [CrossRef]
  18. Molotch, N.P.; Brooks, P.D.; Burns, S.P.; Litvak, M.; Monson, R.K.; McConnell, J.R.; Musselman, K. Ecohydrological controls on snowmelt partitioning in mixed-conifer sub-alpine forests. Ecohydrology 2009, 2, 129–142. [Google Scholar] [CrossRef]
  19. Roth, T.R.; Nolin, A.W. Forest impacts on snow accumulation and ablation across an elevation gradient in a temperate montane environment. Hydrol. Earth Syst. Sci. 2017, 21, 5427–5442. [Google Scholar] [CrossRef] [Green Version]
  20. Varhola, A.; Coops, N.C.; Bater, C.W.; Teti, P.; Boon, S.; Weiler, M. The influence of ground- and lidar-derived forest structure metrics on snow accumulation and ablation in disturbed forests. Can. J. For. Res. 2010, 40, 812–821. [Google Scholar] [CrossRef]
  21. Davis, R.E.; Hardy, J.P.; Ni, W.; Woodcock, C.; McKenzie, J.C.; Jordan, R.; Li, X. Variation of snow cover ablation in the boreal forest: A sensitivity study on the effects of conifer canopy. J. Geophys. Res. Space Phys. 1997, 102, 29389–29395. [Google Scholar] [CrossRef]
  22. Dickerson-Lange, S.E.; Lutz, J.A.; Martin, K.A.; Raleigh, M.S.; Gersonde, R.; Lundquist, J.D. Evaluating observational methods to quantify snow duration under diverse forest canopies. Water Resour. Res. 2015, 51, 1203–1224. [Google Scholar] [CrossRef]
  23. Essery, R.; Bunting, P.; Rowlands, A.; Rutter, N.; Hardy, J.; Melloh, R.; Link, T.; Marks, D.; Pomeroy, J. Radiative Transfer Modeling of a Coniferous Canopy Characterized by Airborne Remote Sensing. J. Hydrometeorol. 2008, 9, 228–241. [Google Scholar] [CrossRef] [Green Version]
  24. Lawler, R.R.; Link, T.E. Quantification of incoming all-wave radiation in discontinuous forest canopies with application to snowmelt prediction. Hydrol. Process. 2011, 25, 3322–3331. [Google Scholar] [CrossRef]
  25. Veatch, W.; Brooks, P.D.; Gustafson, J.R.; Molotch, N.P. ‘Quantifying the effects of forest canopy cover on net snow accumulation at a continental, mid-latitude site’. Ecohydrology 2009, 2, 115–128. [Google Scholar] [CrossRef]
  26. Bauwens, S.; Bartholomeus, H.; Calders, K.; Lejeune, P. Forest Inventory with Terrestrial LiDAR: A Comparison of Static and Hand-Held Mobile Laser Scanning. Forests 2016, 7, 127. [Google Scholar] [CrossRef] [Green Version]
  27. Calders, K.; Origo, N.; Burt, A.; Disney, M.; Nightingale, J.; Raumonen, P.; Åkerblom, M.; Malhi, Y.; Lewis, P. Realistic Forest Stand Reconstruction from Terrestrial LiDAR for Radiative Transfer Modelling. Remote Sens. 2018, 10, 933. [Google Scholar] [CrossRef] [Green Version]
  28. Dassot, M.; Constant, T.; Fournier, M. The use of terrestrial LiDAR technology in forest science: Application fields, benefits and challenges. Ann. For. Sci. 2011, 68, 959–974. [Google Scholar] [CrossRef] [Green Version]
  29. Sankey, T. Decadal-scale aspen change detection using Landsat 5TM and lidar data. Appl. Veg. Sci. 2012, 15, 84–98. [Google Scholar]
  30. Liang, X.; Hyyppä, J.; Kaartinen, H.; Pyörälä, J.; Lehtomäki, M.; Holopainen, M.; Kankare, V.; Luoma, V.; Saarinen, N.; Chen, L.; et al. International benchmarking of terrestrial laser scanning approaches for forest inventories; joint project of eurosdr and isprs; part i: Objective, datasets, evaluation criteria and methods. EuroSDR 2019, 2019, 1–29. [Google Scholar] [CrossRef]
  31. Næsset, E. Practical large-scale forest stand inventory using a small-footprint airborne scanning laser. Scand. J. For. Res. 2004, 19, 164–179. [Google Scholar] [CrossRef]
  32. Reutebuch, S.E.; Andersen, H.-E.; McGaughey, R.J. Light Detection and Ranging (LIDAR): An Emerging Tool for Multiple Resource Inventory. J. For. 2005, 103, 286–292. Available online: https://academic.oup.com/jof/article-abstract/103/6/286/4598654 (accessed on 30 March 2019).
  33. Sankey, T.; Shrestha, R.; Sankey, J.B.; Hardegree, S.; Strand, E. Lidar-derived estimate and uncertainty of carbon sink in successional phases of woody encroachment. J. Geophys. Res. Biogeosci. 2013, 118, 1144–1155. [Google Scholar] [CrossRef]
  34. Donager, J.J.; Sankey, T.T.; Sankey, J.B.; Sanchez Meador, A.J.; Springer, A.E.; Bailey, J.D. Examining forest structure with terrestrial lidar: Suggestions and novel techniques based on comparisons between scanners and forest treat-ments. Earth Space Sci. 2018, 5, 753–776. [Google Scholar] [CrossRef] [Green Version]
  35. Kongoli, C.; Romanov, P.; Ferraro, R. Snow cover monitoring from remote-sensing satellites: Possibilities for drought assessment. In Remote Sensing of Drought: Innovative Monitoring Approaches; CRC Press: Boca Raton, FL, USA, 2012; pp. 359–386. [Google Scholar] [CrossRef]
  36. Nolin, A.W. Recent advances in remote sensing of seasonal snow. J. Glaciol. 2010, 56, 1141–1150. [Google Scholar] [CrossRef] [Green Version]
  37. Sankey, T.; Donald, J.; McVay, J.; Ashley, M.; O’Donnell, F.; Lopez, S.M.; Springer, A. Multi-scale analysis of snow dynamics at the southern margin of the North American continental snow distribution. Remote Sens. Environ. 2015, 169, 307–319. [Google Scholar] [CrossRef]
  38. Alonzo, M.; Andersen, H.-E.; Morton, D.C.; Cook, B.D. Quantifying Boreal Forest Structure and Composition Using UAV Structure from Motion. Forests 2018, 9, 119. [Google Scholar] [CrossRef] [Green Version]
  39. Belmonte, A.; Sankey, T.; Biederman, J.A.; Bradford, J.; Goetz, S.J.; Kolb, T.; Woolley, T. UAV -derived estimates of forest structure to inform ponderosa pine forest restoration. Remote Sens. Ecol. Conserv. 2019, 6, 181–197. [Google Scholar] [CrossRef]
  40. Shin, P.; Sankey, T.; Moore, M.M.; Thode, A.E. Evaluating Unmanned Aerial Vehicle Images for Estimating Forest Canopy Fuels in a Ponderosa Pine Stand. Remote Sens. 2018, 10, 1266. [Google Scholar] [CrossRef] [Green Version]
  41. Sankey, T.; Donager, J.; McVay, J.; Sankey, J.B. UAV lidar and hyperspectral fusion for forest monitoring in the southwestern USA. Remote Sens. Environ. 2017, 195, 30–43. [Google Scholar] [CrossRef]
  42. Sankey, T.T.; Leonard, J.M.; Moore, M.M. Unmanned Aerial Vehicle-Based Rangeland Monitoring: Examining a Century of Vegetation Changes. Rangel. Ecol. Manag. 2019, 72, 858–863. [Google Scholar] [CrossRef]
  43. Sankey, J.; Sankey, T.; Li, J.; Ravi, S.; Wang, G.; Caster, J.; Kasprak, A. Quantifying plant-soil-nutrient dynamics in rangelands: Fusion of UAV hyperspectral-LiDAR, UAV multispectral-photogrammetry, and ground-based LiDAR-digital photography in a shrub-encroached desert grassland. Remote Sens. Environ. 2021, 253, 112223. [Google Scholar] [CrossRef]
  44. Guerra-Hernandez, J.; Gonzalez-Ferreiro, E.; Sarmento, A.; Silva, J.; Nunes, A.; Correia, A.C.; Fontes, L.; Tomé, M.; Diaz-Varela, R. Using high resolution UAV imagery to estimate tree variables in Pinus pinea plantation in Portugal. For. Syst. 2016, 25, eSC09. [Google Scholar] [CrossRef] [Green Version]
  45. Iizuka, K.; Yonehara, T.; Itoh, M.; Kosugi, Y. Estimating Tree Height and Diameter at Breast Height (DBH) from Digital Surface Models and Orthophotos Obtained with an Unmanned Aerial System for a Japanese Cypress (Chamaecyparis obtusa) Forest. Remote Sens. 2017, 10, 13. [Google Scholar] [CrossRef] [Green Version]
  46. Puliti, S.; Ørka, H.O.; Gobakken, T.; Næsset, E. Inventory of Small Forest Areas Using an Unmanned Aerial System. Remote Sens. 2015, 7, 9632–9654. [Google Scholar] [CrossRef] [Green Version]
  47. White, J.C.; Wulder, M.A.; Vastaranta, M.; Coops, N.C.; Pitt, D.; Woods, M. The Utility of Image-Based Point Clouds for Forest Inventory: A Comparison with Airborne Laser Scanning. Forests 2013, 4, 518–536. [Google Scholar] [CrossRef] [Green Version]
  48. Carr, J.C.; Slyder, J.B. Individual tree segmentation from a leaf-off photogrammetric point cloud. Int. J. Remote Sens. 2018, 39, 5195–5210. [Google Scholar] [CrossRef]
  49. Thiel, C.; Schmullius, C. Comparison of UAV photograph-based and airborne lidar-based point clouds over forest from a forestry application perspective. Int. J. Remote Sens. 2016, 38, 2411–2426. [Google Scholar] [CrossRef]
  50. Wallace, L.; Lucieer, A.; Malenovský, Z.; Turner, D.; Vopěnka, P. Assessment of Forest Structure Using Two UAV Techniques: A Comparison of Airborne Laser Scanning and Structure from Motion (SfM) Point Clouds. Forests 2016, 7, 62. [Google Scholar] [CrossRef] [Green Version]
  51. Harder, P.; Schirmer, M.; Pomeroy, J.; Helgason, W. Accuracy of snow depth estimation in mountain and prairie environments by an unmanned aerial vehicle. Cryosphere 2016, 10, 2559–2571. [Google Scholar] [CrossRef] [Green Version]
  52. Lendzioch, T.; Langhammer, J.; Jenicek, M. Estimating Snow Depth and Leaf Area Index Based on UAV Digital Photogrammetry. Sensors 2019, 19, 1027. [Google Scholar] [CrossRef] [Green Version]
  53. Miziński, B.; Niedzielski, T. Fully-automated estimation of snow depth in near real time with the use of unmanned aerial vehicles without utilizing ground control points. Cold Reg. Sci. Technol. 2017, 138, 63–72. [Google Scholar] [CrossRef]
  54. Niedzielski, T.; Spallek, W.; Witek-Kasprzak, M. Automated Snow Extent Mapping Based on Orthophoto Images from Unmanned Aerial Vehicles. Pure Appl. Geophys. 2018, 175, 3285–3302. [Google Scholar] [CrossRef] [Green Version]
  55. Baker, M.B.; Ffolliott, P.F. Role of Snow Hydrology in Watershed Management. J. Ariz. Nev. Acad. Sci. 2003, 35, 42–47. [Google Scholar] [CrossRef]
  56. Sackett, S.S. Reducing Natural Ponderosa Pine Fuels Using Prescribed Fire: Two Case Studies; USDA Forest Service, Rocky Mountain Forest and Range Experiment Station: Fort Collins, CO, USA, 1980.
  57. Allen, C.D.; Savage, M.; Falk, D.A.; Suckling, K.F.; Swetnam, T.W.; Schulke, T.; Stacey, P.B.; Morgan, P.; Hoffman, M.; Klingel, J.T. Ecological Restoration of Southwestern Ponderosa Pine Ecosystems: A Broad Perspective. Ecol. Appl. 2002, 12, 1418–1433. [Google Scholar] [CrossRef]
  58. Larson, A.J.; Churchill, D. Tree spatial patterns in fire-frequent forests of western North America, including mechanisms of pattern formation and implications for designing fuel reduction and restoration treatments. For. Ecol. Manag. 2012, 267, 74–92. [Google Scholar] [CrossRef]
  59. Reynolds, R.T.; Meador, A.J.S.; Youtz, J.A.; Nicolet, T.; Matonis, M.S.; Jackson, P.L.; DeLorenzo, D.G.; Graves, A.D. Restoring Composition and Structure in Southwestern Frequent-Fire Forests: A Science-Based Framework for Improving Ecosystem Resiliency; USDA Forest Service: Fort Collins, CO, USA, 2013.
  60. Massey, R.; Sankey, T.T.; Yadav, K.; Congalton, R.G.; Tilton, J.C. Integrating cloud-based workflows in continental-scale cropland extent classification. Remote Sens. Environ. 2018, 219, 162–179. [Google Scholar] [CrossRef]
  61. Sankey, T.T.; Massey, R.; Yadav, K.; Congalton, R.G.; Tilton, J.C. Post-socialist cropland changes and abandonment in Mongolia. Land Degrad. Dev. 2018, 29, 2808–2821. [Google Scholar] [CrossRef]
  62. Rouse, J.; Haas, R.; Schell, J.; Deering, D. Monitoring Vegetation Systems in the Great Plains with ERTS. NASA Spec. Publ. 1974, 351, 309. [Google Scholar]
  63. Huete, A. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 1988, 25, 295–309. [Google Scholar] [CrossRef]
  64. Haralick, R.M.; Shanmugam, K.; Dinstein, I. Textural Features for Image Classification. IEEE Trans. Syst. Man, Cybern. 1973, 3, 610–621. [Google Scholar] [CrossRef] [Green Version]
  65. Li, W.; Guo, Q.; Jakubowski, M.K.; Kelly, M. A New Method for Segmenting Individual Trees from the Lidar Point Cloud. Photogramm. Eng. Remote Sens. 2012, 78, 75–84. [Google Scholar] [CrossRef] [Green Version]
  66. Roussel, J.-R.; Caspersen, J.; Béland, M.; Thomas, S.; Achim, A. Removing bias from LiDAR-based estimates of canopy height: Accounting for the effects of pulse density and footprint size. Remote Sens. Environ. 2017, 198, 1–16. [Google Scholar] [CrossRef]
  67. Mohan, M.; Silva, C.A.; Klauberg, C.; Jat, P.; Catts, G.; Cardil, A.; Hudak, A.T.; Dia, M. Individual Tree Detection from Unmanned Aerial Vehicle (UAV) Derived Canopy Height Model in an Open Canopy Mixed Conifer Forest. Forests 2017, 8, 340. [Google Scholar] [CrossRef] [Green Version]
  68. Runge, V.; Hocking, T.D.; Romano, G.; Afghah, F.; Fearnhead, P.; Rigaill, G. Gfpop: An R Package for Univariate Graph-Constrained Change-point Detection. arXiv 2020, arXiv:2002.03646. [Google Scholar]
  69. Ffolliott, P.F. Snowpack dynamics in mountain areas: Research findings in the southwestern United States. In Proceedings of the International Symposium of Mountainous Areas, Shimla, India, 28–30 May 1992; pp. 129–139. [Google Scholar]
  70. Ffolliott, P.F.; Thorud, D.B. A Technique to Evaluate Snowpack Profiles in and Adjacent to Forest Openings. Hydrology and Water Resources in Arizona and the Southwest. 1974. Available online: http://hdl.handle.net/10150/300274 (accessed on 5 March 2021).
  71. Olpenda, A.S.; Stereńczak, K.; Będkowski, K. Modeling Solar Radiation in the Forest Using Remote Sensing Data: A Review of Approaches and Opportunities. Remote Sens. 2018, 10, 694. [Google Scholar] [CrossRef] [Green Version]
  72. Abdollahnejad, A.; Panagiotidis, D.; Surový, P.; Ulbrichová, I. UAV Capability to Detect and Interpret Solar Radiation as a Potential Replacement Method to Hemispherical Photography. Remote Sens. 2018, 10, 423. [Google Scholar] [CrossRef] [Green Version]
  73. Mazzotti, G.; Currier, W.R.; Deems, J.S.; Pflug, J.M.; Lundquist, J.D.; Jonas, T. Revisiting Snow Cover Variability and Canopy Structure Within Forest Stands: Insights from Airborne Lidar Data. Water Resour. Res. 2019, 55, 6198–6216. [Google Scholar] [CrossRef]
  74. Kuhn, M. Building Predictive Models in R Using the caret Package. J. Stat. Softw. 2008, 28, 1–26. [Google Scholar] [CrossRef] [Green Version]
  75. Breiman. Random Forests. Mach. Learn. 2001, 45, 5–32. [Google Scholar] [CrossRef] [Green Version]
  76. Hall, D.; Foster, J.; Salomonson, V.; Klein, A.; Chien, J. Development of a technique to assess snow-cover mapping errors from space. IEEE Trans. Geosci. Remote Sens. 2001, 39, 432–438. [Google Scholar] [CrossRef]
  77. Huang, X.; Deng, J.; Wang, W.; Feng, Q.; Liang, T. Impact of climate and elevation on snow cover using integrated remote sensing snow products in Tibetan Plateau. Remote Sens. Environ. 2017, 190, 274–288. [Google Scholar] [CrossRef]
  78. Metsämäki, S.; Vepsäläinen, J.; Pulliainen, J.; Sucksdorff, Y. Improved linear interpolation method for the estimation of snow-covered area from optical data. Remote Sens. Environ. 2002, 82, 64–78. [Google Scholar] [CrossRef]
  79. Vikhamar, D.; Solberg, R. Snow-cover mapping in forests by constrained linear spectral unmixing of MODIS data. Remote Sens. Environ. 2003, 88, 309–323. [Google Scholar] [CrossRef]
  80. Eker, R.; Bühler, Y.; Schlögl, S.; Stoffel, A.; Aydın, A. Monitoring of Snow Cover Ablation Using Very High Spatial Resolution Remote Sensing Datasets. Remote Sens. 2019, 11, 699. [Google Scholar] [CrossRef] [Green Version]
  81. Ffolliott, P.F.; Hansen, E.A.; Zander’, A.D. Snow in Natural Openings and Adjacent Ponderosa Pine Stands on the Beaver Creek Watersheds; Rocky Mountain Forest and Range Experiment Station: Fort Collins, CO, USA, 1965. Available online: https://books.googleusercontent.com/books/content?req=AKW5Qac6Sx9t2jj2Y_uiea4q0fNgSERB-rKaskNTOTxw0nyENCwk1wNndCEQ3YuiHGYUfN3hsThrV2BZSkuNNpbJ1WnIMeVt5wn8rbgHx_JvhU5qZ9mPxEahSsLyDbZc1x3xZD2sa4nwHfVn-_r85Ddhnxb6Lu_LFExMP18qrLFfaB6LOZrU3CETk6pb2R1r0NaD7S0Cc (accessed on 15 April 2019).
  82. Sankey, T.; Belmonte, A.; Massey, R.; Leonard, J. Regional-scale forest restoration effects on ecosystem resiliency to drought: A synthesis of vegetation and moisture trends on Google Earth Engine. Remote Sens. Ecol. Conserv. 2020. [Google Scholar] [CrossRef]
Remotesensing 13 01036 g001 550
Figure 1. An overview of the study site highlighting the locations and extents of the UAV Structure-from-Motion (SfM) data, terrestrial laser scanner (TLS) data, and field-based validation data. (A) shows an example of the field-measured and TLS validated plots. (B) shows the distribution of all field-measured and TLS validated plots within the SfM data extent, which includes both thinned and unthinned portions of the forest. A majority of our ground-based measurements are distributed across the larger, thinned portion of the study site. (C) shows the UAV and TLS instruments used in this study.
Figure 1. An overview of the study site highlighting the locations and extents of the UAV Structure-from-Motion (SfM) data, terrestrial laser scanner (TLS) data, and field-based validation data. (A) shows an example of the field-measured and TLS validated plots. (B) shows the distribution of all field-measured and TLS validated plots within the SfM data extent, which includes both thinned and unthinned portions of the forest. A majority of our ground-based measurements are distributed across the larger, thinned portion of the study site. (C) shows the UAV and TLS instruments used in this study.
Remotesensing 13 01036 g001
Remotesensing 13 01036 g002 550
Figure 2. UAV multispectral image analysis workflow. First, the original multispectral bands and the calculated indices were stacked into a single raster image, one for each image date (n = 11) and classified using a supervised random forest model to generate binary snow/non-snow classes. From the binary classification from each image date, the pixels with 10 or 11 snow-covered days (out of 11) were grouped to create the boundaries of individual persistent snow patches.
Figure 2. UAV multispectral image analysis workflow. First, the original multispectral bands and the calculated indices were stacked into a single raster image, one for each image date (n = 11) and classified using a supervised random forest model to generate binary snow/non-snow classes. From the binary classification from each image date, the pixels with 10 or 11 snow-covered days (out of 11) were grouped to create the boundaries of individual persistent snow patches.
Remotesensing 13 01036 g002
Remotesensing 13 01036 g003 550
Figure 3. Calculating the individual tree metrics for the final analysis involved using both the field-measured and terrestrial laser scanner (TLS) point cloud-derived datasets. Individual trees detected in the TLS point cloud were compared to field-measured trees for an omission and commission analysis. Their crown height (CH), average crown diameter (CD), and crown base heights (CBH) were compared to assess the accuracy of the TLS point cloud-derived estimates. The same comparison was made between matching trees in the TLS point cloud and Structure-from-Motion (SfM)-derived point cloud datasets, with the only difference being the addition of crown volume (CV).
Figure 3. Calculating the individual tree metrics for the final analysis involved using both the field-measured and terrestrial laser scanner (TLS) point cloud-derived datasets. Individual trees detected in the TLS point cloud were compared to field-measured trees for an omission and commission analysis. Their crown height (CH), average crown diameter (CD), and crown base heights (CBH) were compared to assess the accuracy of the TLS point cloud-derived estimates. The same comparison was made between matching trees in the TLS point cloud and Structure-from-Motion (SfM)-derived point cloud datasets, with the only difference being the addition of crown volume (CV).
Remotesensing 13 01036 g003
Remotesensing 13 01036 g004 550
Figure 4. A persistent snow patch and the trees influencing it at the three different ranges: 1, 1.5, and 2 times the height (CH) of the tree. A tree was included within each respective tree shading influence area (TSIA), if its CH multiplied by 1, 1.5, or 2 was extended into the snow patch and if the bearing of its location (XY) to the patch was within the range of daily solar angles.
Figure 4. A persistent snow patch and the trees influencing it at the three different ranges: 1, 1.5, and 2 times the height (CH) of the tree. A tree was included within each respective tree shading influence area (TSIA), if its CH multiplied by 1, 1.5, or 2 was extended into the snow patch and if the bearing of its location (XY) to the patch was within the range of daily solar angles.
Remotesensing 13 01036 g004
Remotesensing 13 01036 g005 550
Figure 5. The percent of the treated and untreated regions of the study site that are classified SCA for each UAV snow-series image following a storm. This data is grouped by storm, and partitioned by forest treatment condition to illustrate the reductions in SCA within and across the storms. SCA patterns following storm events show a greater SCA reduction in the treated portion of the study site than in the untreated portion. Evident also are the relatively similar initial amounts of classified SCA on the first image date after a storm in the treated portion, contrasting the initial amounts in the untreated portion.
Figure 5. The percent of the treated and untreated regions of the study site that are classified SCA for each UAV snow-series image following a storm. This data is grouped by storm, and partitioned by forest treatment condition to illustrate the reductions in SCA within and across the storms. SCA patterns following storm events show a greater SCA reduction in the treated portion of the study site than in the untreated portion. Evident also are the relatively similar initial amounts of classified SCA on the first image date after a storm in the treated portion, contrasting the initial amounts in the untreated portion.
Remotesensing 13 01036 g005
Remotesensing 13 01036 g006 550
Figure 6. Each UAV image (n = 11) was classified into snow/non-snow pixels (A) and the resulting rasters were assembled into a single composite, in which persistent snow patches were identified (B). The orthomosaic image used as the base imagery is from storm 2, day 3 (1 May 2019) and its binary snow classification (A) was used for the classification accuracy assessment because it contained a nearly even snow/non-snow (49/51%) area distribution.
Figure 6. Each UAV image (n = 11) was classified into snow/non-snow pixels (A) and the resulting rasters were assembled into a single composite, in which persistent snow patches were identified (B). The orthomosaic image used as the base imagery is from storm 2, day 3 (1 May 2019) and its binary snow classification (A) was used for the classification accuracy assessment because it contained a nearly even snow/non-snow (49/51%) area distribution.
Remotesensing 13 01036 g006
Remotesensing 13 01036 g007 550
Figure 7. Variable importance scores calculated for the forest structure metric predictor variables used in fitting the final MARS model framework. The variable importance score ranges from 0 to 100, with higher scores indicating that the predictor was more influential in model construction (n = 500). The chart depicts the frequency with which each forest structure metric was within the respective range of the importance score. The data are partitioned by tree shading influence area (TSIA) size to better illustrate the interaction with forest structure metric. Overall, the most influential forest structure metric predictors are canopy cover (CC), mean solar radiation footprint (SRF) value, and mean trees per hectare (TPH).
Figure 7. Variable importance scores calculated for the forest structure metric predictor variables used in fitting the final MARS model framework. The variable importance score ranges from 0 to 100, with higher scores indicating that the predictor was more influential in model construction (n = 500). The chart depicts the frequency with which each forest structure metric was within the respective range of the importance score. The data are partitioned by tree shading influence area (TSIA) size to better illustrate the interaction with forest structure metric. Overall, the most influential forest structure metric predictors are canopy cover (CC), mean solar radiation footprint (SRF) value, and mean trees per hectare (TPH).
Remotesensing 13 01036 g007
Remotesensing 13 01036 g008 550
Figure 8. A collection of the simple linear relationships reflecting each model’s predictive ability for the TSIA 1.5. Each line represents the prediction accuracy of a single model, and the red points are the value pairs for the predicted vs. observed persistent snow patch sizes (m2). The mean performance statistics for the entire set of prediction models in TSIA 1.5 is R2 = 0.70 and a RMSE = 267 m2.
Figure 8. A collection of the simple linear relationships reflecting each model’s predictive ability for the TSIA 1.5. Each line represents the prediction accuracy of a single model, and the red points are the value pairs for the predicted vs. observed persistent snow patch sizes (m2). The mean performance statistics for the entire set of prediction models in TSIA 1.5 is R2 = 0.70 and a RMSE = 267 m2.
Remotesensing 13 01036 g008
Remotesensing 13 01036 g009 550
Figure 9. A comparison of trees and their shadows with different vertical structure metrics, namely CBH:CH and CV. The perceived difference in forest canopy shading between trees with high CBH:CH and low CV values (Panel A) are compared to those trees with overall larger crowns (Panel B). These canopy shading differences were expected to result in vertical forest structure metrics being significant in the modeling of persistent snow patch size.
Figure 9. A comparison of trees and their shadows with different vertical structure metrics, namely CBH:CH and CV. The perceived difference in forest canopy shading between trees with high CBH:CH and low CV values (Panel A) are compared to those trees with overall larger crowns (Panel B). These canopy shading differences were expected to result in vertical forest structure metrics being significant in the modeling of persistent snow patch size.
Remotesensing 13 01036 g009
Table
Table 1. Three storm events, subsequent image collection dates, and weather conditions. The total snowfall indicates the snow recorded for the storm event, while the mean temperatures and wind speeds are cumulative means between the image dates. Each storm event was imaged several times from full snow coverage on the first day until the snow cover nearly melted off. The individual image dates were selected considering daily weather conditions and study site access. Periods of near-freezing daily temperatures yielded longer image date intervals so that noticeable changes in snow coverage could be observed. In contrast, periods of higher mean temperature values required more frequent image date intervals as the snow melted rapidly. Weather data were recorded at the regional National Weather Service Forecast office, approximately 15 km from the study site and compiled from NOAA’s Climate Data Online service.
Table 1. Three storm events, subsequent image collection dates, and weather conditions. The total snowfall indicates the snow recorded for the storm event, while the mean temperatures and wind speeds are cumulative means between the image dates. Each storm event was imaged several times from full snow coverage on the first day until the snow cover nearly melted off. The individual image dates were selected considering daily weather conditions and study site access. Periods of near-freezing daily temperatures yielded longer image date intervals so that noticeable changes in snow coverage could be observed. In contrast, periods of higher mean temperature values required more frequent image date intervals as the snow melted rapidly. Weather data were recorded at the regional National Weather Service Forecast office, approximately 15 km from the study site and compiled from NOAA’s Climate Data Online service.
Storm DateSnow Fall (cm)Image DateMean Daily High (°C)Mean Daily Low (°C)Mean Wind Speed (m/s)
20 January 201828.922 January 20180.7−17.43.3
24 January 20187.8−15.82.1
26 January 20185.8−12.34.9
29 January 201810.0−9.83.6
27 December 201873.401 January 20190.4−14.73.3
03 January 20195.1−17.13.6
05 January 201911.1−9.21.8
14 January 20196.8−5.72.6
05 February 201993.907 February 20191.9−10.35.9
25 February 20194.1−11.03.8
04 March 201910.7−4.33.6
Table
Table 2. Accuracy of the snow/non-snow classification across the 76-ha study site. The image used to generate the error matrix included an even amount of snow and non-snow pixels, providing the least biased estimation of classification error.
Table 2. Accuracy of the snow/non-snow classification across the 76-ha study site. The image used to generate the error matrix included an even amount of snow and non-snow pixels, providing the least biased estimation of classification error.
Class ValueNon-SnowSnowTotalUser’s Accuracy
Non-Snow2361425094%
Snow3521525086%
Total271229500
Producer’s Accuracy87%94%
Overall Accuracy: 90.2%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Belmonte, A.; Sankey, T.; Biederman, J.; Bradford, J.; Goetz, S.; Kolb, T. UAV-Based Estimate of Snow Cover Dynamics: Optimizing Semi-Arid Forest Structure for Snow Persistence. Remote Sens. 2021, 13, 1036. https://doi.org/10.3390/rs13051036

AMA Style

Belmonte A, Sankey T, Biederman J, Bradford J, Goetz S, Kolb T. UAV-Based Estimate of Snow Cover Dynamics: Optimizing Semi-Arid Forest Structure for Snow Persistence. Remote Sensing. 2021; 13(5):1036. https://doi.org/10.3390/rs13051036

Chicago/Turabian Style

Belmonte, Adam, Temuulen Sankey, Joel Biederman, John Bradford, Scott Goetz, and Thomas Kolb. 2021. "UAV-Based Estimate of Snow Cover Dynamics: Optimizing Semi-Arid Forest Structure for Snow Persistence" Remote Sensing 13, no. 5: 1036. https://doi.org/10.3390/rs13051036

APA Style

Belmonte, A., Sankey, T., Biederman, J., Bradford, J., Goetz, S., & Kolb, T. (2021). UAV-Based Estimate of Snow Cover Dynamics: Optimizing Semi-Arid Forest Structure for Snow Persistence. Remote Sensing, 13(5), 1036. https://doi.org/10.3390/rs13051036

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

Article Metrics

Back to TopTop