Impact of Topographic Correction on PRISMA Sentinel 2 and Landsat 8 Images

Abstract

Over the past decades, remote sensing satellite sensors have significantly increased their performance and, at the same time, differed in their characteristics. Therefore, making the data repeatable over time and uniform with respect to different platforms has become one of the most challenging issues to obtain a representation of the intrinsic characteristics of the observed target. In this context, atmospheric correction has the role of cleaning the signal from unwanted contributions and moving from the sensor radiance to a quantity more closely related to the intrinsic properties of the target, such as ground reflectance. To this end, atmospheric correction procedures must consider a number of factors, closely related to the specific scene acquired and to the characteristics of the sensor. In mountainous environments, atmospheric correction must include a topographic correction level to compensate for the topographic effects that heavily affect the remote signal. In this paper, we want to estimate the impact of topographic correction on remote sensing images based on a statistical analysis, using data acquired under different illumination conditions with different sensors. We also want to show the benefits of introducing this level of correction in second level products such as PRISMA L2C reflectance, which currently do not implement it.

1. Introduction

As the performance of satellite sensors has increased, the development of accurate pre-processing procedures has become a priority to meet the need for many new and powerful remote sensing applications [1,2].

One of the most challenging issues, particularly important in soil and water surface studies, is to make the data repeatable over time and uniform with respect to different sensors and acquisition platforms, i.e., it is important to obtain data that represent the intrinsic characteristics of the observed target, rather than a peculiar state due to the specific illumination field at the time of acquisition. In fact, in addition to sensor characteristics, the signal detected by the sensor depends on the acquisition geometry and illumination conditions at the time of acquisition and is influenced by the presence of the air column and other factors, such as topography, adjacency effects, shadows, multiple reflections, clouds, BRDF and others. AC (atmospheric correction) aims to clean the signal as much as possible of unwanted contributions and to move from the sensor radiance to a quantity more closely related to the intrinsic properties of the target, such as ground reflectance or emissivity.

Over time, many techniques for AC have been considered, ranging from empirical to sophisticated physically based approaches. Many different models have been implemented, sometimes including the possibility of removing one or more unwanted factors.

In their model, Reference [3] introduced a numerical approximation technique to calculate the adjacency effect of the target surround on the remotely sensed signal for Lambertian surfaces and documented the success of the algorithm in reducing the recovered spectral reflectance error in the blue/green regions of the spectrum and reducing the total RMS spectral reflectance error, while [4,5] also introduced the effects of non-Lambertian BRDF. In their AC procedure for Landsat 8 images, Reference [6] obtained excellent results including an empirical technique for removing scattering effects due to thin cirrus in the 0.4–2.5 μm region and [7] evaluated different BRDF archetypes to represent the surface reflectance anisotropy, concluding that the MODIS BRDF product alone is able to capture the major variability in reflectance anisotropy for MODIS snow-free land cover.

Not all of these factors have the same impact on the acquired data. Depending on illumination conditions, acquisition geometry, acquired scene, sensor characteristics, and last but not least, the type of application, one unwanted contribution may be much more impactful than the others. This means that there is no one best correction method, but one method may be more appropriate for your specific needs.

In mountainous environments, for example, topographic effects heavily affect the remote signal [8,9,10]. In fact, the irradiation on a slope varies strongly with the relative angle between the azimuth of the slope and the sun, causing consistent errors in the calculation of the reflectance image if the impact of slopes is not adequately modeled before the application of correction algorithms. Very different approaches to considering topographic effects have been followed in the literature. These include empirical models [11,12,13,14], statistical or semi-physical approaches [15,16,17], and algorithms based on incident light normalization [18,19,20,21,22,23,24,25]. In empirical and statistical approaches, the results are always strongly limited to the study area taken into consideration and normally do not provide optimal results for all types of coverage [14]. Moreover, in most approaches, AC and TC (topographic correction) are treated separately, operating TC directly on the radiance image before AC or on the reflectance image after AC. In [8,26,27,28], the most common empirical and normalization approaches for TC are compared, showing that there is currently no one technique with the best ranking in all cases, as the best ranking of methods can change from one sub-scene to the next within a larger scene. In [29], the performances of semi-empirical, empirical, and statistical approaches are compared. Although the authors acknowledge that it is very difficult to evaluate the applicability of topographic correction methods using a single evaluation strategy, their results consistently suggest that the semi-empirical method performed better than the other methods. Some researchers applied a physically based model to both AC and TC [10,30,31] to retrieve the reflectance in roughness terrains. In these models, however, direct and diffuse incident lights are coupled, preventing modeling and correction of environmental impact on remote sensing data. In [32,33], some iterative procedures take into account multiple reflections to appropriately correct the data acquired in a mountainous and snowy environment, while in [34], a PLC method is adopted to model canopy reflectance in a sloping environment. Finally, some attempts have been made to use a physically based method for simultaneous AC and TC. In [35], the Integrated Radiometric Correction (IRC) method was developed, using information from satellite imagery and DEM data at 5 m resolution. The results, based on visual evaluation of the radiometrically corrected images, show significant improvements for each band and for various band compositions, while the independence between the topographically corrected reflectance and topography is demonstrated by very weak correlation coefficients compared with uncorrected data. In [36], the authors presented a combined atmospheric/topographic correction based on a database of look-up tables of atmospheric correction functions (path radiance, atmospheric transmittance, and direct and diffuse solar flux) computed with a radiative transfer code and demonstrated that an accuracy of 1–3% reflectance units can be achieved by comparing in situ reflectance measurements with surface reflectance spectra derived from airborne images. In [37], good results were obtained with an atmospheric and topographic effects correction method taking into account horizontally varying atmospheric conditions and including the height dependence of atmospheric radiance. Finally, References [38,39] show the performance of an advanced physically based correction method that allows simultaneous correction of atmospheric and topographic effects including second order corrections such as adjacency effects.

In addition to numerous research studies, the AC has been addressed directly by some data providers. In some cases, such as PRISMA and EnMAP, images are distributed with different levels of processing, including atmospherically corrected reflectance images [40,41]. In other cases, the data provider also provides tools for AC, as in the case of SENTINEL-2 with the sen2cor processor [42].

As can be seen from the literature, also AC tools used by data providers implement different approaches and different levels of correction. Such approaches are often satisfactory for general data use but may be deficient in some specific applications or environments. Of course, not all types of unwanted contributions can be corrected, and data providers must make choices based on the intended use of the data. Sentinel products, for example, include topographically corrected images, while the reflectance images included in L2C PRISMA products are obtained without considering topography [40].

In this work, using a physically based approach, we want to estimate the impact of TC on remote sensed imagery by comparing data acquired with different sensors and subjected to different correction procedures. We also want to show the benefits that can be obtained by introducing this level of correction where it is not yet applied as in the case of the PRISMA sensor.

In the next chapter, we present the approach followed to assess the impact of topography on remote sensing data. Then, we present the study area and data considered for this research. In Chapter 3, we show and discuss the main results of data analysis. In the last chapter, we present our conclusions and a brief introduction to future investigations.

2. Materials and Methods

Most physics-based approaches involve the use of an analytical invertible model that relates ground reflectance (ρ) to at sensor radiance (L) through parameters that summarize the columnar absorption and scattering properties of the entire atmosphere [43,44,45,46]. These parameters are generally calculated using software for simulating radiative transfer in the atmosphere such as MODTRAN [47,48] or 6S [49,50]. Equation (1) shows a very simple analytical model that considers the adjacency effect and is used by many AC tools, including ENVI-FLAASH [51] and ImaACor [39].

$$L=A\frac{\rho }{\left(1-S{\rho }_e\right)}+B\frac{{\rho }_e}{\left(1-S{\rho }_e\right)}+L_{path}$$

(1)

where L_path is the path radiance; S is the spherical albedo; ρ_e is the surface reflectance of a target neighborhood; A and B are parameters summarizing atmospheric properties to be calculated by a radiative transfer simulator:

$$A=\frac{{\mu }_s{E}_s}{\pi }\left(e^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_s$}\right.}+t_d\left({\theta }_s\right)\right)\left(e^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_v$}\right.}\right)$$

(2)

$$B=\frac{{\mu }_s{E}_s}{\pi }\left(e^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_s$}\right.}+t_d\left({\theta }_s\right)\right)\left(t_d\left({\theta }_v\right)\right)$$

(3)

where $E_s$ is the solar flux at the top of the atmosphere; ${\theta }_s$ and ${\theta }_v$ are the sun zenith angle and the view zenith angle, respectively; ${\mu }_{s/v}=\cos {\theta }_{s/v}$; $e^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_s$}\right.}$; $\tau$ denotes the optical thickness of the atmosphere; $e^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_s$}\right.}$ and $e^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_v$}\right.}$ are the direct downward and upward atmospheric transmittances, respectively; $t_d\left({\theta }_s\right)$ and $t_d\left({\theta }_v\right)$ are the downward and upward diffuse atmospheric transmittances, respectively. For more details on the quantities involved in Equation (1) refer to [4]. To understand how topography can be handled in a physically based correction process, it is necessary to distinguish between ground elevation compensation and slope harmonization. In fact, the impact of ground elevation and the impact of slopes, although the slopes are related to the derivative of the ground elevation, affect the at-sensor signal differently. Ground elevation affects the at-sensor signal due to the change in the height of the atmospheric column, which in turn affects the absorption and scattering properties of the atmosphere. In the absence of atmosphere, the elevation of the ground alone would not affect the remote sensing signal at all. Therefore, to compensate for ground elevation, it is not necessary to change the analytical model in Equation (1), but it will be sufficient to calculate A and B using the correct height of the atmospheric column in the radiative transfer software. For this reason, ground compensation is generally included in the AC process, even if topographic correction is not actually applied. To perform slope harmonization, however, the model in Equation (1) is not adequate and needs to be modified. In fact, it is necessary to discriminate between the direct and diffuse components of incident light to consider the reduction or intensification of the light field caused by the different solar exposures of the slopes. To this end, the first term of the second member of Equation (1) can be rewritten as follows:

$$A\frac{\rho }{\left(1-S{\rho }_e\right)}\to A_d\rho \cdot \frac{\mathit{cos}({\mathsf{\theta }}_n)}{\mathit{cos}({\mathsf{\theta }}_s)}+\left[A\frac{\rho }{\left(1-S{\rho }_e\right)}-A_d\rho \right]\cdot \frac{\left(\pi -\alpha \right)}{\pi }$$

(4)

where θ_n represents the angle between the sun and the normal to the surface and is set to π/2 for angles greater than π/2 (shadow case); α is the slope angle between the inclined surface and the horizontal plane, and A_d is another parameter that depends, like A and B in Equation (1), on atmospheric properties and which can be estimated by means of a radiative transfer simulator:

$$A_d=\frac{{\mu }_s{E}_s}{\pi }{e}^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_s$}\right.}{e}^{-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{${\mu }_v$}\right.}$$

(5)

For further details on the models related to Equations (1) and (4), see [39,45,46]. In the following, we will refer to SAC (standard atmospheric correction) to indicate atmospheric correction not including topographic correction, while TC will refer only to slope harmonization, since ground compensation is considered, if not explicitly mentioned, part of SAC. Moreover, to simplify the notation, when applying AC and TC, we will talk about ATC (atmospheric and topographic correction).

A strategy must be followed to evaluate and compare the performance/impact of different correction levels. The general idea is to perform AC by increasing the correction level and evaluate the quality of the obtained reflectance images to see if and how much improvement occurred. To evaluate the quality of reflectance images, several approaches can be followed:

• Comparison of reflectance spectra of the corrected image and ground truths.
• Comparison of reflectance spectra, within a single image, related to homogeneous covers with different solar exposure.
• Comparison of images, acquired on the same scene, under different illumination conditions.

All of the above strategies have advantages and disadvantages. Of course, comparison with ground truths is the most direct approach and is commonly followed in validation activities, both for data pre-processing and for most remote sensing applications. An application of this strategy to assess PRISMA performances can be found in [52]. However, this approach has many disadvantages, mainly related to the difficulty of comparing data acquired at very different spatial scales such as those corresponding to satellite or airborne sensors and those collected in situ. This would require adequate in situ sampling and homogeneous areas of at least 3 × 3 pixels at different slopes. The second approach has the great advantage of drastically reducing the effort of in situ activities, since their main role is to identify homogeneous areas rather than collecting spectra or samples. Nevertheless, covers can be strongly conditioned by topography, as they can be strongly influenced by the exposure to the sun and the level of humidity. This makes it difficult to identify sufficiently large homogeneous areas with the same cover at different ground elevations and differently exposed to the sun. From this point of view, the third approach is very powerful as it uses the same scene allowing to deploy the entire image with no need of in situ inspections. Nevertheless, even this approach is not entirely without its drawbacks. In addition to the problems associated with co-orthorectification of the different images, in fact, there is the need of considering different seasons to obtain different illumination on the same area, with all that this entails in terms of possible cover variations and uncertainties due to different modeling of the atmosphere. For this purpose, non-vegetated areas can be considered to minimize seasonal variations in land cover. Furthermore, having the entire image available allows to develop a reliable statistical analysis to discriminate between uncertainties due to atmospheric modeling and topographic effects.

In this work, the third approach was followed, choosing as the test site an arid mountainous area of the Argentine Andes near the border with Chile in the Los Andes department, Salta province at 24.90° S; 67.23° W (Figure 1). This area is a high plateau ranging between 3000 and 4000 m with several mountains exciding 6000 m in elevation [53]. The climate here is very rigorous; the prevailing winds from W-NW and W-SW have speeds between 2 m/s and 20 m/s, blowing for 10–18 h per day, reinforcing kinetic factors and increasing evaporation rates, even during rainy episodes [54]. Due to the extreme climatic conditions and low humidity, vegetation is completely absent with a predominantly rocky cover typical of volcanic areas [55]. Being within the tropical zone, this test site allows maximizing the differences in illumination between images, having the solar zenith moving from north to south over the course of the year.

The sensors considered are PRISMA, Sentinel-2 and Landsat-8. To obtain data with different illumination geometries, we selected one winter and one summer image for each sensor. The choice was conditioned by cloud cover and other quality parameters, and it was not easy to identify valid data for all three sensors acquired over the same area and close in time. Finally, we focused on June for the southern hemisphere winter and December for the southern hemisphere summer. Landsat-8 data were acquired on 27 June 2020 (sun zenith = 32.7; sun azimuth = 35.4; id LC08_L1TP_232077_20200627_20200823_02_T1) and 20 December 2020 (sun zenith = 62.2; sun azimuth = 93.4; id LC08_L1TP_232077_20201220_20210310_02_T1), Sentinel 2 data on 19 June 2020 (sun zenith = 35.6; sun azimuth = 30.3; id S2A_MSIL1C_20200619T142741_N0209_R053_T19JFN_20200619T175432) and 16 December 2020 (sun zenith = 63.1; sun azimuth = 92.9; id S2A_MSIL1C_20201216T142731_N0209_R053_T19JFN_20201216T175850), PRISMA data on 30 June 2020 (sun zenith = 35.6; sun azimuth = 30.3; id PRS_L1_STD_OFFL_20200630144708_20200630144712_0001) and 4 December 2020 (sun zenith = 70.1; sun azimuth = 87.1; id PRS_L1_STD_OFFL_20201204145230_20201204145234_0001).

PRISMA data were downloaded at both L1 (radiance) and L2C (reflectance) levels. The L1 data were orthorectified using the RPC (rational polynomial coefficient) files provided with the image data. For the summer image it was necessary to refine the orthorectification with ground control points identified with the help of google Earth. The L2C PRISMA images are not topographically correct. Sentinel-2 data were downloaded at L1 level (radiance), already orthorectified. The reflectance data (L2A) were instead obtained by applying the ESA send2cor processor. We have applied three different versions of the sen2cor processor, specifically, version 2.8 (plugin of SNAP 8.0) and standalone versions 2.9 and 2.10. Although online L2A products always include TC, sen2cor leaves the choice up to the user. In this way, we were able to produce two versions of the L2A product, one with SAC, the other with ATC. Finally, Landsat-8 data were downloaded at level 1 (radiance) and level 2 (not topographically corrected reflectance). Level 2 data are provided in two versions called L2C1 and L2C2. The entire dataset was also atmospherically corrected using the ImaACor tool both applying SAC and ATC. To characterize the atmospheric column, the MODTRAN standard “tropical winter” and “tropical summer” atmospheric models were used for the summer and winter images, respectively.

The radiometric quantities were simulated based on the characteristics of the sensors and the illumination conditions at the time of the acquisitions [39]. To operate ATC, the Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (SRTM) was adopted.

The PRISMA data were downloaded from the ASI PRISMA portal, while all other sensor data and the SRTM DEM were downloaded from Earthexplorer (https://earthexplorer.usgs.gov/, accessed on 14 December 2021).

For each pair of winter/summer reflectance images obtained at the various levels of correction, a normalized absolute difference (NAD) image was constructed:

$${\mathrm{NAD}}_i=\left|\frac{W_i-S_i}{\left(W_i+S_i\right)/2}\right|$$

(6)

were W_i and S_i are the i-th pixel of the winter and summer images, respectively. The wavelength dependence was omitted for simplicity. To highlight the impact of topography, the scene was divided into classes of slope, and a class spectral MRAD (Mean Relative Absolute Difference) index was defined to measure the distance between winter and summer images within a class of slope and estimate the quality of the applied correction:

$${\mathrm{MRAD}}^{Cs}=100\frac{i}{\#Cs}{\displaystyle \sum }_{iϵCs}{\mathrm{NAD}}_i$$

(7)

were Cs is the class of slope, #Cs is its cardinality, and W_i and S_i are the i-th pixel of the winter and summer images, respectively. The wavelength dependence has been omitted for simplicity.

Comparing the MRAD related to AC with the MRAD related to ATC allows to evaluate the impact of topographic correction on the correction process.

Ideally, for an inclusive correction process, MRAD should be zero over the entire spectrum, since reflectance is an inherent property of matter and should not depend on illumination conditions. In real cases, MRAD has always a finite value due to the incompleteness and approximation of the models adopted in the correction process.

For flat areas, the MRAD is mainly determined by the accuracy with which atmospheric characteristics are modeled, and therefore, it depends on the quality of the SAC procedure, but it should not depend on whether topographic correction is applied. Therefore, the MRAD corresponding to the flat areas represents a limit and a benchmark for the sloping areas, to which we tend after ATC. The degree to which we approach this value gives us a measure of the improvement that can be achieved with TC at different slopes.

Initially, we considered only two classes, which allow us to better highlight the improvement due to TC and its dependence on the wavelength. The first class (flat class) corresponds to flat areas with a slope of less than 3 degrees, the second (steep class) corresponds to very steep slopes of more than 20 degrees.

To measure the improvement achieved by TC, we define, for each class of slope, a spectral relative improvement (RI) coefficient as:

$${\mathrm{RI}}^{Cs}=100\times \frac{{\mathrm{MRAD}}_{SAC}^{Cs}-{\mathrm{MRAD}}_{ATC}^{Cs}}{{\mathrm{MRAD}}_{ATC}^{Cs}}$$

(8)

were the dependence on the bands has been omitted for simplicity.

To get a more complete picture, the improvement introduced with TC was also analyzed with reference to the degree of slope. We classified the whole scene by dividing it into one-degree classes and calculated the MRAD and the RI for each class and each sensor as we did for the FLAT and STEEP classes. Then we averaged on the bands obtaining the non-spectral MRAD and RI as a function of the slope.

Based on these indices, we performed a statistical analysis to evaluate the performance of the considered TC procedures and to evaluate the advantage of introducing this level of corrections for sensors which, like PRISMA, do not implement it.

3. Results and Discussion

As a preliminary result, before delving into the statistical analysis, in Figure 2, we show the reflectance image related to the summer and winter PRISMA images obtained with ImaACor both excluding and applying topographic correction, and the related band-averaged NAD images. The figure shows a clear flattening and decrease in correlation with topography (Figure 1) after applying topographic correction.

To better interpret the results of introducing TC for different sensors, we first present the results of a preliminary study of the impact of ground compensation on SAC procedures. For this purpose, we used ImaACor, which allows to manage the necessary parameters, and PRISMA L1, which allows a detailed spectral analysis as it provides a large number of bands. For this setup, Figure 3 shows plots of MRAD versus wavelength for flat and steep classes for a SAC procedure applied both including and excluding ground compensation.

In Figure 3, it can be seen that ground compensation has an impact mainly on the low wavelengths. In this region of the spectrum, atmospheric dispersion is greatest, and an incorrect air column eight setting causes a significant error. Something similar occurs around 2010 nm, where some CO₂ absorption bands are found. An incorrect characterization of the air column, in fact, results in an incorrect assessment of the content of atmospheric gases, which causes particularly visible errors at their absorption regions.

These results are confirmed by Figure 4, where the experiment was repeated for the Sentinel 2 sensor.

In addition to confirming the behavior of SAP at low wavelengths, Sentinel, having a broad band at 945 nm corresponding to water vapor absorption, allows us to observe for H20 what we observed for CO₂ with PRISMA.

As regards the impact of topographic correction, the first results we show concern the PRISMA sensor. Figure 5 shows MRAD vs. wavelength plots for the L2C data for flat and steep classes, along with MRAD of the ImaACor reflectances obtained from the L1C data applying both SAC and ATC. As expected, the plots show no improvement for flat surfaces after topographic correction inclusion in the correction process. The baseline error obtained for flat surfaces is more than acceptable as it is well below 10% in most of the spectrum. As explained earlier, this error represents the limit we tend to after TC and cannot be reduced even in the case of a perfect TC that can totally exclude any effect due to topography. For wavelengths above 1000 nm, the graphs relating to L2C and ImaACor standard atmospheric correction are very similar for both flat and steep classes, demonstrating the goodness of the atmospheric models adopted for the reflectance calculation. The difference between ImaACor and L2C MRAD at a lower wavelength is probably since L2C products do not consider the ground elevation of the scene, thus overestimating, as explained earlier, the height of the atmospheric column. However, there is also an unjustified and sudden increase in MRAD between 500 and 1000 nm for the flat class that needs to be investigated. When ImaACor ATC is applied, MRAD is still far from the values of the flat surfaces, highlighting the fact that residual topographic effects remain on reflectance images. This can be due to various factors, including multiple reflections between mountain slopes, incorrect sky factor modeling [38], imperfect co-orthorectification of images and DEM, shadows, and others. Despite these residuals, the inclusion of TC results in a marked improvement with an RI which, as shown in Figure 6, ranges from 40% to 50%.

As regards Landsat-8, we only show the results for L2C2 data, as the performance is slightly better than L2C1 and the results are more recent. In Figure 7, as for PRISMA, the MRAD of the L2C2 reflectances for flat and steep classes are shown, along with the MRAD related to ImaACor SAC and ATC reflectances obtained from the L1 radiances.

Even for Landsat 8, the plots show no improvement for flat surfaces after topographic correction inclusion, with a baseline MRAD of about 10%, higher than that of PRISMA. Steep slopes show slightly better results for the L2C2 reflectances than for ImaACor, when SAC is applied. Unlike PRISMA L2C, Landsat 8 standard atmospheric correction procedures consider an average ground elevation [56], and MRAD of L2C2 does not exhibit anomalous behavior at low wavelength. Again, when ATC is applied, ImaACor allows for a marked improvement over the entire spectrum, as confirmed by the RI graph of Figure 6.

Regarding sentinel 2, sen2cor allows to evaluate the impact of TC with a tool other than ImaACor. Here, we show results for processor version 2.10, which is the latest and has slightly better performance. For flat and step classes, Figure 8 shows MRAD plots as a function of wavelength of ImaACor and sen2cor SAC and ATC reflectances obtained from L1 sentinel-2 data. Sen2cor shows anomalously high MRAD values in the water absorption band and in the bands below 600 nm when SAC is applied. As for the PRISMA sensor, this is consistent with an incorrect assessment of the atmospheric column height. This problem disappears when sen2cor ATC is applied, as the ground elevation of each pixel in the image is considered and the height of the atmospheric column is no longer overestimated. Therefore, the sen2cor RI appears higher below 600 nm and in the water absorption band compared to ImaACor RI (Figure 6), as it includes both the improvement due to ground compensation and slope harmonization, while the ImaACor RI includes only the slope harmonization as ground compensation is part of the SAC.

Consistent with previous cases, for flat surfaces, the graphs show no improvements for ImaACor MRAD after TC inclusion, with a baseline error of about 10%. The improvement due to TC, albeit substantial, appears to be less than for the PRISMA and LANDSAT 8 cases, as shown by the RI graphs in Figure 6.

So far, we have shown the results for the flat areas and the steepest slopes in the scene, which are the most significant to highlight the impact of TC. To resume the results of the whole analysis, Figure 9 shows the MRAD of the different sensors and correction tools as a function of the degree of slope. Figure 9a shows the MRAD vs. slopes of L2C PRISMA reflectances and the MRAD of ImaACor reflectances obtained from L1C PRISMA radiances by applying both SAC and ATC. Figure 9b shows the MRAD vs. slope of L2C2 Landsat 8 reflectances and the MRAD of ImaACor reflectances obtained from L1 Landsat 8 radiances by applying both SAC and ATC. Figure 9c shows the MRAD vs. slope of sen2cor and ImaAcor obtained by applying SAC and ATC to L1 Sentinel-2 radiances. Finally, Figure 9d shows the trend of ImaACor RI as a function of the slope for all considered sensors.

As expected, confirming the previous results, for all sensors the MRAD related to ImaACor SAC and ATC are very similar at low slopes, in agreement with the low impact of topographic correction on flat surfaces. Then, the difference between the MRAD related to SAC and ATC procedures increases with the slope up to about 30 degrees, with the MRAD related to ATC always showing lower values.

As regards the Prisma L2C and the sen2cor reflectances, Figure 9a,c show higher MRAD values even at low slope, due to the incorrect setting of the atmospheric column height. Landsat 8 L2C2, on the other hand, taking into consideration the average elevation of the scene, shows MRAD values absolutely consistent with those of ImaACor (Figure 9b).

From Figure 9d, contrary to what one might expect, the minimum MRAD value does not correspond to the lower slope degree of the scene, resulting in negative values below 3 degrees, for some plots. This shift in the minimum is observed for any graph, for any sensor, and for both sen2cor and ImaAcor, regardless of the application of TC. This is probably because the illumination field distortion due to the slope partially compensates for some inaccuracies in the winter and summer atmospheric models. We observed this identical behavior in many other mountainous areas. However, the slope value corresponding to the minimum is always less than 3 degrees and the small discrepancy does not affect the analysis results.

All graphs show a marked increase in performance with the slope when TC is applied, up to a maximum that is reached between 20 and 30 degrees. The greatest impact of TC is observed for the PRISMA sensor, with an improvement exceeding 40% between 21 and 33 degrees. All sensors worsen when approaching the major slopes present on the scene. Among the reasons for this trend are certainly the presence of shadows and the imperfect orthorectification of the mountain ridges which make it difficult to model the topography and harmonize the illumination field at steep slopes.

4. Conclusions

We conducted a study to evaluate the impact of topography on remotely sensed data and to evaluate the benefit of introducing topographic correction in L2 products that, like PRISMA L2C, do not implement it. The study was based on the analysis of images acquired in mountainous environment by PRISMA, Landsat 8 and Sentinel 2 sensors and used two AC tools: ImaACor and sen2cor. To assess the impact of topographic correction, a statistical comparison was made between reflectance images of the same scene under different illumination conditions. Two indices were adopted, MRAD to evaluate the difference between the reflectance images as an estimator of the error associated with AC procedures, and PI to measure the improvement obtained by introducing topographic correction into the correction process. The analysis was first carried out by considering the flat areas and the steepest slopes in the scene, focusing on the spectral behavior of the MRAD and PI indices. Next, the scene was divided into 1-degree class of slope to analyze the indices as a function of the slope.

The results show an excellent agreement, for any wavelength and any degree of slope, in the case of the Landsat 8 sensor between the ImaAcor reflectances obtained by a SAC procedure and the L2C2 reflectances, which do not implement topographic correction. This demonstrates the consistency of the adopted atmospheric models. On the other hand, higher MRAD values are observed for the PRISMA L2C and Sentinel sen2cor reflectances, especially at low wavelengths and at the water absorption bands, probably due to an incorrect assessment of the height of the atmospheric column.

ImaACor RMAD values for flat areas are about 10% for all sensors. As expected, these values do not change after the introduction of TC. As we have seen, these values represent, for non-flat areas, a limit and a benchmark to aim for after topographic correction.

For all sensors, and for both ImaACor and sen2cor, a clear improvement is observed when TC is applied. In particular, ImaAcor showed an improvement of more than 40% for PRISMA L2C reflectances between 22 and 33 degrees of slope. Although the improvement is absolutely noteworthy, there is still room for improvement, as the MRAD corresponding to flat areas is still far from being achieved for all slopes. In future studies, we will investigate the possibilities of approaching this limit, by considering elements not yet introduced in ATC models, such as BRDF, sky-factor, and multiple reflections between the surfaces of the slopes.