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.

1 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 structure-from-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.

2 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 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².

3 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 (Smith et al., 2016). 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.

FIGURE 1

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.

TABLE 1

TABLE 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.

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:

$R_{norm}=\frac{2\left(R\right)-R_{\min }}{R_{\max }-R_{\min }}-1(1)$

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:

$CorrectionValue=\left\{\underset{k\left(I_{mean}-I_{shad}\right)R_{norm}if{R}_{norm}\le 0}{k\left(I_{illu}-I_{mean}\right)R_{norm}if{R}_{norm}>0}\right\}(2)$

where k is an arbitrary shadowing coefficient greater than 0 (Figure 1, Step 6). Pixels with R > R_norm are given a correction 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:

$L_n=L\frac{\cos \theta +C}{\cos i+C}(3)$

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² = 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.

4 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⁻² and standard deviation of 60.9 WH/m⁻². Within the flight footprint, 40.3% of the area had below average solar radiation.

FIGURE 2

FIGURE 2. (A) Map of Greenland showing the location of UAV imagery collection shown as yellow star. (B) Orthomosaic UAV imagery before any shadow correction. The location of inset maps C-E is shown with a red rectangle and the location of the ground control points are shown with yellow dots. (C) Zoomed in imagery without any shadow correction. (D) Imagery corrected with a traditional C-correction method. (E) Imagery corrected with the area solar radiation based method.

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 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).

FIGURE 3

FIGURE 3. (A) Supervised classification map with coloration indicating the classification of the original vs. area solar radiation corrected imagery. Black box indicates the location of the inset map (panel C). (B) Proportion of each class. Bare ice to bare ice is not shown but makes up the remainder of the imagery extent (91%). (C) Zoomed in view of the classification map showing one of the larger supraglacial streams in the scene.

FIGURE 4

FIGURE 4. (A) The average likelihood of a correct classification for each pixel and the scene albedo as a function of the k value applied to the area solar radiation correction. The likelihood of a correct classification is a sigmoidal function of the ArcGIS Maximum Likelihood Classification output uncertainty value. The albedo was calculated using a weighted average of the albedo for clean ice (0.55), shallow water (0.26), and concentrated cryoconite deposits (0.09) from Ryan et al. (2018) weighted by the proportion of pixels in each class. The blue dot indicates that albedo derived from the C-Correction classification. (B) The image extent classified as sediment (brown) or water (blue) as a function of the shadow coefficient (k) applied to the area solar radiation correction method.

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%.

5 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 over-classification 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.

6 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 the manuscript, and edited the manuscript based on co-author 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.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are thankful for the support of Polar Field Services, UNAVCO, Sarah W. Cooley, Rohi Muthyala, Laurence C. Smith, and Dino Guthrie. We acknowledge the Greenlandic people and their Elders both past and present, as well as future generations who continue to be stewards of this land.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsen.2021.690474/full#supplementary-material

Supplementary Figure S1UAV imagery with differing values for the k coefficient. The value of k = 2.1 was selected for analysis in this study due to the areas of overexposure on north facing slopes seen in images with k > 2.1.

Supplementary Figure S2Same as the k = 3.4 image of Supplementary Figure S1 with yellow ovals indicating areas of overexposure and blue-shifting that informed the decision to use a k value of 2.1 for the analysis.

Supplementary Table S1Parameters used to generate the DEM in Agisoft Metashape Pro. Camera parameters were automatically calibrated to minimize ground control point error.

References

Adeline, K. R., Chen, M., Briottet, X., Pang, S. K., and Paparoditis, N. (2013). Shadow Detection in Very High Spatial Resolution Aerial Images: A Comparative Study. ISPRS J. Photogrammetry Remote Sensing 80, 21–38. doi:10.1016/j.isprsjprs.2013.02.003

CrossRef Full Text | Google Scholar

Akar, A. L. P. E. R., Gökalp, E. R. T. A. N., Akar, Ö. Z. L. E. M., and Yılmaz, V. (2017). Improving Classification Accuracy of Spectrally Similar Land Covers in the Rangeland and Plateau Areas with a Combination of WorldView-2 and UAV Images. Geocarto Int. 32, 990–1003.

CrossRef Full Text | Google Scholar

Bakker, P., Schmittner, A., Lenaerts, J. T. M., Abe-Ouchi, A., Bi, D., Broeke, M. R., et al. (2016). Fate of the Atlantic Meridional Overturning Circulation – Strong Decline under Continued Warming and Greenland Melting. Geophys. Res. Lett. doi:10.1002/2016GL070457

CrossRef Full Text | Google Scholar

Corripio, J. G. (2003). Vectorial Algebra Algorithms for Calculating Terrain Parameters from DEMs and Solar Radiation Modelling in Mountainous Terrain. Int. J. Geographical Inf. Sci. 17, 1–23.

CrossRef Full Text | Google Scholar

Dare, P. M. (2005). Shadow Analysis in High-Resolution Satellite Imagery of Urban Areas. Photogrammetric Eng. Remote Sensing 71, 169–177.

CrossRef Full Text | Google Scholar

Dumont, M., Brun, E., Picard, G., Michou, M., Libois, Q., Petit, J.-r., et al. (2014). Contribution of Light-Absorbing Impurities in Snow to Greenland’s Darkening since 2009. Nat. Geosci. 7, 509–512. doi:10.1038/NGEO2180

CrossRef Full Text | Google Scholar

Enderlin, E. M., Howat, I. M., Jeong, S., Noh, M.-J., van Angelen, J. H., and van den Broeke, M. R. (2014). An Improved Mass Budget for the Greenland Ice Sheet. Geophys. Res. Lett. 14, 866–872. doi:10.1002/2013GL059010.Received

CrossRef Full Text | Google Scholar

Hochreuther, P., Neckel, N., Reimann, N., and Humbert, A. (2021). Fully Automated Detection of Supraglacial Lake Area for Northeast Greenland Using Sentinel-2 Time-Series. Remote Sensing 13, 1–24.

CrossRef Full Text | Google Scholar

Igneczi, A., Sole, A. J., Livingstone, S. J., Leeson, A. A., Fettweis, X., Selmes, N., et al. (2016). Northeast Sector of the Greenland Ice Sheet to Undergo the Greatest Inland Expansion of Supraglacial Lakes during the 21st century. Geophys. Res. Lett. 43, 9729–9738. doi:10.1002/2013GL058740.Received

CrossRef Full Text | Google Scholar

Leidman, S. Z., Rennermalm, Å. K., Muthyala, R., Qizhong, G., and Overeem, I. (2020). The Presence and Widespread Distribution of Dark Sediment in Greenland Ice Sheet Supraglacial Streams Implies Substantial Impact of Microbial Communities on Sediment Deposition and Albedo Geophysical Research Letters. Geophys. Res. Lett. 48, 1–12. doi:10.1029/2020GL088444

CrossRef Full Text | Google Scholar

Liang, Y.-l., Colgan, W., Lv, Q., Steffen, K., Abdalati, W., Stroeve, J., et al. (2012). A Decadal Investigation of Supraglacial Lakes in West Greenland Using a Fully Automatic Detection and Tracking Algorithm. Remote Sensing Environ. 123, 127–138. doi:10.1016/j.rse.2012.03.020

CrossRef Full Text | Google Scholar

Lindberg, F., Jonsson, P., Honjo, T., and Wa, D. (2015). Solar Energy on Building Envelopes–3D Modelling in a 2D Environment. Solar Energy 115, 369–378. doi:10.1016/j.solener.2015.03.001

McDonald, E. R., Wu, X., Caccetta, P., and Campbell, N. (2000). Illumination Correction of Landsat TM Data in South East NSW. Proc. Tenth Australas. Remote Sensing Conf. 22, 19–23.

Google Scholar

McGrath, D., Steffen, K., Rajaram, H., Scambos, T., Abdalati, W., and Rignot, E. (2012). Basal Crevasses on the Larsen C Ice Shelf, Antarctica: Implications for Meltwater Ponding and Hydrofracture. Geophys. Res. Lett. 39, 1–6. doi:10.1029/2012GL052413

CrossRef Full Text | Google Scholar

Moustafa, S. E., Rennermalm, A. K., Smith, L. C., Miller, M. A., Mioduszewski, J. R., Koenig, L. S., et al. (2015). Multi-modal Albedo Distributions in the Ablation Area of the Southwestern Greenland Ice Sheet. Cryosphere 9, 905–923. doi:10.5194/tc-9-905-2015

CrossRef Full Text | Google Scholar

Noh, M. J., and Howat, I. M. (2015). Automated Stereo-Photogrammetric DEM Generation at High Latitudes: Surface Extraction with TIN-Based Search-Space Minimization (SETSM) Validation and Demonstration over Glaciated Regions. GIScience and Remote Sensing 52, 198–217. doi:10.1080/15481603.2015.1008621

CrossRef Full Text | Google Scholar

Nolin, A. W., and Payne, M. C. (2007). Classification of Glacier Zones in Western Greenland Using Albedo and Surface Roughness from the Multi-Angle Imaging SpectroRadiometer (MISR). Remote Sensing Environ. 107, 264–275. doi:10.1016/j.rse.2006.11.004

CrossRef Full Text | Google Scholar

Olson, M., and Rupper, S. (2019). Impacts of Topographic Shading on Direct Solar Radiation for valley Glaciers in Complex Topography. The Cryosphere 13, 29–40.

CrossRef Full Text | Google Scholar

Otukei, J. R., and Blaschke, T. (2010). Land Cover Change Assessment Using Decision Trees, Support Vector Machines and Maximum Likelihood Classification Algorithms. Int. J. Appl. Earth Observation Geoinformation 12, 27–31. doi:10.1016/j.jag.2009.11.002

CrossRef Full Text | Google Scholar

Pavelsky, T. M., and Smith, L. C. (2008). RivWidth: A Software Tool for the Calculation of River Widths from Remotely Sensed Imagery. IEEE Geosci. Remote Sensing Lett. 5. doi:10.1109/LGRS.2007.908305

CrossRef Full Text | Google Scholar

Rennermalm, A. K., Wood, E. F., Weaver, A. J., Eby, M., and Déry, S. J. (2007). Relative Sensitivity of the Atlantic Meridional Overturning Circulation to River Discharge into Hudson Bay and the Arctic Ocean. J. Geophys. Res. Biogeosciences 112, 1–12. doi:10.1029/2006JG000330

CrossRef Full Text | Google Scholar

Ryan, J. C., Hubbard, A. L., Box, J. E., Todd, J., Christoffersen, P., Carr, J. R., et al. (2015). UAV Photogrammetry and Structure from Motion to Assess Calving Dynamics at Store Glacier , a Large Outlet Draining the Greenland Ice Sheet. The Cryosphere 9, 1–11. doi:10.5194/tc-9-1-2015

CrossRef Full Text | Google Scholar

Ryan, J. C., Hubbard, A., Stibal, M., Irvine-Fynn, T. D., Cook, J., Smith, L. C., et al. (2018). Dark Zone of the Greenland Ice Sheet Controlled by Distributed Biologically-Active Impurities. Nat. Commun. 9, 1–10. doi:10.1038/s41467-018-03353-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarabandi, P., Yamazaki, F., Matsuoka, M., and Kiremidjian, A. (2004). Shadow Detection and Radiometric Restoration in Satellite High Resolution Images. IEEE Int. IEEE Int. IEEE Int. Geosci. Remote Sensing SymposiumIGARSS ’04. Proc. 6, 3744–3747. doi:10.1109/IGARSS.2004.1369936

CrossRef Full Text | Google Scholar

Shahtahmassebi, A., Yang, N., Wang, K., Moore, N., and Shen, Z. (2013). Review of Shadow Detection and De-shadowing Methods in Remote Sensing. Chin. Geographical Sci. 23, 403–420. doi:10.1007/s11769-013-0613-x

CrossRef Full Text | Google Scholar

Sisodia, P. S., Tiwari, V., and Kumar, A. (2014). Analysis of Supervised Maximum Likelihood Classification for Remote Sensing Image. Int. Conf. Recent Adv. Innov. Eng. 1–4.

CrossRef Full Text | Google Scholar

Smith, L. C., Andrews, L. C., Pitcher, L. H., Overstreet, B. T., Rennermalm, Å. K., Cooper, M. G., et al. (2021). Supraglacial River Forcing of Subglacial Water Storage and Diurnal Ice Sheet Motion. Geophys. Res. Lett. 48. doi:10.1029/2020GL091418

CrossRef Full Text | Google Scholar

Smith, M., Carrivick, J., and Quincey, D. (2016). Structure from Motion Photogrammetry in Physical Geography. Prog. Phys. Geogr. 40, 247. doi:10.1177/0309133315615805

CrossRef Full Text | Google Scholar

Soenen, S. A., Peddle, D. R., and Coburn, C. A. (2005). SCS+C: A Modified Sun-Canopy-Sensor Topographic Correction in Forested Terrain. IEEE Trans. Geosci. Remote Sensing 43, 2148–2159. doi:10.1109/TGRS.2005.852480

CrossRef Full Text | Google Scholar

Tarko, A., Bruin, S. D., and Bregt, A. K. (2018). Comparison of Manual and Automated Shadow Detection on Satellite Imagery for Agricultural Land Delineation. Int. J. Appl. Earth Observation Geoinformation 73, 493–502. doi:10.1016/j.jag.2018.07.020

CrossRef Full Text | Google Scholar

Teillet, P. M., Guindon, B., and Goodenough, D. G. (1982). On the Slope-Aspect Correction of Multispectral Scanner Data. Can. J. Remote Sensing 8, 84–106.

CrossRef Full Text | Google Scholar

Tolt, G., Shimoni, M., and Ahlberg, J. (2011). A Shadow Detection Method for Remote Sensing Images Using VHR Hyperspectral and LIDAR Data. Int. Geosci. Remote Sensing Symp. (Igarss), 4423–4426. doi:10.1109/IGARSS.2011.6050213

CrossRef Full Text | Google Scholar

van den Broeke, M., Smeets, P., Ettema, J., and Munneke, P. K. (2008). Surface Radiation Balance in the Ablation Zone of the West Greenland Ice Sheet. J. Geophys. Res. Atmospheres 113. doi:10.1029/2007JD009283

CrossRef Full Text | Google Scholar

van der Veen, C. J., Ahn, Y., Csatho, B. M., Mosley-Thompson, E., and Krabill, W. B. (2009). Surface Roughness over the Northern Half of the Greenland Ice Sheet from Airborne Laser Altimetry. J. Geophys. Res. 114, 264–275. doi:10.1029/2008JF001067

CrossRef Full Text | Google Scholar

Yang, K., Smith, L. C., Chu, V. W., Pitcher, L. H., Gleason, C. J., Rennermalm, A. K., et al. (2016). Fluvial Morphometry of Supraglacial River Networks on the Southwest Greenland Ice Sheet. GIScience & Remote Sensing 53, 459–482. doi:10.1080/15481603.2016.1162345

CrossRef Full Text | Google Scholar

Yang, K., Smith, L. C., Karlstrom, L., Cooper, M. G., Tedesco, M., As, D. V., et al. (2018). Supraglacial Meltwater Routing through Internally Drained Catchments on the Greenland Ice Sheet Surface. Cryosphere Discuss., 1–32.

Google Scholar

Yi, D., Zwally, H. J., and Sun, X. (2005). ICESat Measurement of Greenland Ice Sheet Surface Slope and Roughness. Ann. Glaciology 42, 83–89. doi:10.3189/172756405781812691

CrossRef Full Text | Google Scholar