Effect of Atmospheric Corrections on NDVI: Intercomparability of Landsat 8, Sentinel-2, and UAV Sensors

Abstract

Sentinel-2 and Landsat 8 satellites constitute an unprecedented source of freely accessible satellite imagery. To produce precise outputs from the satellite data, however, proper use of atmospheric correction methods is crucial. In this work, we tested the performance of six different atmospheric correction methods (QUAC, FLAASH, DOS, ACOLITE, 6S, and Sen2Cor), together with atmospheric correction given by providers, non-corrected image, and images acquired using an unmanned aerial vehicle while working with the normalised difference vegetation index (NDVI) as the most widely used index. We tested their performance across urban, rural, and vegetated land cover types. Our results show a substantial impact from the choice of the atmospheric correction method on the resulting NDVI. Moreover, we demonstrate that proper use of atmospheric correction methods can increase the intercomparability between data from Landsat 8 and Sentinel-2 satellite imagery.

1. Introduction

The normalised difference vegetation index (NDVI) is an outstanding index for remote sensing and is popular for its general-purpose usefulness across many different disciplines. NDVI was originally designed for qualitative evaluation of green vegetation conditions [1,2], and that is still its primary application area [3,4,5,6,7]. NDVI’s wide-ranging applications today nevertheless encompass diverse disciplines such as land cover classification and its change [8,9,10], geological hazards monitoring [11], and archaeological surveys [12]. NDVI is thus a fundamental variable in ecological applications [4,13], and it may be useful also for the correct estimation of land surface temperature [14,15]. This broad usage makes precise estimation of NDVI crucial for many applications.

Satellite images are among the most frequently used sources for NDVI data [16,17]. Landsat 8 and Sentinel-2 constitute today an unprecedented, freely available source of medium-high resolution multispectral imagery. With these satellites, we are gathering more data than ever before and with quality superior to that collected by earlier satellites [18,19]. On the other hand, spaceborne spectral imaging used for NDVI estimation is influenced by many factors. The sensor parameters (e.g., gain and offset), solar irradiance, and the sun–earth geometry are usually known and provided together with satellite images and thus may be included reliably in the corrections process. However, satellite images are also affected by absorption and scattering mechanisms of the atmosphere, which are usually unknown [20,21]. Thus, proper atmospheric correction is needed to convert the original top of the atmosphere (TOA) spectral values into an atmospherically unaffected bottom of atmosphere (BOA) spectra [16,22].

There are two basic types of approach to atmospheric correction: (i) absolute methods, in which the final product is corrected into the spectral image of surface; (ii) relative methods that result only in mutually comparable spectral images [23]. Moreover, there are two common ways of performing absolute corrections, and these are the subject of this study: (i) image-based methods, such as dark object subtraction (DOS), quick atmospheric correction (QUAC), atmospheric correction for OLI ‘lite’ (ACOLITE), or Sen2Cor, which perform atmospheric correction using only the information contained within the image itself and (ii) atmospheric modelling methods (or physical methods), such as fast line-of-sight atmospheric analysis of hypercubes (FLAASH) or the second simulation of satellite signal in the solar spectrum (6S) method, in which extra parameters for simulating the effect of the atmosphere are needed [23,24,25]. Each method is unique and provides a different result.

In ideal circumstances, a satellite spectral image after atmospheric correction could be compared to a direct ground spectral measurement to verify the result. In such a situation, evaluation of different atmospheric correction methods and their settings can lead to finding the perfect one [25,26]. For this purpose, close-range remote-sensing techniques may be applied for acquisition in near-ground spectral measurements, such as when utilising the increasingly popular unmanned aerial vehicle (UAV) technologies over small areas while benefitting from their rapid-to-deploy sensor flexibility and reasonable costs [26,27]. UAVs avoid most of the atmospheric layer, where constituents causing radiation distortions, and thus, they can be used as a source of validation data. Most users, however, have to work with images for which direct spectral measurement is not available due to temporal or spatial distance [25,26]. In the absence of direct spectral measurement, proper selection of atmospheric correction methods is crucial.

The selection of proper atmospheric correction has been widely tested for the older generation of satellites. The precursors of today’s Landsat 8, for example, were tested many times. In many investigations, simple methods, such as DOS, have been preferred because in most cases, authors have assumed that all tested atmospheric correction methods would come to fairly good results [22,23,28,29]. Nevertheless, the selection of appropriate atmospheric correction methods for current satellites, such as Landsat 8 and Sentinel-2, and their intercomparability, have not yet been fully described [24]. Our objective, therefore, is to contribute to filling this knowledge gap by testing common atmospheric correction methods on Landsat 8 and Sentinel-2 while observing their effects on NDVI. We used UAV data as a source of ground truth NDVI measurement. Finally, we tested whether atmospheric correction could improve intercompatibility between Sentinel-2 and Landsat 8 in NDVI estimation.

2. Materials and Methods

2.1. Study Site

The area of study is located on the campus of the Czech University of Life Sciences Prague, Prague, Czech Republic. The study locality, which encompasses 33.4 ha, was divided into three different land types: (i) urban, (ii) rural, and (iii) vegetated (Figure 1). The urban area consists of mostly artificial surfaces (i.e., buildings and parking lots). The rural area is made up mainly of garden cottages within surrounding gardens plus managed arboretum, the portion of green surface there exceeding 50%. The vegetated area comprises mostly grassland and bushes. The heterogeneous environment in the study locality results in the occurrence of mixed pixels within processed satellite images. This corresponds to a situation of the sort most commonly encountered, where almost no truly pure pixels of a single material are available on the study locality.

2.2. Geospatial Imagery Data

2.2.1. Satellite Imagery Collection

We used Landsat 8 (path 191 row 25) and Sentinel-2A (sensing orbit number 122) satellite images. The acquisition date of the satellites was the same as that of the UAV flight, 20 June 2017. The time difference between satellite images was only 10 min (10:00:31 for Sentinel-2A and 9:50:43 for Landsat 8, both in the GMT time zone). Both images were cloud free over study locality and were downloaded in the TOA (Level-1C for Sentinel-2 and Level 1 for Landsat 8) and also BOA (Level-2A for Sentinel-2 and Level 2 for Landsat 8). We applied the atmospheric correction methods on the whole satellite scenes, and the results were clipped on our study locality afterwards. The spatial resolution of bands used for NDVI calculations was 10 m for the Sentinel-2 product and 30 m for the Landsat 8 product. In total, 396 pixels from the Landsat satellite, 3456 from the Sentinel satellite, and more than 3 million from the UAV were used. The red- and near-infrared (NIR) band destinations of all sensors used (include UAV) are presented by the spectral response curve below (Figure 2).

2.2.2. UAV-Borne Data Acquisition and Processing

UAV-based images were collected using a Swiss-made senseFly eBee classic UAV. This is a fixed-wing vehicle with a take-off weight of 0.8 kg. The vehicle was equipped with a 4-channel AirinovMultiSPEC 4C multispectral camera that captures imagery in (i) green with a band peak of 550 nm, (ii) red 660 nm, (iii) red edge 735 nm, and (iv) near infrared 790 nm. The spectral response curves for RED and NIR are shown in Figure 2. The flight mission was performed on the same day as were the satellite flights (20 June 2017) using 80% overlaps both laterally and longitudinally and with processing by eMotion 2 ground station software. The flight altitude was 120 m above the ground. The study site was surveyed using a global navigation satellite system (GNSS) device in real-time kinematic (RTK) mode to collect ground control points for data post-processing. UAV-borne imagery was processed using Pix4DMapper 4.4.10 image matching software. The orthomosaic was built with a ground sampling distance of 10 cm. Outputs were also corrected into the values of surface reflectance using values of on-board irradiance sensor and calibration target with known albedo. A detailed description of data collection and data processing can be found in [34]. The data are available from the corresponding author, upon reasonable request.

2.3. Atmospheric Correction Algorithms

We used atmospherically uncorrected images in top of the atmosphere (TOA) format and atmospherically corrected BOA satellite images as supplied by providers. In addition, we included atmospheric correction methods provided by third parties—namely, we processed Landsat 8 images using QUAC, FLAASH, DOS, ACOLITE, and 6S. The Sentinel-2 images were processed using the same algorithms plus Sen2Cor, which cannot be used for Landsat data. Used atmospheric correction methods are described below.

2.3.1. Quick Atmospheric Correction (QUAC)

QUAC is an atmospheric correction that requires only approximate specification of sensor band locations. It uses a scene approach, and thus, it is faster than are corrections with first principle radiative models. The QUAC principle assumes that the average spectral curve of several (>10, typically 50) diverse materials from an image should have the same spectral signature as precalculated ′universal′ signature derived by averaging diverse collection of reflectance spectra from the spectral library. If there is a difference between the average library spectrum and the average from the observed endmembers spectrum, it represents an effect of the atmosphere [35]. We ran QUAC in ENVI software version 5.5 [36] on images with cloudy pixels masked out. Clouds were identified and masked out by visual inspection. All of the clouds (which covers less than 1% of images) were outside of our locality but still on the satellite image. As the input for this correction can be radiance, reflectance, or uncalibrated units, we did not perform any radiometric correction beforehand.

2.3.2. Dark Object Subtraction 1 (DOS)

DOS methods assume that non-zero signal values over supposedly zero-value dark shaded pixels are atmospheric scattering signals. This extra signal must be subtracted from the particular band [25]. Several DOS methods exist, and in the literature, these are often termed DOS1, DOS2, DOS3, and DOS4 [23]. We used the simplest DOS1 method available in the Semi-Automatic Classification Plugin version 5.3.11 within QGIS software v. 3.8 [25,37]. The method is fully automatic and requires no additional input settings.

2.3.3. Atmospheric Correction for OLI ‘lite’ (ACOLITE)

ACOLITE processor is a stand-alone application for atmospheric correction of Landsat (5, 7, and 8) and Sentinel-2 (A/B). This method is primarily used for atmospheric correction over water bodies. Nevertheless, it can be used satisfactorily also for atmospheric correction over land, as shown for example in [38]. The correction process consists of two steps: First, Rayleigh correction is performed using a lookup table generated by 6SV atmospheric correction. Second, aerosol correction is performed based upon an assumption of zero signal in SWIR bands over water bodies [39]. We ran ACOLITE version 20170718.0, in which only accurate pressure of 1008 hPa was set up according to data from the nearby Prague–Ruzyně meteorological station (9 km distant from the study area).

2.3.4. Fast Line-of-Sight Atmospheric Analysis of Hypercubes (FLAASH)

FLAASH is a first-principles atmospheric correction which uses database of precalculated results from the MODTRAN radiation transfer model [40]. MODTRAN uses built-in models and user-specified climatology data to construct vertically stratified and horizontally homogenous atmosphere models to estimate atmospheric influence. In addition to the MODTRAN model, FLAASH incorporates several other features, including correction for adjacent pixels or adjustable spectral polishing for suppression of spectral artefacts remaining after atmospheric correction [41,42]. We ran FLAASH in ENVI version 5.5. The altitude was set to 285 m above sea level in accordance with the study locality, and the atmospheric model was set up according to the manual table as a tropical aerosol model (temperature during sensing period was 29 °C at the locality). The aerosol model was rural, aerosol retrieval 2-band [K-T], and initial visibility was set to 40 km. The setting was set up using the manual as the best suitable for environmental conditions in the study locality.

2.3.5. Second Simulation of Satellite Signal in the Solar Spectrum (6S)

The 6S algorithm enables simulating signals observed by satellites between 0.25 and 4.0 µm of the spectral range. It uses a realistic model of the atmosphere and accounts for the main absorbing gases and aerosols. For a detailed description, see [43]. For this study, 6S atmospheric correction was run with the i.atcorr module in GRASS GIS 7.4.0 [44]. The atmospheric model was set up the same as in FLAASH correction as tropical. As there is no rural aerosol model as in FLAASH, we chose the best fitting equivalent-continental aerosol model. The optical depth was acquired from Aerosol Robotic Network (aeronet.gsfc.nasa.gov) as 0.25 (value for our locality in Prague was set as an average between stations Leipzig (0.35) and Vienna (0.15)). The remaining parameters were set as recommended in the manual for specific satellites.

2.3.6. Sen2Cor

Sen2Cor is a processor to perform atmospheric, terrain, and cirrus correction for Sentinel-2 data developed by the European Space Agency. The atmospheric correction part is based on a precalculated set of lookup tables [45]. A detailed description of the workflow process can be found in the user manual [46]. The processor was controlled via a SNAP standalone tool adapter. All the settings were left as default but for a few exceptions. Altitude was set for our study site (0.285 km) and, because that site is flat, we did not perform bidirectional reflectance distribution function (BRDF) correction.

2.4. Data Analysis

The normalised difference vegetation index was calculated in ENVI software [36] according to the following equation:

$NDVI=\frac{{\mathsf{\rho }}_{\mathrm{nir}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{nir}}+{\mathsf{\rho }}_{\mathrm{red}}}$

where ρ_nir and ρ_red are the near-infrared and red bands of each sensor.

As Landsat 8, Sentinel-2, and UAV data all have different resolutions, direct pixel-to-pixel comparison across all the products was not possible. To overcome this problem, we first evaluated the performance of all BOA and TOA images by comparing the median values of Red and NIR bands (Figure 3) and NDVI (Figure 4, Table 1) across the entire study locality. The pixel values were extracted in ArcGIS [47] and statistically evaluated in R [48].

In order to analyse the effect of land cover type (urban, rural, and vegetated), we used the same approach to analyse all combinations of sensors (Figure 5, Table 2), and only pixels fully inside one or another of the classification boundaries were used.

Finally, we made a statistical comparison of each NDVI value as determined by different methods and sensors. First, we calculated a semivariogram of NDVI values in order to estimate the necessary distance for sampling independent pixels. This came to 50 m for Sentinel-2 and 90 m for Landsat 8. The data points to be used for statistical purposes were sampled at a minimal distance using a regular grid sampling design in ArcGIS. The statistical comparison of NDVI data sets was then performed using correlation matrices (cor function) and pairwise Wilcoxon rank-sum tests adjusted with Bonferroni correction for multiple comparisons (pairwise.wilcox.test function) in R (Appendix A and Appendix B).

3. Results

Six atmospheric correction methods were computed and, together with the producer’s atmospheric correction (Level 2/2A, hereinafter Level2), non-calibrated images (TOA), and UAV-borne imagery, were compared with one another.

To determine the effect of atmospheric correction on raw data, we investigated changes in satellite RED and NIR bands used for NDVI calculation. Each correcting method, with the exception of QUAC, preserved a histogram of a shape similar to that seen for the TOA (Figure 3). The values corrected by the QUAC method were significantly shifted in comparison to the others, and especially in the RED band of the Landsat 8 image and NIR band of the Sentinel-2 image.

Median values, together with an interquartile range of NDVI values for each tested method, are presented in Table 1. Each atmospheric correction method increased the NDVI values, in contrast to those from TOA images. This increase is greater in the case of Landsat 8 data than Sentinel-2. Figure 4 reveals that the QUAC method resulted in higher NDVI values, in contrast to those from other atmospheric correction methods. Homogeneity was generally higher for the Sentinel-2 values, in contrast to the Landsat 8 values (Figure 4). The differences between median NDVI values for DOS1, FLAASH, ACOLITE, 6S, and Level2 were less than 0.020 for Landsat 8 image and 0.039 in the case of the Sentinel-2 image.

The results show the level to which NDVI values are sensor dependent. Median NDVI estimated by TOA Landsat 8 values was lower by 0.133, in contrast to Sentinel-2. After atmospheric correction, the results were reversed in the sense that the Landsat 8 NDVI values were higher than were the NDVIs estimated from the Sentinel-2 image. All atmospheric correction methods resulted in smaller differences between NDVI from Landsat 8 and Sentinel-2 images. The smallest difference was achieved by FLAASH atmospheric correction. UAV values were higher than all other tested values. The values closest to the UAV median were achieved by Landsat 8 QUAC atmospheric correction, which leads to distinctly higher NDVI, compared to those from other atmospheric correction methods.

We also calculated the correlation between the resulting NDVI images from each atmospheric correction method. All atmospheric correction methods had a high correlation with one another, as well as with the TOA product, as Pearson’s r value was never less than 0.98 in the case of Landsat 8 data, and not less than 0.99 in the case of Sentinel-2. A correlation matrix is presented in Appendix A.

The similarity of results from different atmospheric correction methods was also tested statistically. The results show that median NDVI from most of the atmospheric correction methods was substantially different for the study locality, but some atmospheric methods were statistically indistinguishable. In the case of NDVI from Landsat 8 image, the Level2 product was not statistically different from those from FLAASH, ACOLITE, and 6S, and that for 6S did not differ from that for DOS1. In the case of NDVI values from the Sentinel-2 image, that for Level2 was not different from those for ACOLITE and DOS1, and the value for DOS1 did not differ from that for Sen2Cor. Complete results may be found in Appendix B.

The effect of landscape type on median NDVI values is described next (Figure 5). NDVI values in corrected BOA bands changed across landscape types. In general, however, those values followed the trends described above. TOA values resulted in lower NDVI values. The QUAC method and UAV-borne data resulted in high NDVI values, in contrast to other atmospheric correction methods. Heterogeneity across the methods was more demonstrable in the case of Landsat 8 values. In both Landsat 8 and Sentinel-2 data, greater dispersion across the NDVI values was obvious for the urban landscape type, while across the rural landscape, the NDVI values were quite homogenous, and especially so for Landsat 8.

The intercomparability of results is crucial also between different satellite types. Table 2 summarises the differences between Sentinel-2 and Landsat 8 using different atmospheric correction methods. Overall, small differences were achieved when using FLAASH and 6S atmospheric corrections. On the contrary, QUAC, DOS1, the Level2 product, and products without atmospheric correction (TOA) resulted in larger NDVI differences between these two satellite platforms.

4. Discussion

We tested six atmospheric correction methods and evaluated their effect on the resulting NDVI as observed by Sentinel-2 and Landsat 8. In comparing TOA and BOA products, we found that atmospheric correction has a substantial impact on the resulting NDVI and also its intercomparability between satellite sensors.

First, we investigated the impact of atmospheric correction methods on values of RED and NIR bands (Figure 3) used for NDVI calculation. While the results of methods correcting Landsat 8 were homogenous, substantial variations were found across Sentinel-2 data, especially for the NIR band. The remaining values seemed about equivalent. Especially the 6S, DOS1, FLAASH, and ACOLITE, as well as the producer-provided (Level2), corrections gave very similar results for both Sentinel and Landsat data. The exception was the QUAC method, which results in low values of Landsat 8 RED band in contrast to other atmospheric correction methods and high NIR values in the case of Sentinel-2.

The resulting NIR and RED bands further impact the difference between NDVI from BOA and TOA images. That was especially the case of Landsat 8 imagery, where the average NDVI for TOA was 0.306, while the average NDVI from corrected images was 0.553. This 81% increase for Landsat 8 is in contrast to a report of Hadjimitsis et al. [29], who found the average difference between NDVI from corrected BOA images and TOA Landsat 8 images to be only 18% in the area of southern Cyprus. In our case, we found an increase of just 18% (0.497 versus 0.421) in the case of Sentinel-2 images. The increase of reflectivity in the near-infrared spectrum and decreasing of reflectivity in the visible spectrum were also observed by Valdivieso [38] after several atmospheric corrections of Sentinel-2 data in Spain. This is in line with Xie [49], who observed an increase of NDVI after the 6S atmospheric correction of the Landsat-7 image over forest, grassland, and desert areas in China.

When using the same sensor but different atmospheric correction methods (DOS1, FLAASH, ACOLITE, 6S, SEN2COR, but excluding QUAC), the similarity found between the NDVI values indicates that these are comparable in most cases, regardless of the correction method applied. This can be seen in Appendix A, where a strong correlation is shown between each pair of tested atmospheric corrections. Despite this strong correlation, most of the atmospheric correction methods still lead to statistically different results (Appendix B). The maximum difference between medians of resulting NDVI from different atmospheric corrections of Landsat 8 images was 0.019, and it was 0.039 in the case of Sentinel-2 images. Even though these differences may appear small, they may nevertheless be crucial for many applications. Gouveia et al. [50] studied drought assessment; for example, they estimated stressed vegetation to be represented by NDVI anomalies on the magnitude of just −0.025 NDVI. Although Nazeer et al. [22] revealed that atmospheric correction methods may perform variously across different land cover types, our results indicate that the effect of atmospheric correction remains similar across different land cover types when results from a single type of sensor are considered. We found the impact of different land cover types on NDVI to appear when different sensors were involved, as is discussed below.

Another aspect affected by atmospheric correction is the distinctiveness of the resulting raster of NDVI. The broader histogram of NDVI allows us to distinguish more detailed differences in the image (Figure 4). In the case of Landsat 8, FLAASH correction yields the largest range of values, and thus, NDVI results in the largest IQR (Table 1). In the case of the Sentinel-2 Level 2 product, Sen2Cor and also FLAASH resulted in NDVI with the broadest histograms after atmospheric correction. On the contrary, the QUAC method and product without correction (TOA) resulted in a narrow variety of NDVI values. Such a result could negatively affect the sensitivity of NDVI and make it more difficult to distinguish different vegetation types and/or conditions.

Our results also demonstrate that sensor intercomparability can be enhanced by applying the proper atmospheric correction method. In our case, the overall difference in median NDVI between TOA images of Landsat 8 and Sentinel-2 was 0.133. The average difference between medians after atmospheric corrections decreases to 0.043 NDVI. The best result, with a difference between median NDVI values of only 0.013, was achieved by the FLAASH method. We did not assume the smaller difference after atmospheric correction to be due to different atmospheric conditions in each sensing period, because the time span between the Landsat 8 and Sentinel-2 satellite observations was only 10 min in stable, clear sky conditions. The explanation could be that even with similar RED and NIR band destinations the spectral response function is unique for each sensor (see Figure 2). Hence, each sensor is influenced by the atmosphere differently as the atmospheric effect is constantly changing due to the manifestation of different atmospheric gasses and particles in different wavelengths. This unique sensitivity of the sensor to the atmospheric effect can be minimised by atmospheric correction methods.

We generally observed smaller differences between NDVI values estimated from different sensors in the rural environment, followed by vegetated areas (see Table 2). The biggest differences between sensors were in urban areas. The smallest difference was achieved by FLAASH, followed by 6S. This is in line with results from Ke et al. [16], who found that FLAASH and 6S were in better agreement with MODIS spectrometer compared to DOS. Surprisingly, in our study, the vegetated areas resulted in larger differences with FLAASH than did rural and even urban areas, which is contradictory to the findings of Ke et al. [16].

We also tested the capability of the UAV sensor as a source of ground truth data. The differences in band specification (Figure 2) show that the UAV-mounted sensor has a lower spectral band position of the NIR band, which resulted in a higher NDVI across all landscape types (Figure 5). In general, we found a lower variation of NDVI in UAV sensors, compared to satellite imagery. Higher standard deviation values in satellite imagery than those of the UAV in the case of the SAVI index were also confirmed by Messina et al. [51]. The overall difference in median NDVI between the average atmospherically corrected and UAV product was 0.184 in the case of Landsat 8, and 0.229 in the case of Sentinel-2 NDVI. This is a meaningful difference for many studies, and authors tend to use UAVs as ground truth data, which is especially the case of precision agriculture [52,53,54]. Considering the different range of NIR band of tested platforms, we found the shift to the shorter wavelengths of the UAV band significant, and thus, we conclude, contrary to Ryu et al. [55], the differences in measured values observed at our locality limit the possibility for using UAV as a collector of ground truth NDVI values for satellites. This limit could be overcome if the spectral band adjustment factors for UAV cameras were known [56] at least for the most common satellites. Unfortunately, these factors are not usually provided by vendors.

Further validation of our results under different atmospheric conditions would be beneficial for a deeper understanding of atmospheric correction method performance. Unfortunately, it is unique to have both satellites (Landsat 8 and Sentinel-2) scanning one area at one moment. It happened only four times, during the vegetation season (May–September) 2017, over our study area and only once with clear sky conditions.

5. Conclusions

The proper use of atmospheric correction methods is crucial and has a significant impact on NDVI estimation. The results from most of the tested atmospheric corrections (DOS1, FLAASH, ACOLITE, 6S, SEN2CORE) are with some limitations, comparable with one another when just one sensor is used. It also has been proven that atmospheric correction methods lead to NDVI results being more comparable between Landsat 8 and Sentinel-2. The smallest difference of 0.013 median NDVI between sensors was achieved by the FLAASH atmospheric correction method. FLAASH also yielded results well comparable to those from other atmospheric correction methods, and it produced the NDVI raster with the best distinctiveness. Our results suggest that the intercomparability between satellite sensors can be enhanced by proper choice and application of atmospheric correction methods.

Author Contributions

D.M. with I.M. conceived the ideas. D.M. designed the methodology, processed satellite data, and analysed all data. J.K. was responsible for UAV data acquisition and processing. D.M. and J.K. led the writing of the manuscript together with the significant supervision of S.L.-C.M. and I.M. All authors contributed critically to the drafts and gave final approval for publication. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Technology Agency of the Czech Republic (Environment for Life program), project: Possibilities for updating map layers of NATURA 2000 biotopes using advanced remote sensing methods, Grant Number: SS01010046.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A

Table A1

Pearson’s r correlation between NDVI values acquired by all atmospheric correction methods, TOA, and UAV.

Landsat 8	UAV	TOA	DOS1	6S	ACOLITE	FLAASH	QUAC	Level2
UAV		0.8878	0.8821	0.8901	0.8906	0.8909	0.8904	0.8882
TOA	0.8878		0.9807	0.9887	0.9903	0.9915	0.9889	0.9875
DOS1	0.8821	0.9807		0.9963	0.9959	0.9955	0.9966	0.9959
6S	0.8901	0.9887	0.9963		0.9999	0.9997	0.9998	0.9992
ACOLITE	0.8906	0.9903	0.9959	0.9999		0.9999	0.9998	0.9993
FLAASH	0.8909	0.9915	0.9955	0.9997	0.9999		0.9996	0.9993
QUAC	0.8904	0.9889	0.9966	0.9998	0.9998	0.9996		0.9989
Level2	0.8882	0.9875	0.9959	0.9992	0.9993	0.9993	0.9989
Sentinel-2	UAV	TOA	DOS1	6S	ACOLITE	FLAASH	QUAC	Level2	Sen2Cor
UAV		0.8994	0.8991	0.8988	0.8993	0.8996	0.8988	0.8990	0.8993
TOA	0.8994		0.9985	0.9982	0.9992	0.9995	0.9978	0.9982	0.9988
DOS1	0.8991	0.9985		0.9995	0.9995	0.9996	0.9996	0.9991	1.0000
6S	0.8988	0.9982	0.9995		0.9998	0.9996	0.9994	0.9996	0.9995
ACOLITE	0.8993	0.9992	0.9995	0.9998		0.9999	0.9993	0.9995	0.9996
FLAASH	0.8996	0.9995	0.9996	0.9996	0.9999		0.9991	0.9995	0.9997
QUAC	0.8988	0.9978	0.9996	0.9994	0.9993	0.9991		0.9990	0.9996
Level2	0.8990	0.9982	0.9991	0.9996	0.9995	0.9995	0.9990		0.9992
Sen2Cor	0.8993	0.9988	1.0000	0.9995	0.9996	0.9997	0.9996	0.9992

Appendix B

Table A2

p-values of pairwise Wilcoxon rank-sum tests with Bonferroni adjustment for multiple comparisons of NDVI values resulting from the tested atmospheric correction methods applied to Landsat 8 and Sentinel-2 imagery.

Landsat 8	TOA	Level2	QUAC	FLAASH	ACOLITE	6S
Level2	<0.001
QUAC	<0.001	<0.001
FLAASH	<0.001	0.098	<0.001
ACOLITE	<0.001	1.000	<0.001	<0.001
6S	<0.001	0.069	<0.001	<0.001	<0.001
DOS1	<0.001	0.002	<0.001	<0.001	<0.001	1.000
Sentinel-2	TOA	Level2	QUAC	FLAASH	ACOLITE	6S	DOS1
Level2	<0.001
QUAC	<0.001	<0.001
FLAASH	<0.001	<0.001	<0.001
ACOLITE	<0.001	0.610	<0.001	<0.001
6S	<0.001	<0.001	<0.001	<0.001	<0.001
DOS1	<0.001	1.000	<0.001	<0.001	<0.001	<0.001
Sen2Cor	<0.001	<0.001	<0.001	<0.001	0.010	<0.001	1.000

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Location and detail of the study area. Three landscape types are distinguished by colour. The map corresponds to WGS 1984 Web Mercator Auxiliary Sphere projection.

Figure 2. Spectral response functions of RED and NIR bands of sensors used and atmosphere transmittance. Sources: ESA [30], NASA [31], User Manual multiSPEC 4C camera [32], and Modtran [33].

Figure 3. Histogram with highlighted median values of non-corrected TOA and corrected BOA (Level-2, QUAC, FLAASH, ACOLITE, 6S, DOS1, Sen2Cor) bands for NDVI over the study area.

Figure 4. Comparison of different atmospheric corrections and visualisation of resulting histogram changes for each method. The vertical line indicates the median.

Figure 5. Comparison of NDVI values for three landscape types calculated after using various methods for correcting atmospheric effect on surface NDVI and when using Landsat 8 and Sentinel-2 satellite images. Horizontal lines representing median, 25th, and 75th percentiles, dots are outliers (further than 1.5 times the IQR from 1Q and 3Q, respectively), and whiskers symbolise minimum and maximum without outliers.

Table 1

Resulting NDVI median, interquartile range (IRQ), and differences between Landsat 8 and Sentinel-2 across study area for each method applied, together with UAV values measured on the same area.

	Landsat 8		Sentinel-2		Difference
	Median	IQR	Median	IQR	Median	IQR
TOA	0.334	0.148	0.468	0.279	−0.133	−0.130
Level2	0.584	0.216	0.534	0.363	0.049	−0.147
QUAC	0.717	0.176	0.635	0.277	0.082	−0.101
FLAASH	0.564	0.270	0.551	0.329	0.013	−0.060
ACOLITE	0.573	0.236	0.530	0.296	0.043	−0.060
6S	0.583	0.246	0.563	0.303	0.019	−0.057
DOS1	0.579	0.224	0.525	0.302	0.055	−0.078
UAV	0.784	0.193	0.782	0.229
Sen2Cor			0.529	0.355

Table 2

Difference (Landsat 8–Sentinel-2) of median NDVI over study locality and summary of all differences across different landscape types.

	DOS1	6S	ACOLITE	FLAASH	QUAC	Level2	TOA	Ʃ
Rural	0.045	0.017	0.043	0.014	0.067	0.013	−0.167	0.031
Urban	0.097	0.058	0.074	0.046	0.150	0.146	−0.037	0.535
Vegetated	0.073	0.046	0.069	0.049	0.087	0.040	−0.142	0.223