Cast Shadow Detection to Quantify the Aerosol Optical Thickness for Atmospheric Correction of High Spatial Resolution Optical Imagery

Abstract

The atmospheric correction of optical remote sensing data requires the determination of aerosol and gas optical properties. A method is presented which allows the detection of the aerosol scattering effects from optical remote sensing data at spatial sampling intervals below 5 m in cloud-free situations from cast shadow pixels. The derived aerosol optical thickness distribution is used for improved atmospheric compensation. In a first step, a novel spectral cast shadow detection algorithm determines the shadow areas using spectral indices. Evaluation of the cast shadow masks shows an overall classification accuracy on an 80% level. Using the such derived shadow map, the ATCOR atmospheric compensation method is iteratively applied on the shadow areas in order to find the optimum aerosol amount. The aerosol optical thickness is found by analyzing the physical atmospheric correction of fully shaded pixels in comparison to directly illuminated areas. The shadow based aerosol optical thickness estimation method (SHAOT) is tested on airborne imaging spectroscopy data as well as on photogrammetric data. The accuracy of the reflectance values from atmospheric correction using the such derived aerosol optical thickness could be improved from 3–4% to a level of better than 2% in reflectance for the investigated test cases.

1. Introduction

Aerosol scattering is affecting the image quality of airborne and satellite imagery in the visible to near infrared range of the optical spectrum. The physical compensation of this effect requires knowledge of the total aerosol amount. The total aerosol load is typically quantified by the vertical aerosol optical thickness (AOT) or aerosol optical depth (AOD) as an integral value of the aerosol extinction in the vertical atmospheric column [1]. We use the term ‘AOT’ in this paper to describe the effective remote sensing aerosol optical thickness at a wavelength of 550 nm derived from airborne, downward looking system. This AOT is found by inverting the measurable aerosol scattering effect in the remote sensing imagery using radiative transfer calculations. On the other hand, the AOD is usually derived from ground based measurements in a different geometric situation.

Recent optical scanners provide reliable measurements of the at-sensor radiance at high spatial resolutions and with high sensitivity even in dark areas. Such data can be compensated for atmospheric effects by a physical inversion [2] using radiative transfer codes of the atmosphere such as MODTRAN5^® [3,4]. The atmospheric compensation may be done by well established software codes such as for instance ATCOR [5] or FLAASH [6]. The output of such an inversion is the bottom of atmosphere reflectance quantity, which can be approximated by the HDRF (hemispheric directional reflectance factor [7]) for directly illuminated pixels. For shaded pixels, the output is theoretically a bihemispherical reflectance (BHR), but, in the real world, shadow areas are often only partially corrected to a pseudo-reflectance due to missing information about the irradiance field. For the characterization of the surface, the goal should be to retrieve spectral albedo values (i.e., the bihemispherical reflectance, BHR) from each image pixel. This requires the correction of bidirectional reflectance distribution function (BRDF) effects [8] but also the physical correction of cast shadow areas by quantifying the diffuse irradiance. For these steps, the irradiance field is to be known, which is strongly related to the scattering by aerosols in the atmosphere.

Quantifying the aerosol amounts is an important step for fully automatic atmospheric compensation techniques as the AOD is used for both the calculation of the atmospherically backscattered radiance, as well as for the modeling of the direct and diffuse irradiance on the observed pixel. There are currently various approaches to get the aerosol scattering and haze distribution from the imagery: historically, the analysis of dark object signatures (DOS) has first been used for aerosol characterization [9,10]. The DOS method statistically searches for the darkest pixels in each image band and subtracts this value from the original data. This approach is strongly affected by the surface reflectance characteristics, which led to further evolutions where iterative approaches are used to improve the results [11]. Atmospheric correction is also a challenge for ocean color remote sensing, especially in coastal waters, where the assumption of ‘black’ water reflectance in the near infrared (NIR) is not fulfilled. Reinersman et al. 1998 [12] and Amin et al. 2014 [13] describe a method where the aerosol properties can be estimated from cloud shadow and sunlit areas over water assuming homogeneous areas.

The dark dense vegetation approach (DDV) considers vegetation reflectance properties for aerosol detection [14,15,16,17]. It uses the known correlation between the short wave infrared (SWIR) or the NIR band to the red band in vegetation to get an estimate of the expected red reflectance. The aerosol amount is then found as the observed offset in the red band. The relative accuracy of this method was found to be ≈$\pm 20\%$ for MODIS data [18]. This method meanwhile is a well established algorithm used in automatic atmospheric compensation routines [19]. Aerosols also have been successfully retrieved from multiangular satellite imagery [20,21] at comparable accuracies [22]; an approach which is not applicable to nadir geometry standard imagery. More recent developments involve simultaneous multidimensional optimization processes between surface reflectance and atmospheric parameters [23] or by taking into account external data sources like aerosol information from the GOME or the ATSR-2 sensor [24,25]. An approach leading to comparable accuracies in AOD (≈$\pm 0.07$ in AOT) is the use of multitemporal satellite observations [26] that requires a time series of bottom of atmosphere reflectance as a reference for the investigated site.

The goal of this work is to present an atmospheric compensation approach by means of a new operational aerosol optical thickness retrieval method for high spatial resolution optical remote sensing data, including airborne imaging spectroscopy. The method has been further evolved from first results presented by Schläpfer et al. 2017 [27]. The idea is to make use of the cast shadow areas within an image in order to analyze pixels illuminated by diffuse irradiance only. Approaches for aerosol retrieval from cast shadows in single high spatial resolution imagery have been investigated earlier by Thomas et al. 2013 [28]. Useful types of shadows for this kind of analysis are shadows of buildings and trees and terrain structures, whereas cloud shadows are not suitable due to the highly variable irradiance. Therefore, the method will only be applicable to data where a sufficiently high number of full cast shadow pixels are present, limiting the spatial resolution to approximately 5 m in typical remote sensing situations. The method tries to evaluate the spectral signature of the diffuse irradiance directly. Therefore, it does not require a specific ground reflectance type, in contrast to the DDV approach. Nevertheless, the high variability of shadow signatures [29] has to be taken into account by spatial averaging. High spatial resolution instruments can be supported by such a method as the AOT can potentially be retrieved from a single image only, avoiding the problems that would arise from co-registration issues for multiangular or multitemporal approaches.

The presented approach requires at first the accurate detection of cast shadows from within the image itself. The AOT can then be derived by estimating the diffuse irradiance term over the shaded pixels in an iterative way, using the modeling approaches implemented in atmospheric compensation routines. The combined shadow detection and aerosol retrieval method shall be suitable for 4-band photogrammetric data (with blue, green, red, and NIR-band) as well as for airborne imaging spectroscopy data. Even though the application of the method could potentially be extended to satellite imagery, this option is not analyzed here. Tests are performed accordingly on airborne optical data only in this paper.

2. Cast Shadow Detection Method

Cast shadow detection can be achieved by various ways, either from geometrical information or by spectral signature analysis, or by a combination of both [30,31]. For high resolution instruments and photographs, a wide range of shadow detection algorithms have been described so far [32], some of which are based on color transformation techniques [33].

Approaches based on geometric models usually involve the use of ray tracing procedure and highest accuracy terrain description. The use of digital surface models (DSM) has been found to be problematic for high spatial resolution instruments. This is not only due to the limited availability of suited models: for natural objects such as forests, the assumption of a continuous geometric surface does not hold true, which often results in erroneous mapping of shadow borders. With the use of LIDAR techniques, this may be solved in the future [34].

Matched filter and linear mixing procedures use the spectral statistics of the image to derive the shadow fraction [35]. Statistics are taken from non-shaded reflectance references to find large scale shadows. This method is well-suitable for cloud shadow detection [36], but its applicability is not given for high resolution situations such as settlements and tree shadows. The use of region growing approaches to get the correct shadow area borders [37] is a viable addition to these methods. Methods directly using the shadow properties are either based on thresholds to make use of the relative darkness of shadows [38,39], or based on the color properties of shadows [40,41]. Using a blackbody radiation model may further increase the reliability of shadow detection [42].

For the method presented herein, the focus is on image indexing using a minimum number of four spectral bands in order to retrieve the cast shadows in cloud-free conditions. The goal is to find a radiometrically stable method for detection of cast shadows over various natural surface covers, including water bodies, as well as over artificial surfaces such as asphalt. The at-sensor reflectance in four standard optical bands serves as an input. The bands are selected in the blue at 450 nm, the red at 550 nm, the green at 670 nm and the infrared at 780 nm, respectively. The method shall be applicable to single imagery directly in an operational and effective way, avoiding parameter adjustments based on single image statistics as far as possible. The method has been further developed from an initial concept presented by Schläpfer et al. 2013 [43]. It has to be noted that the signatures of shadows vary significantly depending on adjacent reflecting objects and the observed shaded pixels themselves [44,45]. These variations are to be minimized by taking the signatures of the darkest objects in an image into account.

The use of cast shadows for aerosol analysis is limited by the availability of fully shaded pixels in the image and thus by the spatial resolution. Assuming an average solar zenith angle of ${30}^{\circ }$ and typical heights of buildings and trees between 10 and 30 m, the maximum spatial sampling interval for this method is theoretically between 6 and 17 m per pixel. Considering the spatial blurring due to the optical point spread function (PSF) and the variability of the solar zenith angle, a maximum sampling interval of 5 m is considered being the appropriate upper limit to apply this method. This limit may be larger in special cases only if, e.g., spatially distributed terrain shadows are available in an image.

2.1. Cast Shadows over Land Surfaces

Our analysis suggests that land and water surfaces have to be treated differently for shadow detection. The detection of shadows over natural and man-made land surfaces is to be treated by one single index. First, the data is normalized to an apparent reflectance value from calibrated at-sensor radiance as:

$${\rho }_{app}=\frac{\pi L_s}{E_{0}cos\left({\theta }_o\right)},$$

(1)

where $L_s$ is the at-sensor radiance, $E_0$ is the top of atmosphere irradiance and ${\theta }_o$ is the solar zenith angle. A simulation of at-sensor reflectance for typical objects in the visible/near infrared (VNIR) spectral range has been performed from a spectral library of a total of 85 spectra. It contains a subset of the ASTER spectral library [46] with human-made and soil targets and simulated vegetation spectra for Leaf Area Index (LAI) = 1 and LAI = 4 and three levels of chlorophyll a + b, from PROSPECT/SAIL simulations [47]. At-sensor radiance simulations were done using MODTRAN5 for sensor altitudes of 2 km, 5 km and top of atmosphere and for a variable visibility at levels of 5, 10, 20, 40, and 80 km (corresponding to an AOT of 0.01 to 0.78). Visibility is used as a proxy for aerosol optical thickness for practical reasons. The direct illumination was scaled with a shadow fraction from $f_{sh}=0\%$ to 100% in steps of 1%, where $f_{sh}=100\%$ is fully illuminated. The water vapor amount was set to 2 cm columnar water vapor, the solar zenith was set to ${30}^{\circ }$ and the standard rural aerosol model was used. The simulation follows a simplified forward formulation of at-sensor radiance for clear sky conditions:

$$L_s=L_{atm}+\frac{\rho \left(E_{dir,t}cos\left({\theta }_o\right)f_{sh}+E_{dif,t}\right)}{\pi (1-s\rho )},$$

(2)

where $L_{atm}$ is the atmospheric back scattered path radiance, $E_{dir,t}$ and $E_{dif,t}$ are the direct and diffuse irradiances from the ground transmitted to the sensor, $f_{sh}$ is the shadow fraction defined from 0 to 1 with a value of 0 for cast shadow areas, s is the spherical albedo, and $\rho$ is the surface reflectance. Parameters other than $f_{sh}$ and $\rho$ are taken from the ATCOR look-up tables, which have been compiled using MODTRAN 5.4 [4].

Two typical samples of simulated spectra are shown in Figure 1 for concrete and vegetation. The apparent reflectance values increase towards the blue spectral range as the path scattered radiance is increased and the bottom of atmosphere reflectances are constant or do typically slightly decrease towards the blue. A similar shape can be found for most cover types of the library except for water and blue man-made materials. After normalizing the apparent reflectance to the red spectral band (right graphs of Figure 1), it can be shown that the relative increase of the observed at-sensor reflectance in the blue is more pronounced for shaded areas than for directly illuminated surfaces. This ratio between the blue and the red spectral band shall be exploited for shadow detection. Furthermore, it has to be considered that vegetation signals on average show higher ratio values between the blue and the red bands than spectrally flat surfaces such as soils due to the strong chlorophyll absorption in the red spectral range by plants.

The shadow fraction $f_{sh}$ parameter is now to be calculated from this red-to-blue spectral feature. The respective apparent reflectance values are used as input. They are simply named ${\rho }_b,{\rho }_g,{\rho }_r$, and ${\rho }_n$ for the blue, green, red, and NIR spectral bands, respectively. In a first step, the ratio between the red and the blue spectral band is calculated. This index increases with the portion of direct illumination: low index values show shadow areas, whereas high values are directly illuminated pixels. The blue apparent reflectance is corrected to account for the influence of chlorophyll absorption, using the positive difference in apparent reflectance between NIR and red ${({\rho }_n-{\rho }_r)}_{>0}$. Setting the minimum of this correction factor to zero inhibits its unintended application to non-vegetated pixels:

$$i_{sh}={\displaystyle \frac{{\rho }_r+k_n{\left({\rho }_n-{\rho }_r\right)}_{>0}}{{\rho }_b}}.$$

(3)

The factor $k_n$ is used to weight the vegetation correction factor and is to be found empirically. It is derived by minimizing the variation of $i_{sh}$ as of Equation (3) depending on surface cover types, using the spectra of the above-mentioned simulated data set. An optimal value of $k_n=0.1$ has been found from this analysis.

When calculating the index from Equation (3) for various aerosol amounts, a systematic decrease of the index value with increasing aerosol amounts is observed (see Figure 2, left). This can be accounted for by using dark object signatures in the blue, which are mainly ruled by aerosol scattering (see dotted lines in Figure 2). The dark signature ${\rho }_{b,dark}$ for the blue band is found from the image by using the average apparent reflectance of the darkest pixels in the image from 0.1% of all image pixels (this percentile is increased to 1% for small images). Relating these dark object signatures, an aerosol-amount-dependent trend can be found and applied to the index as a normalization function such that the shading index for land surfaces may be written as:

$$i_{lsh}={\displaystyle \frac{{\rho }_r+k_n{\left({\rho }_n-{\rho }_r\right)}_{>0}}{{\rho }_b}}{\displaystyle \frac{1}{a{e}^{b{\rho }_{b,dark}}}}.$$

(4)

The parameters a and b from the correction function are found by an exponential fit between ${\rho }_{b,dark}$ and the raw index $i_{sh}$ on all spectra of the simulated data set, on three sensor altitudes (2 km, 5 km, and satellite level). Values of $a=1.58$ and $b=-0.04$ have been found leading to the lowest variation of $i_{lsh}$ with respect to the visibility. The relative deviation from mean of the index derived from the full set of reference spectra as described above with the full variation of visibilities could be reduced from 32% to 14% by applying this correction. Applying Equation (4) with optimized correction parameters to the simulated data leads to a more stable average index as shown in Figure 2, right. However, some dependency on the flight altitude remains. This index can now be used for shadow detection by defining appropriate thresholds as described in Section 2.3 below.

2.2. Cast Shadows over Water Bodies

The above index is applicable for most land situations, but it can not be used to derive shadows over water bodies. Specific methods have to be developed to find shaded areas over water in order to exclude false alarms in water. An automatic approach for shadow and water delineation for imaging spectroscopy data was shown by Bochow et al. [48], and the feasibility for cloud shadow detection in a marine environment has been reported in earlier studies [13,49]. As our method is restricted to four red-green-blue (RGB) and NIR spectral bands, we propose a method to quantify the cast shadows over water by using the same irradiance effect as described in Section 2.1. However, using the red-blue ratio does not allow for a good shadow detection over water. The blue appearance of water reflectances is interfering with the aerosol signatures in this wavelength range and needs to be accounted for. It has been found that an extrapolation of the signatures in the green-red bands towards the blue spectral range is a good alternative to find a shadow-sensitive index for water bodies. A corresponding index is calculated as a relation between an extrapolated blue signature to the measured blue apparent reflectance, while both parameters are corrected by subtracting the blue dark objects signature:

$$i_{wsh,0}={\displaystyle \frac{{\rho }_{b,ext}-{\rho }_{b,dark,ext}}{{\rho }_b-{\rho }_{b,dark,ext}}},$$

(5)

where

$${\rho }_{b,ext}={\rho }_g-{\displaystyle \frac{{\rho }_r-{\rho }_g}{{\lambda }_r-{\lambda }_g}}({\lambda }_g-{\lambda }_b).$$

(6)

The dark object signature for the blue band can not be taken from the blue band directly as the natural blue reflectance of water may be quite high. Therefore, this signature ${\rho }_{b,dark,ext}$ is estimated as a weighted combination between the blue dark signature and the red dark signature as:

$${\rho }_{b,dark,ext}=0.8{\rho }_{r,dark}+0.2{\rho }_{b,dark}.$$

(7)

Results of this normalization to the red-green trend line on two typical water surfaces are shown in Figure 3. A good discrimination of shade levels in the blue spectral range can be achieved (compare also the water area in Figure 4). For this index, the aerosol influence also needs to be accounted for, but in a different way than for land surfaces. It has been found that the index depends mainly on flight altitude rather than on aerosol optical thickness, which may be explained by the relation between atmospheric path radiance and diffuse ground irradiance that changes with altitude. An exponential decrease of the index value with flight altitude h has been derived from the simulated data set, which can be used for a height-dependent index equation as:

$$i_{wsh}={\displaystyle \frac{i_{wsh,0}}{a{e}^{-bh}}}.$$

(8)

The values of the exponential function are derived from the simulated data as $a=1.18$ and $b=-0.4{\mathrm{km}}^{-1}$. The flight altitude h is given in km.

2.3. Index Combination

The indices can be combined based on a land/water mask. The water mask is derived by a classification rule that uses the following threshold values:

$$\begin{array}{l}\left[({\rho }_n<({\rho }_{nir,dark}+1\%)\mathrm{AND}({\rho }_r-{\rho }_b<-3\%)\mathrm{AND}({\rho }_n-{\rho }_r<3\%)\right]\\ \mathrm{OR}\left[(2\times {\rho }_g>({\rho }_r+{\rho }_n)+4\%)\mathrm{AND}({\rho }_n<7\%)\right].\end{array}$$

Using a transition range of 1% for each reflectance threshold, a smooth weighting factor w between water ($w=0$) and land surfaces ($w=1$) is derived. Offsets are applied to both indices to place them in a range between 0 and 1 such that $i_{lsh}^{\prime }=(i_{lsh}-0.3)>0<1$ and $i_{wsh}^{\prime }=(i_{wsh}-0.6)>0$. The two re-scaled indices are then combined into one layer for later use in image analysis:

$$i_{sh}=w{i}_{lsh}^{\prime }+(1-w)i_{wsh}^{\prime }.$$

(9)

The parameter $i_{sh}$ is finally scaled to a shadow fraction $f_{sh}$ between 0 (full cast shadow) and 1 (no cast shadow) using a lower limit for full cast shadow areas and an upper limit for full illumination. It has been found that the lower limit of $i_{sh}$ has to be set to values between 0.3 and 0.55. This threshold depends on the instrument spectral response, the radiometric calibration accuracy of the instrument, but also on the average flight altitude. In desert areas where no water is expected, the detection of shadows on water is not useful and may lead to erroneous results. Therefore, a simplified mode without the second index is applied in such situations by setting the parameter w to one. A sample of the whole cast shadow detection process is illustrated in Figure 4. The indices for land and water are calculated and then combined into one fractional shadow value using the water mask. Application of a final threshold leads to the cast shadow mask.

For topographic correction in the ATCOR software, the direct illumination on an image pixel is scaled by the cosine of the local solar incidence angle cos($\beta$), which is typically derived from digital elevation data. As the image derived shadow fraction $f_{sh}$ is also related to the direct illumination, it can be combined with the incidence angle in an illumination factor $f_{ilu}$ as the smaller value of both factors per image pixel:

$$f_{ilu}=min\left[f_{sh};cos\left(\beta \right)\right].$$

(10)

This fractional illumination factor may now be used for correction of illumination effects in the atmospheric compensation routine, for cast shadow correction, but also for aerosol detection as described hereafter.

3. Aerosol Signature Retrieval

Shaded areas are mainly illuminated by the diffuse irradiance, which strongly depends on the aerosol distribution and multiple scattering effects. The irradiance on directly illuminated surfaces, on the other hand, is driven by a combination of aerosol absorption and scattering effects, often cancelling each other. We propose a method that makes use of this difference in the irradiance description to derive the amount of aerosols from imagery containing a decent portion of fully shaded areas. Due to high variability of shadow signatures [29,44], the method is evaluating average aerosol contents over image patches. Note that there is a strong dependency between water vapor retrieval and aerosol detection, which should be considered if detecting water vapor over shaded areas [50,51]. For this reason, this method uses image bands with negligible water vapor absorption.

3.1. Theoretical Background

The shadow based aerosol optical thickness retrieval (SHAOT) is performed on an image patch of selected size between 50 m and 2 km diameter, assuming constant aerosol contents within the patch. The method uses an atmospheric compensation model and minimizes the difference between the retrieved reflectance values in cast shadow areas to adjacent directly illuminated areas. The reflectance retrieval formulation is using the concepts and formulations as used in the ATCOR model [52]. The model includes a Look-Up-Table (LUT) of all relevant parameters required for atmospheric correction, pre-calculated using MODTRAN 5.4. Using these parameters and the acquisition geometry, the bottom of atmosphere $\rho (x,y)$ is derived from the calibrated at-sensor radiance $L_s$ as:

$$\rho (x,y)=\frac{\pi \left[\left(d^2{L}_s-L_p({\theta }_v,{\theta }_o)\right)(1+q)-q\overline{L_{s,dir}}\right]}{E_{dir,t}{f}_{ilu}+E_{dif}+E_{ter}},$$

(11)

where the path radiance term $L_p=L_{p,single}+L_{p,mult}$ is derived from the ATCOR LUT depending on ground altitude z, sensor altitude h, observation zenith angle ${\theta }_v$, solar zenith angle ${\theta }_o$, and the relative azimuth angle $\varphi$. The factor d is the relative sun–earth distance. The adjacency radiance $\overline{L_{s,dir}}$ has to be considered in this retrieval as the aerosol scattering is the key driver to the signal over dark areas. The adjacency factor $q={\tau }_{dif}/{\tau }_{dir}$, i.e., the relation between diffuse and direct transmittance is used for this purpose in Equation (11). The parameters ${\tau }_{dif}$ and ${\tau }_{dir}$ are also read from the ATCOR LUT. The adjacency factor is applied in the numerator of the equation in one single iteration by using the averaged first order direct at-sensor radiance of the image patch $\overline{L_{s,dir}}=\overline{d^2{L}_s-L_p}$ [53]. Note that the parameters q and $E_{dif}$ have significant impacts on the modeled directly reflected radiance in shaded pixels. Therefore, these parameters can not easily be approximated or neglected in the calculation.

The total irradiance is calculated using direct irradiance $E_{dir,t}$ and the diffuse irradiance $E_{dif,t}$ transmitted to the sensor level for a 100% reflecting target. These values are also available in the ATCOR LUT. The factor $f_{ilu}$ is the combined illumination and shadow fraction as described in Equation (10). The diffuse term $E_{dif}$ is defined according to [54] as:

$$E_{dif}=E_{dif,t}\left[{\tau }_s\frac{f_{ilu}}{cos\left({\theta }_0\right)}+(1-{\tau }_s)V_{sky}\right],$$

(12)

where $E_{dif,t}$ is the generalized diffuse irradiance on the ground transmitted to the sensor level, ${\tau }_s$ is the sun-to-surface transmittance, and $V_{sky}$ is the skyview factor for the given pixel in terrain (i.e., the fraction of the visible sky hemisphere as seen from a pixel). Furthermore, the irradiance from the adjacent terrain on a pixel $E_{ter}$ is described as:

$$E_{ter}=\left(E_{dir,t}cos{\theta }_o+E_{dif,t}\right)\left(1-V_{sky}\right)\overline{{\rho }_{ter}},$$

(13)

where $\overline{{\rho }_{ter}}$ is approximated as the spatially averaged reflectance of all pixels in the image patch after applying Equation (11) with $E_{ter}=0$ as a first order correction.

Note that effects of gaseous absorption can be neglected in these formulations as only spectral bands in atmospheric windows are used, avoiding interference of absorption with aerosol scattering. Thus, the irradiance for a shaded pixel with $f_{ilu}=0$ is given as:

$$E_{shade}=E_{dif}+E_{ter}=E_{dif,t}\left(1-{\tau }_s\right)V_{sky}+\overline{{\rho }_{ter}}\left(E_{dir,t}cos{\theta }_o+E_{dif,t}\right)\left(1-V_{sky}\right).$$

(14)

This equation can be used as irradiance term in Equation (11) for reflectance retrieval in cast shadow areas.

The effect of variable aerosol amounts on the atmospheric correction of cast shadow areas is depicted in Figure 5. It shows the results of atmospheric correction including cast shadow removal based on the displayed cast shadow mask and using Equation (11). The results vary from under-corrected shadows for large aerosol amounts to overcorrection at lower aerosol amounts. An optimum visibility parameter can be found where the corrected cast shadow brightness is in agreement with the brightness of the directly illuminated areas. Thus, the atmospheric correction may be iterated with the visibility parameter until the brightness of the corrected pixels agrees to a reference of directly illuminated pixels.

We use the visibility parameter as defined in MODTRAN for this retrieval for practical reasons. The visibility parameter is transformed into vertical AOT using the MODTRAN-based relationship between vertical optical thickness and visibility (VIS) :

$$AOT(VIS,z)=e^{a\left(z\right)+b\left(z\right)ln\left(VIS\right)},$$

(15)

where a(z), b(z) are coefficients obtained from a linear regression of $ln\left(AOT\right)$ versus $ln\left(VIS\right)$ on the various height levels z [52].

3.2. AOT Retrieval Procedure on Image Patch

The following method is implemented within the framework of the ATCOR-4 atmospheric compensation package [52]. The procedure works iteratively on image patches of typical sizes between 200 × 200 and 2000 × 2000 pixels and spatial diameters between approximately 50 m to 2 km. The visibility optimization is started for a single image patch if at least 100 pixels of bright reference and 300 pixels of shaded area are found within the patch. The larger number of dark pixels is required to improve the statistical reliability of the highly sensitive dark pixel information while computationally intensive single pixel modeling [55] can be avoided. Within this image patch, the MODTRAN-based visibility parameter is used as a tuning factor for the average aerosol amount. For each image patch, the following steps are performed (compare also Figure 6):

• Create the data subset patch, containing at-sensor radiance image, shadow fraction, skyview factor and subset of atmospheric LUT for average ground altitude.
• Derive patch-specific information such as the ground altitude and the view zenith angle.
• Find masks for cast shadow pixels and bright neighbor reference pixels (compare also last image in Figure 7 and see comment below).
• Set the start values for the iteration, starting with a visibility of 80 km.
• Find the optimal visibility by applying the atmospheric compensation including shadow correction iteratively:

  (a)

  Get atmospheric parameters and adjacency weighting factor for given visibility.

  (b)

  Calculate direct and diffuse irradiance.

  (c)

  find at-sensor radiance corrected for path radiance while correcting the adjacency influence.

  (d)

  Find terrain irradiance using average first estimate of neighboring reflectance.

  (e)

  Calculate reflectance value for whole image patch.

  (f)

  Compare brightness of pixels in cast shadow areas to illuminated bright areas, and

  (g)

  stop iteration if the average difference between reference pixels and shadow corrected pixels drops below 0.05% or after 30 steps.

• Store the derived optimal visibility value and its location within the image.

The pixels used as a reference for the iteration are found by shifting the mask of the cast shadow areas in direction of the local solar azimuth angle by a number of pixels, which is given in dependency of the spatial resolution $s_{pix}$ in [m] as $n_{pix}=min[(20-s_{pix});6]$, i.e., a number of 20 pixels for resolutions below 1 m and a minimum number of six pixels for lower resolutions. Pixels that are shaded by more than 50% are excluded from the list of reference pixels. The output of the retrieval on an image patch is one single visibility value, which is later transformed into an AOT value using Equation (15).

3.3. Application to Large Scale Imagery

The application of the method to large scale imagery is done by a moving window technique using the window sizes as described in the previous section. The complete AOT retrieval scheme is shown in Figure 6. It uses the outputs from the shadow detection procedure as described above including a skyview factor map, which takes the reduced skyview factor within cast shadow areas into account. The aerosol detection then follows the below steps:

• Read meta data and calibration data for the sensor.
• Load spectral band closest to a predefined reference wavelength (typically at 550 nm) and transform to at-sensor radiance.
• Load scan angle, shadow fraction, and skyview file for current image.
• Loop over moving window using a the pre-defined tile size.
• Load specific atmospheric LUT for flight altitude of instrument, solar zenith angle and average ground altitude and reference band for the tile.
• Calculate visibilities for the tiles (i.e., for the image patches) using the procedure described in the previous section.
• Omit outliers, triangulate valid retrievals, and interpolate visibility values in order to cover the full image extent.
• Derive the height-dependent AOT map using the digital elevation model, the visibility index, and the MODTRAN-based relation between aerosol optical thickness and visibility according to Equation (15).

The result of this procedure is a map of the observed vertical aerosol optical thickness at the reference wavelength (i.e., at 550 nm), which later can be used as an input to the atmospheric compensation step.

3.4. Sensitivity Analysis

The sensitivity of the aerosol retrieval on various parameters has been tested by analyzing the retrieved AOT values from two image subsets of the ADS system (compare Section 4.3 below). The subsets have been selected in the center of a settlement and in a sparsely forested area. The impact of six critical parameters on the AOT retrieval results is summarized in

Table 1. Sensitivity of aerosol retrieval on various error sources for 0.25 m resolution data, given as percentage error in aerosol optical thickness (AOT). The (+) sign is added whenever systematic increases are observed, whereas the other values are statistical variation levels.

Error Source	Forest	Settlement	Remarks
overestimate cast shadow	(+) 13.5%	(+) 51.7%	includes directly illuminated areas
underestimate cast shadow	5.8%	6.4%	no clear trend
vary reference shift distance	0.4%	1.8%	in steps of 1 m
vary reference angle	6.0%	0.9%	in steps of ${90}^{\circ }$
neglect skyview factor	(+) 5.2%	(+) 5.3%	set to 100%
neglect adjacency effect	(+) 30.0%	(+) 52.4%	no adjacency correction

. The cast shadow detection is shown to be crucial; if the cast shadow areas are too large, the aerosol amount is significantly overestimated. This is due to the fact that bright areas are included in the shadow area and lead to a distortion of the statistics. Too small cast shadow areas, on the other hand, may lead to smaller errors on a 6% level. An error level of up to 6% is also found if analyzing the impact of the selection of the reference reflectance areas. This has been tested by varying the distance at which the reference areas are searched and by using angles differing from the solar azimuth for reference area calculation. The latter introduces larger errors for forest areas. This may be explained by the reflectance difference of trees in comparison to the forest undergrowth: if staying in the direction of the shadow, the reference is more likely of the same type as if using at a differing angle. A last test is done by omitting radiometric elements in the retrieval: omitting the skyview factor influence leads to an increased aerosol amount by about 5% and omitting the adjacency correction leads to aerosol AOT overestimations at a level of over 30%. It can be concluded that the correct treatment of the adjacency effect and the accurate shadow detection are the most important factors for the retrieval accuracy.

4. Validation Results

Validation of this method is shown on two types of imagery: the airborne photogrammetric camera system Leica Airborne Digital System (ADS) [56] and the Airborne Prism Experiment (APEX) imaging spectrometer [57]. Further validation results on Hyspex imagery have already been published earlier [27]. The applicability of the method on satellite data has been tested, but the results are still under investigation.

The method has been implemented in the IDL (interactive data language) programming language and is embedded in the ATCOR atmospheric compensation software, which is commercially available. Hereafter, some results for shadow detection, aerosol retrieval and subsequent atmospheric correction are shown. For each instrument, the thresholds for shadow detection have to be determined prior to the processing. They have been set to a value of 0.33 for the APEX imager and to 0.51 for the ADS system (compare discussion in Section 5). All further processing is done fully automatically and is performed in an operational way for a series of images.

4.1. Validation of Cast Shadow Detection

The cast shadow detection method is validated for both data sets using subsamples of 120 × 120 pixels each (details regarding the data sets can be found in the subsequent sections). A shadow mask is interpreted manually on the two subsets. This mask is used as a reference for shadow detection. For comparison purposes, shadow detection is performed by four methods for each sample: a simple near infrared threshold, the scaled image brightness according to Nagao (1979) [38], the Hue-Saturation-Intensity (HSI) conversion method according to Tsai [33], and the herein presented method. Spatial results of the methods are shown in Figure 8 and Figure 9. Shadow detection results are comparable between the methods. The simple NIR threshold is showing most errors as reported in earlier studies [32] and the Nagao method improves the results. Both the HSI and the SHAOT visually agree best with the reference map.

A quantitative comparison was done by calculating the classification accuracy by the kappa value between the reference ‘shade’ to ‘no shade’ classification and the results of the four methods. The thresholds for all methods have been optimized manually in order to achieve maximum kappa values. Results of this comparison are shown in

Table 2. Classification accuracy (kappa values) on two sample images for four various shadow detection methods.

Test Case	NIR Threshold	Nagao Method	HSI Inversion	SHAOT Method
ADS (airborne digital system)	69%	82%	71%	85%
APEX (airborne prism experiment)	46%	62%	75%	76%

NIR: near infrared; HSI: Hue-Saturation-Intensity

. For the ADS sample, the SHAOT method shows a better overall agreement than the other methods, whereas, for the APEX case, the classification based on the HSI inversion and the SHAOT method are on the same level of accuracy.

The validation shows a performance of the herein presented method, which is well comparable to established shadow detection techniques, and may even be on a higher accuracy level. More validation work would be required to show the absolute accuracy of the retrieval method. It is to be remembered, however, that this shadow detection method is meant as a stable shadow detection method to be used for aerosol scattering retrieval. The major difference to the other tested methods is the fact that this method uses radiometric indices based on apparent reflectance data that allows a higher stability and less dependency on the individual image statistics of the resulting shadow distribution.

4.2. Results on APEX Imaging Spectroscopy Data

The APEX imaging spectrometer provides data in 299 contiguous spectral bands covering the spectral range from 380 to 2500 nm, with an useful range between approx. 450 and 2350 nm [58]. The investigated data sets were acquired on 26 June and 29 June 2010, over the Laegern test site in Switzerland. For 26 June, the acquisition time was 2:59 p.m. UTC, at a flight altitude of 4.428 km a.s.l. and a solar zenith angle of $48.1^{\circ }$. For 29 June, the acquisition time was 10:07 a.m. UTC, flight altitude at 4.470 km a.s.l. and the solar zenith angle was $29.3^{\circ }$. Data were processed to at-sensor radiance within the APEX Processing and Archiving Facility [59] using laboratory calibration coefficients [60]. Validation is done by analysis of the bottom of atmosphere reflectances in comparison to ground reference spectra. Furthermore, some of the AOT results could be compared to AOD values retrieved from Aeronet measurements [61].

The shadow detection routine successfully retrieved the shadows using a lower threshold value of 0.33 of the combined index as of Equation (9). The results are shown in Figure 10. The shown fractional illumination $f_{ilu}$ is the terrain shading as derived from a smoothed digital surface model (DSM) combined with the image derived fractional cast shadow. Cast shadows from building and trees are well detected and also most shades over water areas can be discriminated. Problems occur in turbid waters and with blue artificial objects, which may be wrongly detected as shadows by the routine.

The absolute aerosol optical thickness is evaluated in comparison with in situ aerosol optical depth measurements. The Aeronet station ‘Laegern’ is located in the center of the given scene at a ground altitude of 735 m. It delivered data for both dates: on 26 June, the Aeronet aerosol optical depth at 550 nm was at $0.21\pm 0.02$. On 29 June, the value of the AOD was at $0.51\pm 0.05$. The AOT from imagery was at 0.25 (SHAOT) and 0.19 (DDV) for the first day and AOT = 0.49 (SHAOT) and 0.39 (DDV) for the second day. The agreement of the SHAOT results is within the expected accuracy of AOT retrieval, whereas the DDV method underestimates the aerosol amount for the second day. Note that, without applying the term of adjacency correction as described in Equation (11), an overestimation of the AOT by about 20% was observed at first. Also be aware that we compare sun photometer based AOD retrievals to image based AOT, which are acquired in different geometric situations such that the natural variability within a guessed range of 1 km may lead to significant differences.

Using the retrieved aerosol optical thickness distribution for atmospheric compensation leads to correction results as displayed in Figure 11. The outputs are calculated in raw geometry while making full use of the digital elevation data with the topographic atmospheric correction module in ATCOR-4. Visual comparison shows under-correction artifacts over the forests in the reflectance outputs of the DDV based correction, whereas the SHAOT results do not exhibit this problem.

For spectral evaluation, ground reference spectra are compared to the bottom of atmosphere reflectance outputs. The spectra have been acquired in-field on the day of overflight using the ASD field spectrometer, at a solar zenith angle ranging from $34.9^{\circ }$ to $44.2^{\circ }$ as described in [62]. In Figure 12, the difference between the DDV and the SHAOT based atmospheric compensation is depicted. The aerosol amounts found from shadow analysis are slightly higher than estimated from DDV, specifically for the 29 of June. This leads to lower reflectance values for dark objects in the blue spectral range and a better agreement to most ground references. For bright objects, only small differences are visible in the blue as aerosol scattering and absorption cancel each other at high reflectance levels. For dark spectra, some differences appear in the near infrared, where the adjacency effects from the dominant vegetation areas are insufficiently corrected with the DDV results, whereas the correction is satisfactory with the SHAOT amount. Furthermore, an offset of 2–3% in reflectance can be observed for dark objects in the DDV based results. Remaining differences between the two acquisition dates may be attributed to incidence BRDF effects due to the variable solar zenith angle between $29.3^{\circ }$ and $48.1^{\circ }$. For targets with little BRDF variations such as asphalt and dark roof, the spectra agree well between acquisition dates, specifically for the SHAOT results. Note that the sand spectra was acquired on a sports field with disturbed sand surface, which is affected by self shading effects and therefore exhibits incidence BRDF variations. Correlation analysis shows an improved average correlation coefficient between ground reference spectra and image spectra from 0.895 for the DDV correction to 0.923 for the SHAOT based correction. Remaining differences in the blue spectral range are still to be discussed, as they may be related to sensor calibration issues.

4.3. Results on Airborne Digital Photogrammetry Data

The Leica ADS camera scans the earth in four spectral bands at a total FOV angle of ${46}^{\circ }$ at spatial resolutions up to 10 cm using the latest ADS-100 cameras. For large scale atmospheric correction of such data, aerosol detection and correction is of interest, specifically in rugged terrain where aerosol amounts often vary significantly due to local meteorological conditions. The Swiss Federal Institute of Topography (swisstopo) uses ADS-80 and ADS-100 systems for regular mapping applications in Switzerland. Atmospheric compensation has been applied to the data to increase the spectral consistency in mountainous areas and to make multitemporal data better comparable. Experiences on these data showed that the DDV approach does not provide reliable aerosol distribution maps for such high resolutions; the approach failed in some situations. This led to the development of the herein presented approach. Due to the high spatial resolution of the ADS system, cast shadows can be found in almost all data acquisition situations at mid-latitude illumination conditions. Figure 7 shows a sample of ADS processing, with shadow detection results on top and the used aerosol retrieval masks below. It also illustrates the reference pixel areas that have been used for comparison between cast shadow correction results and the bright reference areas while doing the iteration on the aerosol amounts. The tentative cast shadow correction result is showing well fitting corrections within the settlement, whereas overcorrections in some forest areas can be observed. As the main focus of this paper is the aerosol retrieval, we do not follow up on the shadow correction any further at this point, although the correct removal process is an important further research topic [63,64].

For evaluation of the method in a critical situation, a large scale sample of ADS-80 data has been selected consisting of three scenes that were acquired on 19 March 2015 over the Swiss Jurassic mountains. The solar zenith angle was at $48.9^{\circ }$. The data was flown at an altitude of 2.82 km a.s.l. with an ADS-80 sensor at a spatial resolution of 0.25 m and is covering an area of 23.8 × 2.1 km, resulting in a data size of 6.4 GB for the four spectral band data. Fully automatic atmospheric compensation was run with both, the DDV and with the SHAOT method. The standard atmospheric processing takes 56 min on a 2.6 GHZ Intel Quad-Core i7 processor system with solid state disks, while the processing time increases to 98 min for the shadow based processing. The standard atmospheric correction uses a default visibility of 50 km with a rural aerosol model and takes the large scale terrain effects into account, including adjacency effects and terrain irradiance effects. The results are shown in Figure 13. Aerosol scattering effects are mostly removed by the SHAOT method while some blurring remains for the standard atmospheric correction with an automatically chosen visibility of 50 km (last image). As the data was acquired before the vegetation growing season, dark vegetation pixels are sparse and the DDV method for aerosol retrieval fails for most of the image tiles; i.e., the DDV correction results are worse than without atmospheric compensation and the three image tiles are inconsistently corrected (not shown in figure). Therefore, the standard visibility of 50 km had to be chosen as fall-back solution for the standard processing.

The SHAOT method leads to reasonable AOT results for this area ranging between 0.2 and 0.6. The retrieved AOT are visibly higher over the lake, which may be related to the difficult discrimination between water and shadow pixels in the turbid lake. In addition, some unnatural spatial AOT variations are visible that are attributed to the method inherent statistical error. However, using this aerosol content allows a successful atmospheric compensation over land (compare Figure 13). Aerosol variability from partially snowy mountain areas at altitudes of 1300 m a.s.l. down to the city area at an altitude of 300 m a.s.l. could be mostly removed. No absolute validation of AOT is possible for this data as no Aeronet data is available for the time and location of data acquisition. However, the visual comparison between the SHAOT result and the standard atmospheric correction clearly shows an improved correction of the areas at lower ground altitudes.

The spectral variability is analyzed by comparing the statistics of four image patches of about one square kilometer each, covering similar agricultural surfaces types (red frames in Figure 7). The statistical results for the subsets and the two spectral bands that are mostly affected by aerosol scattering in the blue and the NIR are shown in

Table 3. In-scene consistency analysis between two image patches on the blue and the NIR band of ADS imagery, comparing apparent at-sensor reflectance values (${\rho }_b$ and ${\rho }_n$) to the standard reflectance products (${\rho }_{b,std}$ and ${\rho }_{n,std}$) and the dark dense vegetation (DDV) and SHAOT based reflectance outputs.

	Altitude [m]	${\mathbf{AOT}}_{\mathbf{DDV}}$	${\mathbf{AOT}}_{\mathbf{SHAOT}}$	${\mathit{\rho }}_{\mathit{b}}$	${\mathit{\rho }}_{\mathit{n}}$	${\mathit{\rho }}_{\mathit{b},\mathbf{std}}$	${\mathit{\rho }}_{\mathit{n},\mathbf{std}}$	${\mathit{\rho }}_{\mathit{b},\mathbf{ddv}}$	${\mathit{\rho }}_{\mathit{n},\mathbf{ddv}}$	${\mathit{\rho }}_{\mathit{b},\mathbf{shaot}}$	${\mathit{\rho }}_{\mathit{n},\mathbf{shaot}}$
Sample 1	705	0.51	0.36	7.8	23.4	6.7	24.7	5.3	25.7	5.0	25.2
Sample 2	609	0.14	0.44	8.7	28.0	7.7	29.6	3.9	29.3	5.3	30.7
Sample 3	524	0.14	0.44	8.8	26.6	7.7	28.2	3.8	27.9	5.2	29.2
Sample 4	425	0.16	0.50	9.6	22.4	8.6	23.7	8.6	23.5	5.6	24.5
Mean				8.7	25.1	7.7	26.6	5.4	26.6	5.3	27.4
Mean Deviation				0.5	2.2	0.5	2.4	1.6	2.0	0.2	2.6

. The atmospheric compensation based on the SHAOT-derived AOT values leads to a higher consistency between the image samples in the blue spectral range, where the natural variation of the average reflectance is knowingly low for land surfaces. On the other hand, the contrast is enhanced in the NIR. This can be explained with the natural variation of near infrared vegetation signatures between the samples at various ground altitudes. Such natural variations are more pronounced after atmospheric compensation as blurring effects by the atmosphere are removed. The DDV based results lead to the worst agreement due to the above-mentioned reasons. Note that, due to failure of the DDV method, default values for aerosol amounts have been chosen automatically for samples 2 and 3, leading to the large variation.

Incomplete correction of aerosol scattering mostly affects the reflectance outputs of atmospheric correction on dark objects. Therefore, the reflectance outputs over asphalt samples (3 × 3 pixels) and the lake (20 × 20 pixels) are compared to reference spectra from the John Hopkins University (JHU) / Aster spectral library for the asphalt and to a typical oligotrophic lake spectrum provided with the ATCOR code. The locations of the samples are indicated in Figure 13 as yellow circles for the lake and in green for the asphalt. The results of this validation are shown in Figure 14. The levels of reflectance for the asphalt target is at 8.0% for the reference. With the SHAOT method, an average value of 9.1% is achieved, whereas the standard correction leads to a level of 10.9% reflectance. Furthermore, the typically increasing reflectance with wavelength is only visible after the SHAOT correction. In the water target, the average reflectance in the four bands is at 2.4% for the reference while it is at 5.4% for the standard correction. The SHAOT results shows an average reflectance of 4.0%.

For cross-comparison, the image subset containing the sample spectra was also corrected using the FLAASH atmospheric correction code. The DDV method of FLAASH (named ‘K-T method’) was applied with various parameter configurations. The processing resulted in AOT values between 0.68 and 0.85 depending on the chosen method parameters, whereas the AOT from SHAOT was at 0.56. The overestimation of the AOT values in FLAASH leads to negative reflectances in the blue spectral range in water and other dark areas. The same effect was also observed if applying the ATCOR-based DDV method to the validation subset with an AOT value of 0.78. As a consequence, the average reflectance for the water targets was at 0.9% and the average for the asphalt targets was at 6.5% with very low or negative values in the blue spectral band for both surface types. It is to be noted that no in situ spectra of water are available for this data set and there may be a bias in the selected reference. In addition, the relatively low reflectance values of ADS in the blue spectral band are still to be investigated. However, these results indicate an improvement of the reflectance retrieval over dark objects from an accuracy of about 3% to below 2% in reflectance if the SHAOT based aerosol amounts are used.

5. Discussion

The sensitivity analysis has shown potential error sources and indicates that the cast shadow detection and the treatment of the adjacency effect are the most critical elements. The presented cast shadow detection algorithm performs on the same or even better level of accuracy as known cast shadow detection routines in clear sky situations. The method has been designed for operational usability and therefore mainly relies on spectral indices rather than using image statistics. However, the selection of the thresholds for shadow detection still depends on the type of sensor and aerosol scattering properties and needs to be adjusted for optimal results, depending on flight altitude and sensor system. After this initialization step, the thresholds can be kept constant, e.g., for a complete flight campaign. It is to be noted that, other than that, no parameters are to be tuned in the operational method.

The presented method of shadow detection fails on cloud shadows as it relies on spectral criteria derived from clear sky situations, which are not valid for shadows below clouds. A combination of matched filter approaches with the presented approach would need to be developed to support combined cloud shadow and cast shadow situations. Such a method should be developed for fully automatic processing in a way similar to the cloud shadow correction described by Zhu [65]. The use of at-sensor irradiance measurements may also be useful as it was shown for cloud shadows [66]. Measurements of the irradiance field could potentially also be used to increase the overall accuracy of cast shadow detection, specifically for low flight altitudes.

Comparison of the airborne AOT outputs to the AOD values of a local Aeronet station showed a reasonable relative agreement on the basis of the imaging spectroscopy data, with a tendency to overestimation. The achieved relative accuracy for the imaging spectrometer case is within the reported accuracy of Aeronet measurements of $\pm 10\%$. The spatial structure of the aerosol retrieval shows a unnatural variation on a level of approximately $\pm 10$ to $\pm 15\%$, specifically for the large scale ADS case (compare Figure 13). These uncertainties can be explained by the reference reflectance variations, the shadow mask inaccuracies, the adjacency effect model, and the irradiance model. The diffuse irradiance is difficult to approximate in shaded areas as a flat neighborhood is assumed for each pixel and the skyview factor can only be approximated. The direct irradiance from facing adjacent buildings or trees is not modeled by the atmospheric compensation approach. More sophisticated 3D modeling of the surface-to-surface radiative interactions would be required to analyze these radiative interactions for shadowed areas in more details. In addition, the radiative transfer model may be revisited to improve the results. It was observed that a higher adjacency factor q or a larger multiple scattering path radiance term (compare Equation (11)) would lead to lower AOT results that are better comparable to Aeronet based AOD measurements. Therefore, the calculation of the factor q which weights the strength of the adjacency effect may need to be revisited.

Our preliminary tests of the SHAOT method on satellite data (not presented herein) have shown various problems: low spatial resolution limits the availability of cast shadow areas, band to band mis-registration occasionally leads to inaccurate shadow signatures, and the quantization error reduces the sensitivity for the pixels in low light conditions. However, some of the resulting AOT distribution maps from satellite imagery have still been suitable for atmospheric compensation and allow for the correction of variable aerosol amounts for satellite data even if no vegetation is present in the imagery. Such samples have been found for satellite imagery containing small scale terrain shadows at a resolution of 10 m. These tests have also shown some disadvantages of the method: the method fails, if not enough shade pixels are found, if clouds or thick haze are present, and if the radiometric resolution and calibration is not good enough to measure the dynamics within shaded pixels. Further research is required to validate the usability and limits of this method for these more challenging situations.

6. Conclusions

An atmospheric correction workflow using a new aerosol optical thickness retrieval procedure has been developed. The method uses spectral signatures in cast shadow pixels in combination with an atmospheric compensation approach. The procedure is advantageous for the quantification of the aerosol scattering effect in regions where not enough dark dense vegetation areas can be found and for spatial resolutions of 5 m and below, where dark vegetation signals are highly variable. The method has been tested on high spatial resolution airborne photogrammetric data and on airborne imaging spectroscopy data. For the two investigated test cases, the retrieved aerosol amounts are reasonable. According to the sensitivity analysis and the image samples, the relative accuracy of retrieved AOT is on a level of 10–15%. Using these values for subsequent atmospheric compensation leads to improved results if compared to the use of an average aerosol amount or with the use of the dark dense vegetation aerosol retrieval outputs. Improved spectra were specifically achieved in the blue spectral range but also for the correction of the adjacency effects, visible in the near infrared spectral range. In the test cases, the level of accuracy of bottom of atmosphere reflectance was increased by this method to below 2% from a typical level of 3–4% in reflectance reported for standard atmospheric correction.

The advantages of the method is its applicability to high resolution imagery were no dark dense vegetation is available or where the high variability of vegetation signatures leads to unstable results using the DDV method. Furthermore, it is applicable to all kinds of data containing the blue, green, red and near infrared band. Therefore, it is specifically useful for digital photogrammetric data as well as for high resolution airborne imaging spectroscopy. On the other hand, the method is limited to images with spatial resolutions better than 5 m, relies on the absolute radiometric calibration, and requires high sensitivity in dark image areas. Therefore, it is not suited for 8-bit or 10-bit satellite imagery or to low cost mapping data, acquired with unmanned aerial vehicles. A further disadvantage is the required optimization of the cast shadow detection thresholds, which depends on both the sensor type and the flight altitude.

Further investigations are required to use the retrieved aerosol amounts for an accurate cast shadow removal solely based on physical models. Our results have shown a good average correction of the shadow brightness if compared to the directly illuminated areas, whereas variations within shadow areas due to variability of the terrain irradiance can not be covered by such a generic illumination correction approach. The irradiance variations are very challenging to model and a physical shadow removal should include not only an angular atmospheric irradiance model but also a pixel-wise quantification of the direct neighbor pixels irradiance impacts.

The method is completely implemented within the ATCOR-4 software package, and is ready to be used operationally. It can be engaged as an automatic atmospheric compensation procedure for relatively high resolution airborne optical scanners in clear sky conditions. By design, it is applicable to vegetated areas but also to arid environments and winter-time data acquisitions. Furthermore, the method may be used to retrieve the relative distribution of aerosol amounts, even in highly variable terrain.