Topographic Correction of LAPAN-A3/LAPAN-IPB Multispectral Image: A Comparison of Five Different Algorithms

The study was conducted in South Sulawesi Province, Indonesia, which is located in the middle of Indonesian archipelago in an ancient and currently inactive volcanic region. The volcanic rock formation is geologically known as Baturappe-Cindakko (Tpbv), as shown in Figure 1. The slope area comprised mostly of agricultural land (53%) and forest (41%). Small portions of the study area are also covered with bushes and shrubs (4%), settlements (1.5%), as well as water-body (0.5%) as shown in Figure 1C. The relief is undulating, with slopes varying from 0 to more than 60 degrees (Fig. 2A, 2B, 2D, 2E). A 23×23 km rectangular subset was chosen as the area of interest (AOI). The AOI was specifically selected at mountainous region with different facing slopes to better understand the effect of topographic correction (Fig. 2C and 2F).

Fig. 1

Fig. 2

The LAPAN-A3/LAPAN-IPB (LA3) used in this study were acquired on August 28th 2019 at 09.34 AM local time. The solar azimuth was at 69.4437°, with an elevation of 50.75°. LA3 has four spectral bands, ranging from blue to near infrared. The spatial resolution is 15 meters (Zylshal et al. 2018), and the LA3 imagery came with WGS84 projection. After geometric correction, the LA3 data was then transformed into UTM projection (Zone 50 South).

For the chosen correction methods, terrain elevation information is needed in order to simulate the acquisition lighting conditions. Various digital elevation models (DEMs) are available free of charge and cover almost the entire globe; namely, the Shuttle Radar Topographic Mission (SRTM); the ASTER Global Digital Elevation Map (GDEM); and the ALOS Global Digital Surface Model (AW3D30). The latter is currently the newest freely available dataset, being released in March 2017. It utilizes the Advanced Land Observing Satellite (ALOS), specifically the PRISM (Panchromatic Remote-sensing Instrument for Stereo Mapping), to compute the elevation, which delivers with approximately 30 meters of spatial resolution (Japan Aerospace Exploration Agency 1997, Takaku et al. 2016, JAXA 2017, Takaku and Tadono 2017). The DEM was used to compute the slope and aspect of the terrain; a previous study found that AW3D30 performed the best (Zylshal 2019).

Most previous studies were performed on orthorectified images (Riaño et al. 2003, Gao and Zhang 2009a, Richter et al. 2009, Hantson and Chuvieco 2011, Balthazar et al. 2012, Nichol and Hang 2013, Vanonckelen et al. 2013, 2014, Gao et al. 2016, Phiri et al. 2018). Due to the lack of Rational Polynomial Coefficient (RPC) sensor model parameters in LA3 metadata, for now, the geometric correction was performed using image-to-image correction (Jensen 2005).

A corresponding Landsat-8 OLI Tier 1 data were used as the reference image. Fifty ground control points (GCPs) were selected, and the root mean square error (RMSE) was kept at less than 1 LA3 pixel size. Before the topographic correction was performed, it was necessary to clip both the DEM and LA3 to the same extent. Therefore, the AOI described in Figure 1 was used as the clipping boundary. Slope and aspect were then derived from the DEM, as shown in Figures 2B and 2C, respectively. Figures 2D, 2E, and 2F show the distribution of terrain elevation, slope and aspect in the study area. In this preliminary study, the topographic correction was performed on LA3 digital number without prior radiometric or atmospheric correction. Each resampling procedure used the nearest neighbour algorithm, due to its ability to preserve the original pixel value (Parker et al. 1983, Zylshal et al. 2017). As explained on previous section, the LA3 DN value was used in further analysis.

This initial study focuses on topographic correction using illumination modelling (Pimple et al. 2017). This approach can be categorized into two groups: that which assumes Lambertian conditions, in which the reflectance is independent of the observation angle; and that which considers bidirectional reflectance (non-Lambertian) (Hantson and Chuvieco 2011).

Figure 3 illustrates the illumination modelling (IL) used for the topographic correction. The first step is to calculate the illumination modelling (Riaño et al. 2003) using equation 1 below: (1) IL=cosγi=cosθscosθn+sinθssinθncos{ϕn−ϕs} {\rm{IL}} = \cos {\gamma _i} = \cos {\theta _{\rm{s}}}\cos {\theta _{\rm{n}}} + \sin {\theta _{\rm{s}}}\sin {\theta _{\rm{n}}}\cos \{ {\phi _{\rm n}} - {\phi _{\rm{s}}}\} where: IL

– the local solar illumination angle,

Fig. 3

DEM was used to obtain θ_s and ϕ_s. The simplest method is cosine correction, as proposed by Teillet et al. (1982). This algorithm assumes the ground as a Lambertian surface. It is easier to apply since it does not need any external parameters, and is defined as: (2) ρH=ρT(cosθnIL) {\rho _H} = {\rho _T}\left( {{{ {\it cos}\, {\theta _n}} \over {IL}}} \right) where: ρ_H

– the reflectance of a horizontal surface,

While this method is simpler to apply, several studies have found that for the low illuminated areas where cos γ_i is close to zero, it tends to produced overcorrection (Riaño et al. 2003, Gao and Zhang 2009b). Civco (1989) therefore proposed an alternative version of cosine correction by adding average illumination conditions (equation 3). IL¯ \overline {IL} Denotes the average IL value. Both of these algorithms are independent of the wavelength. One of the most cited non-Lambertian methods was proposed by Minnaert (1941), which was originally used to assess the roughness of the moon's surface, but then more widely adopted in EO data (equation 4).

This approach introduced a k factor, which is commonly referred to as the Minnaert constant, whose value ranges from 0 to 1. A k value of 1 indicates a perfect Lambertian surface. The k factor is calculated for each band by linearization of equation (4), so equation (4) becomes equation (5) (Riaño et al. 2003) (5) ln(ρT)=ln(ρH)+Kln(ILcosθn) \ln \left( {{\rho _T}} \right) = \ln \left( {{\rho _H}} \right) + K\ln \left( {{{IL} \over {{\it cos} \;{\theta _n}}}} \right)

K and ln(ρ_H) are the linear regression coefficients, and ρ_H is constant for the whole image (Riaño et al. 2003). Several studies has shown that the Minnaert algorithm has been able to improve accuracy compared to cosine correction (Justice et al. 1981, Itten and Meyer 1993, Richter et al. 2009). Other studies then modified equation (4) to include slope information (Colby 1991, Riaño et al. 2003), as defined below: (6) ρH=ρTcosθs(cosθnILcosθs)K {\rho _H} = {\rho _T} {\it cos}\, {\theta _s}{\left( {{{{\it cos}\, {\theta _n}} \over {IL\, {\it cos}\, {\theta _s}}}} \right)^K}

Civco (1989) proposed a two stage normalization method. First, a hillshade model is simulated to match illumination conditions based on the Sun's azimuth and elevation at the satellite acquisition time. Each of the original bands is then transformed into a topographically normalized image using equation (7) below: (7) ρH=ρT+((ρTIL¯−ILIL¯)Cλ) {\rho _H} = {\rho _T} + \left( {\left( {{\rho _T}{{\overline {IL} - IL} \over {\overline {IL} }}} \right){C_\lambda }} \right) where IL¯ \overline {IL} is the mean value for the south-facing and north-facing slopes, and C_λ is the correction coefficient for each band λ, calculated as equation (8): (8) Cλ=(Sλ−Nλ)[(Nλ×μN−μkμk)−(Sλ×μs−μkμk)] {C_\lambda } = {{\left( {{S_\lambda } - {N_\lambda }} \right)} \over {\left[ {\left( {{N_\lambda } \times {{{\mu _N} - {\mu _k}} \over {{\mu _k}}}} \right) - \left( {{S_\lambda } \times {{{\mu _s} - {\mu _k}} \over {{\mu _k}}}} \right)} \right]}} where: N_λ

– the mean of the away-facing slope from the sun for the uncorrected image,

The performance of each algorithm was evaluated using two different assessments, namely qualitative inspection, based on visual appearance, and quantitative evaluation. By comparing the same image before and after correction, visual assessment was used to observe whether there was over/under correction. Overcorrection would result in brighter pixels for the sun-shaded compared to the sun-facing slopes (Gao and Zhang 2009b). Under correction, on the other hand, would result in darker sun-shaded slopes compared to the sun-facing ones.

Under- or overcorrection was also evaluated based on the regression fitting line between the local solar illumination angle (cos γ_i) and the pixel value, before and after correction (Vincini and Frazzi 2003, Gao et al. 2013, Zhang and Gao 2011). A successful correction method should be able to reduce the correlation (r) between cos γ_i and the pixel value. For this purpose, 3942 samples were randomly generated. The pixel value before and after correction, as well as cos γ_i for each sample per band, were then taken and fitted into a linear regression. A t-test was then performed to test the significance at 95% confidence level.

For quantitative analysis, changes in the coefficient variation (CV) of the pixel values before and after correction were used. This method has been widely employed to validate topographic correction (Gao and Zhang 2009b, Hantson and Chuvieco 2011). Changes in CV are also referred to as dispersion indices (Gao and Zhang 2009a, Zhang and Gao 2011).

The sample points previously created were then used to calculate the mean and standard deviation value from both before and after image correction. The coefficient variation was then calculated using equation (9) below: (9) CV(%)=δμ×100 CV\left( \% \right) = {\delta \over \mu } \times 100 where δ is the standard deviation, and μ is the mean pixel values from the samples.

The changes in the CV were then calculated using equation (10) (Richards 1995, Vanonckelen et al. 2014, Pimple et al. 2017, Zylshal 2019): (10) CVdifference=CVbefore−CVafter C{V_{difference}} = C{V_{before}} - C{V_{after}}

The CV_difference was calculated for each band and each algorithm and the results compared. Successful topographic correction should be able to reduce the variability within each band; that is, giving a positive value of CV_difference.

Figure 4 shows the images resulting from each of the five topographic correction algorithms applied to the AOI. The original image is shown in Figure 4A using a red-NIR-blue composite to better enhance the areas of dense vegetation in the sloped region. Visually, cosine correction (Fig. 4B) performs worst, with visible overcorrection in the hilly region. Improved cosine correction (Fig. 4C) and two-stage normalization (Fig. 4F) are both able to reduce such overcorrection. However, they exhibit an inverse effect, with some of the valleys appearing as hills. Minnaert correction (Fig. 4D), as well as the modified Minnaert (Fig. 4E) performed the best. Both these methods were able to flatten the hilly region and reduce the topographic effect.

Fig. 4

Figure 5 shows scatter plots and regression lines between the pixel value and local illumination (cos γ_i) over the near-infrared band for the 3,942 samples taken. Figure 5A shows a clear and statistically significant correlation between the pixel value of the NIR band and local illumination (R² = 0.2837). Figure 5B confirms the overcorrection shown in Figure 4B, with an increase in the correlation, as well as the regression line being tilted in another direction. Increased correlation is also found in both Figure 5C and 5D (improved cosine correction and two-stage normalization respectively, albeit at not such high levels as with cosine correction. Only the Minnaert correction (Fig. 5D) and modified Minnaert (Fig. 5E) were able to reduce the correlation coefficient, with the latter performing the best, with the highest correlation coefficient reduction (Table 1). Additional quantitative assessment by calculating the coefficient variation before and after correction is shown in Table 1.

Fig. 5

For the red and green bands, Minnaert correction performed the best, with the highest correlation coefficient reduction (2.52 and 0.53 respectively). For the blue band, all five correction methods were unable to reduce the correlation coefficient. Minnaert correction performed the best, with the lowest increase in CV_difference (−0.16). Observing the number of bands successfully corrected, Minnaert correction and modified Minnaert correction were able to correct three out of the four LA3 bands. The improved cosine correction and the two-stage normalization methods were only able to correct one of the four LA3 bands, while cosine correction performed the worst, with no bands successfully corrected.