Terrain-Based Shadow Correction Method for Assessing Supraglacial Features on the Greenland Ice Sheet

The presence of shadows in remotely sensed images can reduce the accuracy of land surface classifications. Commonly used methods for removing shadows often use multi-spectral image analysis techniques that perform poorly for dark objects, complex geometric models, or shaded relief methods that do not account for shadows cast on adjacent terrain. Here we present a new method of removing topographic shadows using readily available GIS software. The method corrects for cast shadows, reduces the amount of over-correction, and can be performed on imagery of any spectral resolution. We demonstrate this method using imagery collected with an uncrewed aerial vehicle (UAV) over a supraglacial stream catchment in southwest Greenland. The structure-from-motion digital elevation model showed highly variable topography resulting in substantial shadowing and variable reflectance values for similar surface types. The distribution of bare ice, sediment, and water within the catchment was determined using a supervised classification scheme applied to the corrected and original UAV images. The correction resulted in an insignificant change in overall classification accuracy, however, visual inspection showed that the corrected classification more closely followed the expected distribution of classes indicating that shadow correction can aid in identification of glaciological features hidden within shadowed regions. Shadow correction also caused a substantial decrease in the areal coverage of dark sediment. Sediment cover was highly dependent on the degree of shadow correction (k coefficient), yet, for a correction coefficient optimized to maximize shadow brightness without over-exposing illuminated surfaces, terrain correction resulted in a 49% decrease in the area covered by sediment and a 29% increase in the area covered by water. Shadow correction therefore reduces the overestimation of the dark surface coverage due to shadowing and is a useful tool for investigating supraglacial processes and land cover change over a wide variety of complex terrain.

The presence of shadows in remotely sensed images can reduce the accuracy of land surface classifications. Commonly used methods for removing shadows often use multispectral image analysis techniques that perform poorly for dark objects, complex geometric models, or shaded relief methods that do not account for shadows cast on adjacent terrain. Here we present a new method of removing topographic shadows using readily available GIS software. The method corrects for cast shadows, reduces the amount of over-correction, and can be performed on imagery of any spectral resolution. We demonstrate this method using imagery collected with an uncrewed aerial vehicle (UAV) over a supraglacial stream catchment in southwest Greenland. The structure-from-motion digital elevation model showed highly variable topography resulting in substantial shadowing and variable reflectance values for similar surface types. The distribution of bare ice, sediment, and water within the catchment was determined using a supervised classification scheme applied to the corrected and original UAV images. The correction resulted in an insignificant change in overall classification accuracy, however, visual inspection showed that the corrected classification more closely followed the expected distribution of classes indicating that shadow correction can aid in identification of glaciological features hidden within shadowed regions. Shadow correction also caused a substantial decrease in the areal coverage of dark sediment. Sediment cover was highly dependent on the degree of shadow correction (k coefficient), yet, for a correction coefficient optimized to maximize shadow brightness without over-exposing illuminated surfaces, terrain correction resulted in a 49% decrease in the area covered by sediment and a 29% increase in the area covered by water. Shadow correction therefore reduces the overestimation of the dark surface coverage due to shadowing and is a useful tool for investigating supraglacial processes and land cover change over a wide variety of complex terrain.

INTRODUCTION
The use of high resolution imagery has revolutionized our understanding of natural landscapes. Analysis of high resolution satellite (e.g., WorldView satellites) and Uncrewed Aerial Vehicle (UAV)-based imagery allows for classification of features that may be hidden in coarser resolution remote sensing products (e.g., Landsat, MODIS) (Ryan et al., 2015;Akar et al., 2017). This is especially important in glacial terrain where dark features such as water bodies and sediment impact the energy balance of the ice surface. Many supraglacial features including stream channels, cryoconite holes, and fractures have low spatial coverage and often occur at spatial scales below high resolution satellite imagery (Ryan et al., 2015;Yang et al., 2018). Yet, these features have a disproportionately large impact on albedo. Supraglacial streams, for example, are responsible for 12% of the ice albedo variability despite only covering 2% of the ice surface (Ryan et al., 2018). In some regions of the ice sheet, supraglacial streams can also contain expanses of dark sediment deposits covering up to 25% of the channel bed furthering this albedo reduction (Leidman et al., 2020). Understanding the spatial distribution of these features can facilitate a more process-based understanding of satellite-derived albedo measurements (Moustafa et al., 2015). In turn, this improved understanding can aid climate model calculations of surface mass balance (Enderlin et al., 2014;Ryan et al., 2018) and dynamic losses from hydrofracturing and runoff induced ice acceleration (McGrath et al., 2012;Smith et al., 2021). Additionally, knowing the extent of hydrologic features on the surface of Antarctica can help improve predictions of ice shelf instability (McGrath et al., 2012). Climate models show that the areal coverage of supraglacial surface water will increase with climate warming especially in northeast Greenland, exacerbating the effect of ice sheet melting (Igneczi et al., 2016). Understanding the effect of changing land cover is therefore important for assessing the impact of ice sheets on large scale ocean circulation and sea level rise (Rennermalm et al., 2007;Bakker et al., 2016).
One issue that often complicates land surface classification is shadowing. Complex topography and high relief objects can cause large shadows which obscure terrain and reduce classification accuracy. This is especially true in the polar regions where low sun angles cause long and dark shadows on adjacent terrain. Crevassing, supraglacial stream incision, and large-scale subglacially-derived undulations can cause large topographic features to develop on the ice sheet surface in Greenland (Yi et al., 2005;Nolin and Payne, 2007;van der Veen et al., 2009). Shadows cast by these features have similar spectral properties as dark sediment within the visible spectrum which potentially leads to over-classification (Ryan et al., 2018). Even so, the impact of shadowing on melting of the ablation zone is rarely addressed, let alone quantified, in remote sensing studies of the Greenland Ice Sheet.
Several methods have been developed to detect and remove shadows from remotely sensed imagery (Dare, 2005;Adeline et al., 2013;Shahtahmassebi et al., 2013). Identifying shadowed regions by hand can be costly, time consuming, and inaccurate (Tarko et al., 2018). A threshold analysis using an index of spectral bands or an invariant color model is commonly used, but can lead to over-classification of shadows as dark colored objects, requires multispectral imagery, and ignores the causes of the observed shadows (Sarabandi et al., 2004). Traditional shaded relief methods cannot remove cast shadows on adjacent land surfaces and only account for decreased reflectance due to a particular slope not being in direct view of the sun, also known as self-shadowing (Shahtahmassebi et al., 2013). Additionally, shadow relief methods such as cosine functions over-correct in shadowed regions since they assume Lambertian scattering (Soenen et al., 2005). A commonly used shadow relief method is the C-correction. This method is similar to a cosine-correction but adds a coefficient C to both the cosine zenith and cosine incident angle (Teillet et al., 1982). C is defined as the y-intercept divided by the slope of the cosine incident vs digital number linear regression. The C-correction is computationally inexpensive and has shown to perform better than other Digital Elevation Model (DEM) based corrections in some situations (McDonald et al., 2000). However, its effectiveness is dependent on ground cover type and can cause differential correction for each spectral band. Shortcomings with the C-correction method have been remedied with complex methods utilizing a support vector machine (SVM) supervised classifiers trained with line-of-sight analysis of airborne laser scanners (Tolt et al., 2011). These methods have high coherence with visual interpretation however they require high resolution DEMs (0.25 m) and hyperspectral imagery as well as extensive computational processing time. A simpler method for removing cast shadows that uses widely available software is therefore necessary for accurately analyzing remotely sensed imagery.
Here, we develop a relatively simple method to remove cast shadows from remotely sensed imagery. We demonstrate this method using a high resolution DEM derived from structurefrom-motion analysis of UAV imagery. Structure-from-motion analysis of remotely sensed stereopair images has proliferated the availability of high resolution DEMs, especially in the polar regions (Noh and Howat, 2015), however, to our knowledge, these DEMs have not been used to remove shadows and improve classification of supraglacial features or predict shadow locations. We apply this method to imagery collected in 2016 over the Greenland Ice Sheet. The shadow removal algorithm utilizes ArcGIS Pro's solar radiation toolset to identify shadowed regions and performs a supervised classification on the original and corrected imagery to quantify the impact of our shadow removal process on the classification of water, sediment, and bare ice surfaces. Finally, we show how the use of DEM-derived shadow models can improve classification of remotely sensed imagery.

DATA
In order to perform the shadowing correction described below, two datasets are needed; 1) a high resolution nadir image with a known image collection time, and 2) a co-located digital elevation model of the region with a spatial resolution fine enough to adequately represent the topography that is causing the shadowing. For this analysis, imagery was collected at 12:30 West Greenland Summer Time (UTC-2 h) on July 21st, 2016 near the ice edge in southwest Frontiers in Remote Sensing | www.frontiersin.org August 2021 | Volume 2 | Article 690474 Greenland (67.15N, 50.00W, the so-called Point 660, 30 km east of the town of Kangerlussuaq). A 3DR Solo quadcopter UAV mounted with a Canon S110 digital camera (12 MB) collected a total of 1,196 RGB images in RAW format over four flights within 1 hour of solar noon. Four ground control points marked with orange plastic rectangles (0.3 m × 0.4 m) were simultaneously georeferenced with a Trimble R7 survey-grade GPS with a 20 min static occupation to allow for higher positional accuracy. The UAV was flown at a maximum of 100 m above the ice surface with a minimum of 75% image overlap. As a result, the images had a total footprint of 1.3 km 2 .

METHODS
UAV imagery was processed in several steps to produce a corrected and classified image (Box A-E in Figure 1). First, UAV imagery (Box A, Figure 1) was processed with structure-from-motion (using Agisoft Metashape Pro version 1.6.5) to produce a 0.034 m resolution orthomosaic RGB image and a 0.068 m resolution DEM (Box B, Figure 1). These were the highest possible output resolutions for the imagery. Camera calibration parameters and bundle adjustments were automatically set to reduce ground control point (GCP) error accounting for the 0.05 m accuracy of the GNSS and digitization (Table 1). This resulted in a mean GCP georeferencing error of 0.15 m. Ice velocity, derived from the GCP static occupations, showed a negligible effect during the image collection period compared to the georeferencing accuracy. This method of structure-from motion bundle adjustment and DEM generation has proved to result in highly accurate elevation data . The DEM, resampled to 0.68 m resolution with cubic convolution to conserve computing time, was used as the input to the Area Solar Radiation tool in ArcGIS Pro 2.1 to create a radiation raster image (Box C, Figure 1). For the Area Solar radiation tool, we used a transmissivity value of 0.545, which is representative for the atmospheric conditions in this area (van den Broeke et al., 2008). The Area Solar Radiation tool was then set to calculate solar radiation every 6 minutes from 12:12 to 12:48 local time (UTC-2h) July 21st, 2016. The rescaled area solar radiation output image is hereafter referred to as the radiation raster. After the radiation raster was created, a novel shadow correction method (Box D, Figure 1) scaled the values to the difference in brightness between shadowed and illuminated areas. Pixels with below average (below the arithmetic mean) solar radiation (R mean ) were considered to be in shadow used to create a shadow mask ( Figure 1, Steps 1, 2 in Box D). All three of the orthomosaic bands were then averaged together to create a composite black and white image hereafter termed an illumination surface (I, Figure 1, Step 3). This averaging process ensures that the same correction factor is applied to all bands and therefore reduces the potential for blue shifting in the corrected image. Then the average value of the illumination surface (I) within shadowed regions (R < R mean ) and unshadowed regions (R ≥ R mean ) were calculated using a zonal statistic (I shad and I illu respectively, Figure 1, Step 4). The solar radiation raster (R) was normalized to values between −1 and 1 using the equation: where R norm is the normalized radiation raster, R is the radiation value, and R max and R min were the maximum and minimum values of the radiation output, respectively ( Figure 1, Step 5).
Using the values for I illu , I shad , and the mean value of the illumination surface (I mean ), a correction surface was made following equation: where k is an arbitrary shadowing coefficient greater than 0 ( Figure 1, Step 6). Pixels with R > R norm are given a correction FIGURE 1 | Flow chart of processing steps. R min , R mean , and R max are the minimum, mean, and maximum values of the radiation surface (R), respectively. Imean is the mean value of the illumination surface (composite black and white image). I shad and I illu are the mean values of the illumination surface within the shadowed and illuminated areas respectively based on the area solar radiation shadow mask. R norm is the normalized radiation value used as a threshold for determining the amount of correction applied to each pixel. k is an arbitrary scaling factor.
Frontiers in Remote Sensing | www.frontiersin.org August 2021 | Volume 2 | Article 690474 factor scaled to the difference in illumination between I illu (the illumination for unshadowed regions) and I mean (the average value of the illumination surface). Pixels with R ≤ R norm are given a correction factor scaled to the difference between I mean and I shad (the average illumination for shadowed regions). This thresholding of the correction factor based on R norm was applied since shadowing causes a Poisson distribution of radiation values with shadowed pixels decreasing in radiation more than optimally oriented pixels increase in radiation. However, since the amount of correction decreases as R norm of each pixel approaches 0, the precise threshold value has very little impact on the overall correction. The correction surface was then subtracted from each band to create the corrected image (Corrected Original − Correction Surface) (Figure 1, Step 7). This process causes areas with below average radiation to be brightened and areas with above average radiation to be dimmed.
In addition to the novel shadow removal method described above, the original imagery was processed using a traditional C-correction shadow removal algorithm (McDonald et al., 2000) for comparison. The C-correction method uses the equation: where L n is the corrected radiance value, L is the original radiance value, θ is the solar zenith angle, i is the incident angle, and C is equal to a/b where a is the y-intercept and b is the slope of the cosine incident vs. digital number linear regression (Teillet et al., 1982;Soenen et al., 2005). Incident angles were determined using the ArcGIS slope and aspect tool and the solar elevation angle (90°-solar zenith angle), solar azimuth, and solar declination from the National Oceanic and Atmospheric Association (NOAA) Solar Calculator (44.21°, 161.4°, and 20.3°, respectively) (https://www.esrl.noaa.gov/gmd/grad/solcalc/). Digital numbers were averaged for all three bands and plotted against the cosine of the incident angle (R 2 0.31) resulting in a C value of 1. 116 which was used for the C-correction shadow removal method. The original imagery, the radiation based corrected imagery, and the C-corrected imagery was classified using a maximum likelihood supervised classification. This method, easily implemented using ArcGIS Pro, has shown to be a highly accurate way of classifying land cover (Sisodia et al., 2014). A total of 91 training polygons were used; 36 bare ice, 34 water, and 21 sediment polygons for a total of 1,699,961, 290,459, and 292,312 cells each, respectively. Special consideration was taken to make sure that bare ice polygons covered a wide variety of ice surfaces including steep surfaces and weathered ice. Each of the three classes are fairly spectrally unique and have similar non-Lambertian properties therefore variable viewing angles likely have little impact on the classification. An accuracy assessment was then done on 450 randomly generated points (150 from each class) to determine the total accuracy, producer's accuracy, and user's accuracy for the three images. Producer's accuracy is a measure of how often a feature on the ground is accurately identified whereas user's accuracy is how often a classification on the map accurately depicts the actual feature. In other words, producer's error measures errors of omission whereas user's error measures errors of commission. For the sake of most glaciological research, user's accuracy is more relevant. The overall accuracy denotes the number of correctly identified points divided by all 450 points. Lastly, a rough estimate of surface albedo was determined based on the classified images. The albedo value was calculated using the albedo of clean ice (0.55), shallow water (0.26), and concentrated cryoconite deposits (0.09) determined by Ryan et al. (2018) weighted by the spatial extent of each class. The albedo values from Ryan et al. (2018) are based on similar UAV imagery in southwest Greenland under similar atmospheric conditions therefore no additional atmospheric correction is necessary.

RESULTS
Supraglacial streams and sediment were pervasive throughout the study area. Sediment seemed to consolidate into supraglacial floodplains and melt ponds similar to results found by Leidman et al. (2020). Smaller streams tended to have relatively more sediment versus larger and more highly incised channels. The slope for the study area was 16.9°± 10.7°(mean ± standard deviation). This topography resulted in extensive dark shadows throughout the study area, especially along the banks of supraglacial streams ( Figure 2C). The output solar radiation raster had a left skew with a mode of 351 WH/m −2 and standard deviation of 60.9 WH/m −2 . Within the flight footprint, 40.3% of the area had below average solar radiation.
The C-correction method resulted in some shadowed regions being overexposed with reflectance values shifted disproportionately towards more blue hues, i.e. blue-shifting ( Figure 2D). Cast shadows remained dark after the correction which was particularly noticeable within stream channels. Additionally, the contrast of non-shadowed regions decreased and received a browner coloration relative to the original image. This was in contrast to the area solar radiation based shadow removal method which saw no such coloration change. A value of 2.1 was used for the k coefficient in order to maximize shadow correction while minimizing over-exposure (Supplementary 1 | Accuracy of each classification type. The user's accuracy indicates how often a classification on the map accurately depicts the actual feature, also known as errors of commission. Kappa, or the Cohen's kappa statistic, is a measure of the level of agreement between the classification and the assessment points that factors in the possibility of the relationship occurring by chance. Figures S1, S2). The resulting image displayed almost no overexposed areas or blue-shifting. It also removed cast shadows and maintained the coloration of illuminated areas ( Figure 2E). As a result, features within the scene were much more easily identifiable than in the original and C-corrected images. Classification of the ice surface using the original, C-corrected, and area solar radiation corrected imagery showed that bare-ice covers the vast majority of the scene (92, 89, and 93% respectively. Only 3% of cells changed class after the shadow correction was performed. The majority of the changed classification was reclassification of sediment to bare-ice (neon green, Figure 3). Most of the reclassified areas (86%) occurred within shadowed regions indicating that the classification scheme more easily identified surface classes after shadows were brightened. This can also be seen with the change in classification confidence with increased shadow correction. As the k value is increased, the likelihood of a correct classification significantly increases ( Figure 4A). The degree to which classification confidence correlates with k, however, is likely dependent on the accuracy of the training polygons in identifying the spectral properties of each class. Additionally, as k increases, the proportion of cells classified as sediment steadily decreases ( Figure 4B). At a k value of 2.1, the percentage of the image extent covered by sediment decreases by 49% compared to the original imagery. The proportion of cells classified as water also increases although to a lesser extent (+29%). The shadow correction rarely results in pixels changing classification to sediment (0.07% of the time).
The accuracy assessment for the three classified rasters (i.e., original, C-corrected, area solar radiation corrected with k 2.1) shows that, overall, both shadow correction methods reduce the total accuracy of the scene (Supplementary Table S1). However, the accuracy decrease is significantly less for the area solar radiation based method compared to the C-correction method and the accuracy of the area solar radiation method is within the margin of error for the accuracy assessment. Shadow removal resulted in a large increase in the user's accuracy for water pixels but a decrease in user's accuracy of bare ice and sediment classification. The high Cohen's kappa statistic for the original image and the area solar radiation classification indicates the high total accuracy for both classifications is likely not merely the result of chance. This is less so for the C-correction in which the level of agreement for the kappa statistic is considered moderate.
The average albedo of the original, C-corrected, and area solar radiation corrected imagery was determined using a weighted average of each surface class multiplied by the albedo value reported by Ryan et al. (2018) (0.55 for bare ice, 0.26 for shallow water, and 0.09 for concentrated cryoconite deposits, which are likely analogous to the sediments seen in this imagery). Using these values, the original and C-correction imagery both have albedo values of 0.516 ( Figure 4A). After the solar radiation based shadow removal with a k coefficient of 2.1, the scene has an albedo of 0.523. This indicates that the increase in classified bare ice within the scene after shadow removal increased the classification-derived albedo by 1.4%.

DISCUSSION
The application of the area solar radiation tool for removing terrain effects in remotely sensed imagery proved to be an effective means of classifying land cover without the downsides of overexposure, color shifting, or cast-shadow preclusion. This novel approach is able to lessen the exposure of highly illuminated areas and increase the exposure of areas within cast shadows by taking a mechanistic approach to shadow correction and using available topographical data to calculate sun exposure throughout a scene. Additionally, the shadow correction was able to be performed with widely used GIS software, in this case ESRI ArcGIS Pro, with minimal processing. While not tested here, processing steps involving ArcGIS' Area Solar Radiation tool could potentially be implemented using QGIS' SAGA "Potential Incoming Solar Radiation" or SEBE (Solar Energy on Building Envelopes) algorithms in order to make the workflow completely open source (Lindberg et al., 2015). As such, classification can be performed quickly with insignificant loss of accuracy and increased confidence.
The shadow correction algorithm presented here resulted in a decrease in overall accuracy of the classification by 1.7% as well as a and drop in user's accuracy for sediment and bare ice. This is contrary to visual evaluation that suggested that the corrected image represented the distribution of classes closer to expectation. While the user's accuracy decrease is likely more relevant to glaciological studies trying to identify features on the ice surface, a similar small decrease in accuracy is seen in the producer's accuracy (85, 95, and 90% for bare ice, water, and sediment respectively). The accuracy drop for the area solar radiation correction method was less than that of the C-correction method, however, the fact that accuracy is decreasing begs the question of why this is occurring while classification confidence is increasing. Firstly, the classification of the accuracy assessment points into the three classes was based on the manual visual interpretation and is therefore subjected to human errors. Secondly, the correction slightly decreases the exposure of highly illuminated areas and illuminated areas composed nearly 60% of the scene. Thus, it is possible that small accuracy decreases within the illuminated areas could have caused the overall accuracy to decrease despite accuracy increasing within shadowed regions. Additionally, the 450 points used for the accuracy assessment may not have been sufficient for determining the accuracy of a scene composed of 1.13 billion pixels. Potentially, a more intensive classification method such as decision trees or support vector machines could improve this classification accuracy (Otukei and Blaschke, 2010). Qualitatively, the scene after the correction is applied is more representative of how supraglacial features are distributed on the ice sheet based on the authors' experience in the field however, the accuracy assessment suggests that this correction is probably best suited for visual interpretation of a scene and individual feature identification rather than quantitative analysis of bulk land cover change over the entire ablation zone.
The change in the distribution of bare ice, water, and sediment before and after the area solar radiation shadow correction shows that shadows complicate classification on dark surfaces and potentially lead to an over-classification of sediment and other dark surfaces. This overestimation of distributed and cryoconite sediments in shadowed regions was also noted (albeit not quantified) by Ryan et al. (2018). This potential overestimation of sediment cover may result in erroneously low calculations of surface albedo when calculated based on area weighted averages. This underestimation of albedo may also cause assessments of Greenland Ice Sheet melt rates to be slightly overestimated, especially in areas with steep topography or in studies using imagery taken far from solar noon. Analysis of ICESat data has shown that the surface roughness and average slope of the ice surface decrease with elevation in Greenland (Yi et al., 2005) and therefore the errors caused by shadowing likely decreases with elevation as well. Olson and Rupper (2019) showed that shadows on valley glaciers, especially at high latitudes, can decrease incoming solar radiation on the ice surface by over 5 W/m2; however, our work shows that accounting for shadows in remotely sensed imagery is also important for understanding the energy balance of ice sheets and ice caps. Changes in the albedo of the Greenland Ice Sheet by 1% have been related to a change in surface mass balance by 27 Gt/yr (Dumont et al., 2014), therefore errors in albedo on the scale estimated in this study (1.7%) would potentially have substantial impacts on predictions of sea level rise.
Understanding and correcting for shadows is also important for understanding the hydrology of the Greenland Ice Sheet. River widths derived from satellite imagery is a common method for mapping supraglacial stream networks and moulin meltwater inputs (Pavelsky and Smith, 2008;Yang et al., 2016). However, some of the most intense shadows seen in the imagery presented here occurred in highly incised regions near moulins. Therefore, special care should be taken to make sure that cast shadows do not result in incorrect width estimates. Hochreuther et al. (2021) found that cast shadows resulted in an overclassification of supraglacial lake pixels in Sentinel-2 imagery and applied a shadow mask from the R-package "insol" (Corripio, 2003) to mask out these regions. While (Hochreuther et al., 2021) did not use "insol" to improve the classification as we have done here, it shows that the package could be used similarly to the method we described here to perform shadow corrections using open source software. Additionally, (Liang et al., 2012), found that shadows from clouds can lead to a significant number of false positives when using automated classification schema of supraglacial lakes. While cloud shadows are not addressed by the area solar radiation based shadow removal method, this shows that shadows can potentially lead to erroneous assessments of the impact of supraglacial lakes on ice sheet surface mass balance.
The use of the area solar radiation based shadow removal method has some potential limitations. For one, it requires a high resolution digital elevation model. The increasing availability of stereopair imagery and structure-from-motion software though means that this method is becoming increasingly viable, even in remote, and extreme environments such as polar regions. The area solar radiation tool in ArcGIS is also computationally intensive, especially for areas with large spatial extents and may be prohibitively so if applied to a very large area with high spatial resolution. The area solar radiation shadow removal method is able to correct imagery of a wide range of bands; however, it is unclear if this method may cause problems when applied to imagery far from the visible range. The intensity of the shadow correction needed for a given image will depend on the exposure of the camera, local weather conditions, and the extent of the shadows themselves therefore the k coefficient would need to be scaled appropriately via trial and error and may differ from the value of 2.1 that was ideal for the scene used in this study. Additionally, edge effects may come into play for imagery collected in which topography outside the scene causes shadowing within the image extent. Despite these limitations, our method likely can be applied to a wide range of imagery with varying pixel resolution, wavelength ranges, and surface features. While the use of this algorithm was tested over the Greenland Ice Sheet, this method can be applied to any surface and may bring insights when applied to other areas with complex topography such as forests and urban areas. This method however does not preserve reflectance values and therefore cannot be applied for analysis that requires band indexing or image textures to be preserved.
This research shows that shadows can result in the overclassification of dark surfaces and, as a result, lead to an overestimation of solar absorption for a landscape. For the Greenland Ice Sheet, this error translates to shadows being classified as sediment instead of bare ice on the ice sheet surface leading to lower estimates of albedo when albedo is derived from land cover classification. If not accounted for, this error might skew our understanding of sediment transport processes and the impact of climate change on supraglacial hydrology and microbiology. This is especially prescient in Antarctica where the presence of water bodies on the ice surface is closely tied to the stability of ice shelves (McGrath et al., 2012). Increased ponding of water on surface of the Larsen A and B ice shelves lead to increased hydrofracturing and their eventual disintegration. Therefore, accurately classifying the extent of melt ponds on the surface of Larsen C, Amery, and Getz Ice Shelves would likely improve predictions of when these ice shelves may collapse and further contribute to increased dynamic ice losses.

CONCLUSION
Shadows complicate the classification of land surfaces in remotely sensed imagery, especially in the Greenland ablation zone where extensive crevassing and supraglacial stream channels create significant topographic heterogeneity. This is compounded by low sun angles and limited availability of ground-truthed surface classification measurements. We show here a novel method of correcting for this shadowing that models shadow locations based on a digital elevation model. This method can overcome some of the pitfalls associated with other correction methods such as overexposure, color shifting, and underexposed cast shadows without significantly reducing classification accuracy. Imagery not corrected for shadows can lead to an over-estimation of the extent of dark surface features such as supraglacial debris cover and therefore may lead to underestimating albedo and thereby incorrect assessments of the surface energy balance. Based on our results, we recommend that image classification methods applied to the Greenland Ice Sheet apply the shadow correction method presented here with a k value of 2.1, to avoid underestimating dark sediment material by up to 49%. Incorporating this shadow model based correction to high resolution imagery can improve our understanding of the overall processes occurring on the ice sheet surface and help us better predict how these processes will change with increased melting.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.
The DEM and Orthomosaic Imagery from the UAV flight have been submitted to the Arctic Data Center under creative commons public domain usage rights.

AUTHOR CONTRIBUTIONS
SL: Collected field data, ran the structure-from-motion software, developed the algorithm, wrote the majority of Frontiers in Remote Sensing | www.frontiersin.org August 2021 | Volume 2 | Article 690474 the manuscript, and edited the manuscript based on coauthor and editor responses. ÅR: Provided input for scientific ideas and contributed to editing and writing draft manuscripts. RL: Provided input for scientific ideas and contributed to editing. MC: Collected field data and contributed to editing the manuscript.

FUNDING
Funding for this research was provided by the National Science Foundation Graduate Research Fellowship Program and the NASA Cryospheric Science Program Grants #NNX14AH93G and #80NSSC19K0942.