Factors Inﬂuencing the Accuracy of Shallow Snow Depth Measured Using UAV-Based Photogrammetry

: Factors inﬂuencing the accuracy of UAV-photogrammetry-based snow depth distribution maps were investigated. First, UAV-based surveys were performed on the 0.04 km 2 snow-covered study site in South Korea for 37 times over the period of 13 days under 16 prescribed conditions composed of various photographing times, ﬂight altitudes, and photograph overlap ratios. Then, multi-temporal Digital Surface Models (DSMs) of the study area covered with shallow snow were obtained using digital photogrammetric techniques. Next, the multi-temporal snow depth distribution maps were created by subtracting the snow-free DSM from the multi-temporal DSMs of the study area. Then, snow depth in these UAV-Photogrammetry-based snow maps were compared to the in situ measurements at 21 locations. The accuracy of each of the multi-temporal snow maps were quantiﬁed in terms of bias (median of residuals, Q ∆ D ) and precision (the Normalized Median Absolute Deviation, NMAD). Lastly, various factors inﬂuencing these performance metrics were investigated. The results are as follows: (1) the Q ∆ D and NMAD of the eight surveys performed at the optimal condition (50 m ﬂight altitude and 80% overlap ratio) ranged from − 2.30 cm to 5.90 cm and from 1.78 cm to 4.89 cm, respectively. The best survey case had − 2.30 cm of Q ∆ D and 1.78 cm of NMAD; (2) Lower UAV ﬂight altitude and greater photograph overlap lower the NMAD and Q ∆ D ; (3) Greater number of Ground Control Points (GCPs) lowers the NMAD and Q ∆ D ; (4) Spatial conﬁguration and accuracy of GCP coordinates inﬂuenced the accuracy of the snow depth distribution map; (5) Greater number of tie-points leads to higher accuracy; (6) Smooth fresh snow cover did not provide many tie-points, either resulting in a signiﬁcant error or making the entire photogrammetry process impossible.


Introduction
Snowfall in most regions of Korea occurs in the form of shallow snow, a term referring to snow with depth of less than 20 cm. Shallow snowpacks pose significant impacts on ecosystems. For example, shallow snowpacks increase surface albedo, causing the depth of frozen soil to increase [1]. This subsequently leads to greater erosion of soil [2], which contains high concentration of phosphorous and nitrogen, playing important role in plant growth [3,4]. Shallow snowpacks also influence animal habitat and predation by altering food availability [5][6][7]. Therefore, to study and model the behavior of shallow snowpacks, it is important to first measure them precisely.
However, it is difficult to measure shallow snowpacks precisely because of their high spatio-temporal variability [8][9][10]. While efforts have been made to estimate the spatial variability of shallow snow depths [11][12][13], most publicly available snow depth data are provided in the format of a set of sparsely placed point measurements or of low-resolution raster data [14] that are not accurate and precise enough to be used as input for hydrological models aimed at estimating snowmelts, groundwater recharge, and soil erosion.
Unmanned Aerial Vehicles (UAVs), which recently went through dramatic technological advancement, emerged as the spotlight of the field of surveying in the last decades. Accordingly, snow depth surveying utilizing UAV was also a part of the spotlight [31]. The techniques of photogrammetry [32,33] can construct three-dimensional models of snow fields from multiple photographs taken from UAV, so it has been recently applied to estimate the snow depth [34][35][36]. Table 1 summarizes the results of studies that estimated snow depths using UAV-photogrammetry. This study also applied the UAV-photogrammetry method to estimate the snow depth, but is distinctive from the studies summarized in the Table 1 as follows: (1) The lowest NMAD validated based on 21 in situ measurements is 1.78 cm, which is significantly lower than other previous studies; (2) This study analyzed how the accuracy of the estimated snow depth varies with regard to varying conditions of photography, such as altitude, photograph overlap ratio, time of survey, and condition of GCPs.  Figure 2 shows the accumulated snow depth measured by laser snow depth gauge at the Daegwallyeong AWS station during the surveying period. There were three snow events between 26 January and 30 January 2020, and the maximum accumulated snow depth was approximately 31 cm (recorded on 30 January 2020).

Survey Instrument Specification (UAV and GNSS Receiver)
A Phantom 4 Pro 2.0 (produced by DJI) was used to take the base photographs of this study. The drone has the maximum and recommended flight time of 30 min and 15 min, respectively. The maximum allowable wind velocity is 10 m/s [45]. It has the GPS/GLONASS positioning system and the hovering accuracy is ±0.1 m (vertical) and ±0.3 m (horizontal). The camera onboard the drone (DJI FC6310) has an image dimension of 5472 × 3648, a sensor size of 13.2 mm (width) × 8.8 mm (height), and a focal length of 8.8 mm. Each of the photographs are recorded with the positioning metadata, which contains the axis roll, pitch, and yaw of the camera along with the three-dimensional coordinate of the drone.
A Trimble R10 GNSS (Global Navigation Satellite System) receiver was used to measure the three-dimensional coordinates of the GCPs and the locations of in situ snow depth measurement [46]. The instrument uses the technology of Real-Time Kinematic (RTK), Virtual Reference Station Networks (VRSN), and Real Time eXtended (RTX) to secure the measurement precision. The maximum precision is 8 mm (horizontal) and 15 mm (vertical).

Acquisition of the Snow Depth Map
Process of the snow depth map acquisition comprises: (1) placing GCPs to secure the accuracy of the survey; (2) UAV photographing and photogrammetry to obtain the DSM of accumulated snow; and (3) UAV photographing and the photogrammetry to obtain DSM of the ground surface after snow completely melted. The map of the accumulated snow was obtained by subtracting the product of process (3) from that of process (2).

GCPs
Before the snow events, 17 GCPs were chosen, where the coordinates (longitude, latitude, and height) were precisely measured using the GNSS receiver. These GCPs were chosen as the locations with easily identifiable landmarks, such as manholes on roads (Figure 3a,b). Snow on these locations was removed soon after the snow events, either by a researcher or by snow plowing vehicles operated by the local disaster management authority. Furthermore, we installed additional GCPs between the UAV surveying. A black planar object sized 50 cm × 50 cm was placed on the accumulated snow immediately prior to a set of UAV surveying (Figure 3c). Coordinates of an edge of the black planer object were measured using the GNSS receiver precisely considering the vertical location of the sensor in the instrument and the manually measured snow depth ( Figure 3c). The green circles in Figure 1c show the location of the 17 GCPs used by this study.

UAV Surveying Campaigns
The UAV surveying campaigns were repeated 37 times after and between the three snow events with varying conditions of UAV altitude, photograph overlap ratio, and surveying time of day, as summarized in Tables 2 and 3. Figure 5a-d shows the flight path and shooting locations with different overlap ratios and flight altitudes. Figure 5e shows the DSM of the study area viewed from the UAV take-off position.

Photogrammetric Processing to Obtain the DSM
Photographs taken from the UAV were fed into the photogrammetry software tool, ContextCapture, Bentley [47], to produce the DSM and orthomosaic of the study site with accumulated snow. Here, we describe the general process of photogrammetry instead of the one actually employed by the software tool, which is kept confidential. The process is composed of the following three steps:

1.
Camera calibration: this is the process of estimating camera's intrinsic parameters such as its focal length, skew, distortion, and image center. This is first done by identifying common features in multiple photographs of which real-world 3D coordinates are precisely known (such as GCPs). Then, a set of equations [48] are solved to obtain the camera intrinsic parameters. Here, some of the intrinsic parameters are precisely provided by the camera manufacturer (e.g., focal length) and can be directly adopted.
The intrinsic parameters are used to correct the distortion of the photographs and they apply to all photographs, so the accurate 2D coordinates of the pixels in all photographs can be obtained in this step.

2.
Structure from Motion (SfM): in this process, the tie-points, which are another set of features that coexists in a series of photographs, are first identified. Then, an optimization process of which the objective function is composed of a set of collinearity equations [37,49,50] are solved to obtain extrinsic parameters of the camera (e.g., position and orientation of camera corresponding to all photographs) and 3D coordinates of the tie-points. Here, the key to successful implementation is the automated algorithms for tie-point detection, which requires the comparison between all input photographs [51]. 3.
Bundle Adjustment: the GCPs, as opposed to tie-points, have the precisely known pair of 2D and 3D coordinates, which are also included in the objective function of the optimization process. Therefore, accuracy and the spatial configuration of GCPs greatly influences accuracy of photogrammetry. 4.
DSM construction. In this step, the triangular irregular network, which is composed of tie-points, is constructed and textures extracted from photographs are applied to obtain the DSM. Here, the accuracy of the DSMs may be estimated by performing the pixel-wise correlation analysis.
After the snow completely melted, a UAV survey was performed to obtain the DSM of the ground surface of the study site. A greater accuracy would have been achieved if the site was surveyed immediately before the snow event because the ground surfaces slightly deform as snow melts and runs off the ground. However, this was unfortunately impossible because the large snow events, as investigated in this study, are rare in Korea and only after this kind of large snow event occurs, the study area can be fixed. The photogrammetry was performed using the 17 GCPs, shown as yellow stars in Figure 1c. The accuracy of the DSM of the ground surface was validated at seven check points (CPs), shown as X marks in Figure 1c. The RMSE and the NMAD were 2.75 cm and 1.19 cm, respectively. Figure 6a-c show the DSM of the study site (a) with and (b) without accumulated snow, and (c) the snow depth map, respectively. Note that the local distortion of elevation tends to occur at the edge of the DSM (magenta area in Figure 6c). This is primarily because the edges are not photographed from as diverse angles as the internal area. Another reason may be the sparse placement of GCPs. The black line in Figure 6c delineates the study area where in situ measurement and analysis were performed. Most of this internal study area is not influenced by the local distortion of the DSM.

In Situ Snow Depth Measurment and Validation
Immediately after each of the UAV survey campaigns (Table 2), the snow depth was manually measured using a stainless metal ruler at and around the 21 specified locations shown as red triangles in Figure 1c. The in situ snow depths corresponding to each survey day are expressed as box plots in Figure 2. They ranged between 17.5 cm and 36 cm and were consistently greater than the measurement at the AWS stations, while the temporal trends coincide. This degree of discrepancy is not surprising considering the high spatial and temporal variability of rainfall in Korea [52,53]. The estimated snow depth values obtained from photogrammetry were validated based on these in situ snow depth measurements.
The mean of residuals (µ r , Equation (2)), tandard deviation of error (σ r , Equation (3)), and Root-Mean-Square-Error (RMSE, Equation (4)) are the typical metrics to measure the accuracy of the survey. Here, note that the term "accuracy" represents the general closeness of a measurement to the true value, which refers to both bias and precision. The term "bias" represents a measure of how far the expected value of the estimate is from the true value, and the term "precision" represents a measure of how similar the multiple estimates are to each other: where D obs,i and D U AV,i represent the observed in situ snow depth and the snow depth estimated from the UAV photogrammetry at the ith measurement location, respectively, and ∆D i (Equation (1)) represents the measurement error. These statistics (Equations (2)-(4)) can be good accuracy measurement metrics only if the measurement errors (∆D i s) are normally distributed. In addition, the uncertainty analysis results based on these statistics may be highly sensitive to a few measurement outliers, which are hard to filter out due to the ambiguous standard for detection. In such cases, the median of residuals (Q ∆D and Normalized Median Absolute Deviation (NMAD, Equation (5)) can be good alternatives to quantify the bias and the precision, which are less sensitive to outliers [39,54]: Therefore, this study used all five statistics (µ r , σ r , RMSE, Q ∆D , and N MAD) to quantify the bias and precision of snow depth estimation. Figure 7 shows a scatter plot comparing the observed snow depth (x) and the snow depth estimated from the UAV photogrammetry (y) for the dates of 4th (red circles), 6th (green triangles), and 10 (blue squares) February. The results corresponding to 29 January and 1st February were not included because the photogrammetry process could not identify enough tie points due to the extremely smooth snow surface and the corresponding high reflectivity [35]. Thus, the snow map either could not be produced (29 January) or had very low accuracy (1 February). In addition, the results corresponding to 50 m of UAV altitude and 80 percent of photograph overlap are only shown here because this condition produced the lowest NMAD. The results corresponding to all four surveying times (9 a.m., 11 a.m., 1 p.m., and 3 p.m.) are shown in Table 4. The NMAD varied between 3.11 to 4.45 cm. The causes of varying accuracy are discussed in the later sections, because the residuals shown in Figure 7 were caused by a variety of factors, and thus, difficult to explain without isolating each of those factors. Note that the NMAD represents the precision of the measurement but not bias. The RMSE, which is the standard deviation of the residuals (∆D i ), can be a good proxy for overall measurement accuracy reflecting both bias and precision. The RMSE varied between 3.19 cm and 5.97 cm. The causes of varying accuracy are discussed in the later sections, because the residuals shown in Figure 7 were caused by a variety of factors and are thus difficult to explain without isolating each of those factors.   Figure 8 shows the box plots of the (a) Q ∆D and (b) NMAD of the estimated snow depth that varies with four different conditions of UAV flight altitudes and photograph overlap ratios. Low Q ∆D and NMAD is associated with low flight altitude and high photograph overlap. This is because a greater number of tie-points can be identified at low flight altitudes and high photograph overlap.  Figure 9a,b shows the Q ∆D and the NMAD of each survey varying with the number of identified tie points. A greater number of tie points is associated with low Q ∆D and low NMAD of the estimated snow depth. The vertical variation in the scatters is associated with flight time, altitude, and overlap ratio, and the number, accuracy and spatial configuration of GCPs, which we explain in the following sections.

Influence of Spatial Density of GCPs
It is widely known that GCP spatial density affects the accuracy of photogrammetric surveys [55]. This study assumed that the same principle applies to the snow depth survey. To verify this hypothesis, this study repeated the process of intentionally removing the GCPs while performing the photogrammetric processing and recorded NMAD and Q ∆D of the estimated snow depth. Results with the lowest NMAD value (50 m height-80% overlap ratio, 4 February, 9 a.m.) were used for this experiment. Figure 10a,b show NMAD and Q ∆D of the estimated snow depth varying with the number of GCPs. The GCPs were removed evenly over the space so that the remaining GCPs are distributed equally over the space. The number of remaining GCPs varied between 4 and 17, which correspond to the spatial density of 0.4 GCP/km 2 and 1.7 GCP/km 2 , respectively. The magenta area of Figure 6c suggests sudden increases of snow depth, which strongly suggests that this area is associated with the error that is caused by the local deformation of DSM. This is a known issue of photogrammetry, which often happens at the area located outside the spatial networks of GCPs. This study performed an experiment to verify this assumption. The color of the points in Figure 10a,b represent the number of the in situ measurement locations outside the GCP network. The red points mean that most of the in situ measurement locations are located outside the GCP network, and vice versa. The vertical distribution of the color dots did not show a clear trend, which means that there are other important factors influencing the accuracy of the snow photogrammetry than the relative locations between the GCP network and in situ measurement locations. This suggests that the density and spatial distribution as well as type of GCP (i.e., road, manhole, and planar object on snow field) may be more critical factors influencing the accuracy of the snow photogrammetry. For example, most red colors were associated with the GCP network, containing many GCPs placed on snow field, which has relatively low accuracy. In addition, the GCPs were evenly placed over the study area, so no in situ measurement locations were truly further away from the GCP network, which may be another reason that the trend was not detected.

Influence of the Spatial Configuration and the Accuracy of the GCPs
The spatial configuration as well as the spatial density of GCPs is important. Figure 11c shows the kernel density of spatially distributed GCP points; the pinker the area, the more GCPs are influencing the area. Figure 11a shows an interpolated map of snow depth residual of the estimated snow depth corresponding to the measurement with the overall lowest accuracy (50 m height-80% overlap ratio, 10 February, 3 p.m.). Figure 11b shows an interpolated map of CP residual of snow-free surface DSM, which were used to obtain Figure 11d, which is the map of the snow depth residual after CP correction. Figure 11c shows the kernel density of the GCPs. It is shown that the blue area of Figure 11c (area further away from GCPs) and the purple area of Figure 11d (area with greater estimation error) coincides, suggesting that the number of GCPs used for bundle adjustment may influence the accuracy of the photogrammetry.   Significant positive bias is apparent for the surveying time of 9 a.m. and 3 p.m. One possible cause of this time-dependent accuracy is that shadows of the object generally degrade the accuracy of the photogrammetry result, and shadow sizes were greater as the surveying time was further away from noon. For example, the time lapse of the photographs existing between different swaths can be several minutes, and even in those several minutes, the size of the shadows can change significantly in early morning and late evening time [42]. In addition, clear contrast between snow and overlaying shadow induces significant positive bias in the photogrammetry result. Figure 13 shows the orthophoto and the estimated snow depth of a small fraction of the study site for 9 a.m. of 4 February. While the snow depth in the shadow seems not much different from the non-shaded area according to the orthophoto as shown in Figure 13a, the snow depth is significantly greater in the shadow area according to the photogrammetry, as shown in Figure 13c. These areas with high brightness contrast are greater for morning and evening time, which causes greater bias for the overall photogrammetry results.  Figure 14a,b show the identified tie points and the orthophoto images (gray) during the survey of 29 January and 10 February, respectively. Figure 14c,d are the magnified view of the same subset of the images. Color of the points represents the number of photos that identified corresponding pixel location as a tie point. Inspection of survey data for 29 January shows that not as many tie points could be identified as were identified from data for 10 February, because the snow field initially had a very smooth surface on 29 January. Then, wind and uneven melting caused rifts on the field, dramatically increasing the number of tie points for the survey of 10 February.

Issues with Image Saturation
The histograms (Figure 14e) of the orthomosaic subsets (Figure 14c,d) reveal that the colors of 26 January photographs are limited to within a very narrow range of the colors, while the colors of 10 February photographs are significantly more diverse.
This issue of image oversaturation as shown through the case of the 29 January may be resolved through the application of the Contrast Limited Adaptive Histogram Equalization (CLAHE) technique. Figure 15a compares the same subset of the area (shown as Figure 14a) before and after applying the CLAHE filter. The comparison shows that the filter significantly resolves the issue of image saturation. Following this finding, we performed an experiment to see if the different levels of contrast factor of CLAHE filter increases the number of tie points. Figure 15b shows the grayscale histograms of the image of the entire study area corresponding to the contrast factor (CF) of 0 (no filter), 0.001, 0.004, and 0.010, respectively. The number of the identified tie points corresponding to each filter was 30,310 (no filter), 32,789 (CF = 0.001), 33,495 (CF = 0.004), and 11,897 (CF = 0.010), respectively. The number of unusable photographs were 18, 15, 16, and 21, respectively. Even though greater contrast factor slightly increased the number of tie points, the overall photogrammetry accuracy was not enhanced. This is because a smooth snow surface occupies a primary portion of the study area, and the photographs corresponding to this area were not enhanced, even after applying the CLAHE filter, as shown by Figure 15c. However, if the study area had an uneven surface such as variable terrain or many structures, the CLAHE filter may be a good option to improve snow depth estimation accuracy. Figure 14. The location of the tie points that were identified by the program that were used in this study for the date of (a,c) 29 January 2020 (9 p.m.) and (b,d) 10 February 2020 (3 p.m.). (e) Gray histogram of photo taken on 29 January and 10 February at the same position, respectively. Black area in (a) and (b) represents the area where photographs could not be tied due to the lack of tie points.
A way to resolve this issue may be to apply enhanced algorithms for image matching and orientation specialized to resolve image saturation, such as Semi-Global Matching [56,57]. However, a further test may be necessary to see if the algorithm would work for extremely saturated images as the case of this study. Another alternative is to scatter tiny color papers (e.g., 3 cm × 3 cm) using the UAV before taking photographs, which we consider as a future work. However, this approach should be cautiously performed as to not harm the environment and experimental condition (e.g., change of albedo influencing snow melts). Footprints can also be used to create tie-points, but they bring along with them a lot of limitations if the survey field is large or in inclement weather and complex terrain conditions. Alternately, high-performance cameras with more dynamic ranges and higher resolution as well as enhanced image filters may be applied.

Conclusions
This study investigated factors influencing the accuracy of a snow depth distribution map obtained from the UAV-photogrammetry. The results suggest that low flight altitude of UAV and high overlap ratio increases the number of the tie-points, which enhances the accuracy of the snow map. As a greater number of accurate GCPs are spatially well distributed, the accuracy of the snow map increased. A smooth snow surface, evident immediately after a snow event, does not have shadows to be used as tie-points, thus, either inducing significant errors or making the entire photogrammetry process impossible.
Another important finding of this study comes from the low NMAD value corresponding to the best survey case, which is 1.78 cm (RMSE: 2.85 cm). This value is significantly lower than that attained by previous studies, which ranges between 8 cm and 42 cm. This result suggests that the abundancy, accuracy, and spatial configuration of GCPs are the controlling factors for the accurate measurement of snow depth, but not the performance of sensors and drones.