A New Single-Band Pixel-by-Pixel Atmospheric Correction Method to Improve the Accuracy in Remote Sensing Estimates of LST. Application to Landsat 7-ETM+

Abstract

Monitoring Land Surface Temperature (LST) from satellite remote sensing requires an accurate correction of the atmospheric effects. Although thermal remote sensing techniques have advanced significantly over the past few decades, to date, single-band pixel-by-pixel atmospheric correction of full thermal images is unsolved. In this work, we introduce a new Single-Band Atmospheric Correction (SBAC) tool that provides pixel-by-pixel atmospheric correction parameters regardless of the pixel size. The SBAC tool uses National Centers of Environmental Prediction (NCEP) profiles as inputs and, as a novelty, it also accounts for pixel elevation through a Digital Elevation Model (DEM). Application of SBAC to 19 Landsat 7-ETM+ scenes shows the potential of the proposed pixel-by-pixel atmospheric correction to capture terrain orography or atmospheric variability within the scene. LST estimation yields negligible bias and an RMSE of ±1.6 K for the full dataset. The Landsat Atmospheric Correction Tool (ACT) is also considered for comparison. SBAC-ACT LST deviations are analyzed in terms of distance to the image center, surface elevation, and spatial distribution of the atmospheric water content. Differences within 3 K are observed. These results give us the first insight of the potential of SBAC for the operational pixel-by-pixel atmospheric correction of full thermal images. The SBAC tool is expected to help users of satellite single-channel thermal sensors to improve their LST estimates due to its simplicity and robustness.

1. Introduction

Land Surface Temperature (LST) is a key magnitude in the characterization of the surface-atmosphere energy exchanges, evapotranspiration, meteorology, climatology and hydrology [1,2,3,4,5,6]. Remote sensing in the thermal infrared (TIR) allows the estimation of spatially distributed LST. Satellite-derived LST is used in a wide variety of scientific explorations and societal applications, such as water resource management, detection of droughts and heat waves, volcano activity, urban heat islands, etc. [7,8,9,10].

Top-of-atmosphere measurements must be corrected from atmospheric and emissivity effects to convert brightness temperatures into LSTs. TIR-based science has advanced significantly over recent decades, and different atmospheric correction methods have been proposed depending on the features of the sensor. Split-window methods apply to those sensors provided by at least two spectral channels, usually at 11 and 12 × µm. This has been the correction method traditionally applied to NOAA/Advanced Very High Resolution Radiometer (NOAA/AVHRR, [11,12,13,14]), EOS Terra-Aqua/Moderate Resolution Image Spectroradiometer (EOS/MODIS, [15,16,17]), Envisat/Advanced Along Track Scanning Radiometer (Envisat/AATSR, [17,18]) or the Meteosat Second Generation/Spinning Enhanced Visible and Infrared Imager (MSG/SEVIRI, [19,20,21]) or Suomi National Polar-Orbiting Partnership/Visible Infrared Imaging Radiometer Suite (s-NPP/VIIRS, [22,23]). When the sensor is provided with three or more bands in the TIR domain, Temperature and Emissivity Separation (TES) method [24] can be applied. Once corrected from atmospheric effects, at-ground TIR radiances are used as input in the TES approach. TES method was originally developed to be applied to the EOS-Terra/ASTER [24,25]. LST products of MODIS (MOD21) and s-NPP/VIIRS are also based on TES technique [26,27].

Sensors provided with a single channel in the TIR domain need a simulation of the atmospheric effects from an estimation of the atmospheric water vapor and air temperature profiles. In such cases, Single Channel (SC) techniques are required. Different sources can be used for the atmospheric profiles needed. Launching a local atmospheric sounding close in time to the sensor overpass is the most accurate solution, although rarely possible. Currently, atmospheric profiles can be obtained from satellite sounders such as the EOS-Aqua/Air Infrared Sounder (EOS-Aqua/AIRS). MODIS product MOD07 is also available. However, not always the sensor of interest is concurrent with the sounder overpass and the results are not desirable [28]. Finally, some meteorological models forecast atmospheric profiles with a fixed periodicity and very low spatial resolution (3–6 h and 0.25°–2° latitude/longitude). Li et al. [29] compared the results obtained between the National Center for Environmental Prediction (NCEP) and MOD07 for Huan-Jing 1B /Infrared Scanner (HJ-1B/IRS) data in two sites in China, covering water, bare soil, wheat and corn, and their results yielded differences around ±1.2 K for both sources. More recently, Meng and Cheng [30] conducted an evaluation on several reanalysis products (MERRA, NCEP, ERA and JRA) using 32 L8-TIRS images over Gansu province, China. Almost no difference in terms of LST is observed by these authors with RMSE values within ±2 K. Similar results were obtained by Duan et al. [31], over sands, grasslands and snow, using NCEP data in conjunction with MODTRAN 5 code to correct L8-TIRS data from atmospheric effects.

The Atmospheric Correction Tool (ACT) developed by Barsi et al. [32,33] directly provides the atmospheric correction parameters at a given geographical point for band 6 of Landsat 5-TM or L7-ETM+, and band 10 of L8-TIRS. The ACT tool uses NCEP data. Coll et al. [34] showed a good accuracy of the ACT tool using L7-ETM+ data over full covered rice fields. These authors compared modeled values with ground measurements of LST, showing differences within ±1 K.

Tardy et al. [35] presented the software LANDARTs. This is a Python-based tool using ERA reanalysis product profiles and MODTRAN 4 radiative transfer code [36] to provide pixel-by-pixel atmospheric correction parameters for Landsat images. These authors compared two derived Landsat LST images with two ASTER LST products close in time (time acquisition delay of 13–30 min between platforms). LST differences were mostly within ±2.5 K.

The aim of the present work is to introduce the new Single-Band Atmospheric Correction (SBAC) tool. SBAC provides pixel-by-pixel atmospheric correction parameters regardless of the pixel size. NCEP operational Model Global Tropospheric Analyses (NCEP FNL, ds. 083.2, [37]) are used as inputs. Moreover, as a novelty SBAC accounts for the surface elevation through a digital elevation model (DEM). In this work, we focus on the application of SBAC to L7-ETM+ sensor, although extension to other single channel sensors is automatic and will be explored in further works.

The SBAC tool is introduced in Section 2, together with a description of the study sites and the ground measurements and satellite images used for its validation. Results showing the performance of the SBAC tool, compared to some ACT modeling scenarios, are shown in Section 3. Finally, the main findings of this work are summarized in the conclusion section.

2. Materials and Methods

2.1. Atmospheric Correction in the TIR Domain

The radiance emitted by a surface at temperature (T) with emissivity (ε_i), in a thermal band i is expressed as ε_i B_i(T), where B_i(T) is the blackbody radiation described by Planck’s law. The ‘at’-sensor radiance in band i (L_i) follows the Radiative Transfer Equation (RTE):

$$L_i\left(T_i\right)={\tau }_i\left({\epsilon }_i{B}_i\left(T\right)+\left(1-{\epsilon }_i\right)L_i^{\downarrow }\right)+L_i^{\uparrow }$$

(1)

where T_i is defined as the brightness temperature and corresponds to the temperature for what the Plank’s Law functions equals the ‘at’-sensor radiance. The atmospheric correction parameters τ_i, L^↑_i, L^↓_i, are the atmospheric transmittance, upwelling radiance and the downwelling radiance, respectively. The transmittance has a strong inverse dependence with water vapor content whereas downwelling and upwelling radiances depend on both the atmospheric water vapor and the air temperature vertical profiles [34]. Atmospheric profiles of these two variables, together with some additional information on the aerosol distribution and the viewing geometry, can be used as inputs in radiative transfer codes such as MODTRAN to obtain the atmospheric correction parameters required in Equation (1).

Lambertian behavior in terms of reflectance is commonly assumed in land surface studies. Thus, the downwelling radiance in Equation (1) is defined as L^↓ = (F^↓_hem)/π, where F^↓_hem is the hemispherical downwelling irradiance coming from the surroundings and the atmosphere, which is, in turn, reflected by the surface to the sensor. Garcia-Santos et al. [38] evaluated different methods to estimate this parameter in the field. Galve et al. [17] calculated the downwelling irradiance as:

$$F_{hem}^{\downarrow }={{\displaystyle \int }}_0^{2\pi } d\phi {{\displaystyle \int }}_0^{\raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$2$}\right.} F^{\downarrow }\left(\theta \right) sin\theta cos\theta d\theta$$

(2)

where F^↓ (θ) is the downwelling irradiance at a zenith angle (θ ) since horizontal homogeneity of the atmosphere is assumed, and independent on azimuth angle (ϕ). F^↓ (θ) is simulated with a radiative transfer model at different viewing angles. Wan and Dozier [39] proposed six angles (11.6°, 26.1°, 40.3° and 53.7°), named as Gaussian angles, for their good distribution. For completeness Galve et al. [17] added nadir (0°), 65°, 70°, 80°, 85° and 89° in order to improve the description at larger angles.

Radiance emitted by the surface can be extracted from Equation (1):

$$B_i\left(T\right)=\left(\left(\frac{L_i\left(T_i\right)-L_i^{\uparrow }}{{\tau }_i}\right)-\left(1-{\epsilon }_i\right)L_i^{\downarrow }\right)/{\epsilon }_i$$

(3)

Land surface temperature, T, can be finally obtained by inverting Planck’s function as:

$$T=\frac{k_2}{Ln \left(1 + \frac{k_1}{B_i \left(T\right)}\right)}$$

(4)

where k₁ = 666.09 Wm⁻² sr⁻¹ µm⁻¹ and k₂ = 1282.7 K for Landsat 7-ETM+ band 6. Coll et al. [34] showed that using the band-averaged Equation (3) instead of the spectral version yields and overestimation of LST between 0.15 K and 0.3 K depending on the bandwidth, atmospheric humidity and spectral variations of surface emissivity [40].

As shown in Equation (1), emissivity correction is also required to derive LSTs. When the sensor is provided with more than two channels in the TIR region, the TES method [24] can be used. However, when only one or two TIR bands are available, emissivity values can be assigned or estimated using VNIR information. In this work we used the methodology proposed by [41,42] to obtain emissivity values in vegetated surfaces from the proportion of vegetation cover (Pv):

$${\epsilon }_i={\epsilon }_{i,v} P_v+{\epsilon }_{i,s}\left(1-P_v\right) \left(1-1.74P_v\right)+1.7372 P_v \left(1-P_v\right)$$

(5)

where ε_i,v and ε_i,s are the emissivity assigned to full vegetation cover (0.985), and bare soil (0.960) respectively. P_v was obtained from the Normalized Difference Vegetation Index (NDVI, $i=\left({\rho }_{NIR}-{\rho }_R\right)/\left({\rho }_{NIR}+{\rho }_R\right)$), being ${\rho }_{NIR}$ and ${\rho }_R$ reflectance of near infrared and red bands, as follows:

$$P_v=\frac{\left(1-\raisebox{1ex}{$i$}\!\left/ \!\raisebox{-1ex}{$i_s$}\right.\right)}{\left(1-\raisebox{1ex}{$i$}\!\left/ \!\raisebox{-1ex}{$i_s$}\right.\right)+K\left(1-\raisebox{1ex}{$i$}\!\left/ \!\raisebox{-1ex}{$i_v$}\right.\right)}$$

(6)

where i_s and i_v are the NDVI assigned to bare soil and full vegetation cover of the scene. Duchemin et al. [43] monitored wheat phenology and irrigation in central Morocco and present an exponential relation between LAI and NDVI where extreme values for NDVI are defined. Then infinitely-dense canopy is defined with i = 0.93, and i = 0.14 is defined for a dry bare soil. Asrar et al. [44] determined that NDVI values for extreme bare soils conditions are within 0.1–0.2. In this work, values of i_s = 0.15 and i_v = 0.91 were obtained and assigned to bare soil and full vegetated pixels, respectively. In addition, $K=\left({\rho }_{NIR,v}-{\rho }_{R,v}\right)/\left({\rho }_{NIR,s}-{\rho }_{R,s}\right)$, where ${\rho }_{NIR,j}$ and ${\rho }_{R,j}$ means the reflectance in near infrared and red spectral ranges, respectively, in the full covered pixels (j = v) and bare soil pixels (j = s) and yields values from 2 to 9.

2.2. SBAC Method

The SBAC method is introduced as a new tool that uses the atmospheric profiles extracted from the NCEP as a basis. In particular, the global tropospheric analyses product [36] is the input used to derive the atmospheric correction parameters. NCEP profiles are provided on a 1° × 1° longitude/latitude grid every 6 h (00:00, 06:00, 12:00, and 18:00 UTC). The atmospheric correction magnitudes for a given point and observation time can be obtained by spatial and temporal interpolation from the surrounding NCEP profile sites and the two closest times.

In the SBAC method, the concept is extended to automatically provide the atmospheric correction parameters on a pixel-by-pixel basis for the whole scene, taking as well into account the elevation of each pixel according to a DEM co-registered with the satellite data. For this end, we used the ASTER Global DEM with spatial resolution of 30 m provided by the U.S. Geological Survey’s Land Processes Distributed Active Archive Center (LP DAAC; http://lpdaac.usgs.gov). In the ACT tool [33], NCEP profiles are spatially and temporally interpolated to obtain the profile at the indicated point and time from which the atmospheric transmittance and radiance are derived. However, this procedure is not operational for a pixel-by-pixel atmospheric correction for a full satellite image. The SBAC method works on the basis of the atmospheric correction parameter values calculated for the 1° × 1° NCEP nodes and times, and also accounting for the pixel elevation, as described below. We employed the narrow-band radiative transfer code MODTRAN 5.2 [36,45]. NCEP atmospheric profiles include geopotential height, temperature and relative humidity at 26 fixed pressure levels from 1000 to 10 hPa. Data were downloaded from http://rda.ucar.edu/datasets/ds083.2. A tridimensional grid is stablished (see Figure 1) covering the extent of the sensor scene with a horizontal resolution of 1° × 1° (corresponding to the location of the NCEP data) and 13 altitudes from sea level to the maximum height in the scene (namely z₀ = 0, 50, 100, 150, 200, 300, 500, 750, 1000, 1500, 2000, 3000 and 5000 m). The z₀ values were selected in order to provide enough vertical resolution, with greater detail in the lower layers where atmospheric water vapor is concentrated. NCEP profiles are distributed at fixed pressure levels that are not always associated with a constant altitude. Thus, the NCEP profiles corresponding to each 1° × 1° node have to be modified to achieve a set of 13 adapted profiles in which the different atmospheric levels are adjusted to the 13 predefined altitudes. The modified profiles are generated as follows:

(a)

When the prescribed z₀ value is located between two NCEP layers, the corresponding air temperature and humidity are obtained by linear interpolation from the original NCEP levels, all the upper NCEP layers remaining unchanged and all the layers below the prescribed altitude z₀ being ignored.

(b)

When the prescribed altitude z₀ is below the lowest level of the NCEP profile, we take the NCEP altitude as the lowest z₀ in the 1° × 1° node and do not consider any z₀ value below it.

The 13 modified profiles for each 1° × 1° node are used as inputs in MODTRAN 5 to calculate the spectral atmospheric transmittance and upwelling radiance at nadir, together with downwelling radiance following Equation (2). The spectral magnitudes are then integrated to the thermal infrared band using its spectral response function.

Considering the geolocation (x,y) of a given pixel and its altitude (z, DEM), the SBAC method selects the atmospheric profiles for the 4 closest 1° × 1° nodes and the two closest z₀ altitude levels for the two times closest to the sensor overpass. For the 4 selected 1° × 1° nodes and the two closest z₀ altitudes, the transmittances and radiances at the sensor acquisition time are obtained by linear interpolation between the two closest times. Next, the values corresponding to the x, y coordinates of the pixel at the two closest z₀ altitudes are obtained by inverse distance weighted interpolation from the values for the 4 closest nodes. Thus, the interpolated value for a generic magnitude Q is given by:

$$Q=\frac{{{\displaystyle \sum }}_{k=1 }^4\frac{Q_k}{{\mathrm{d}}_{\mathrm{k}}^2}}{{{\displaystyle \sum }}_{k=1}^4 \frac{1}{d_k^2}}$$

(7)

where Q_k is the value for node k, and d_k is the horizontal distance between the pixel location and node k [46]. Finally, linear interpolation between the two closest z₀ levels yields the atmospheric correction parameters for the given pixel.

In order to be operational, we set a grid with a spatial resolution of 0.01° degrees (~1 km) latitude/longitude and depending on the spatial resolution desired, we resize the atmospheric correction parameters using a linear interpolation.

2.3. Study Sites and Ground Measurements

Two different datasets were used in this work in order to validate the SBAC method, covering a variety of surface and environmental conditions. The first site is located in a semiarid area in central Spain. Ground measurements of LST were taken in “Las Tiesas” experimental farm (39°03′35″N, 2°06′W) during the growing seasons of 2015 and 2016. This has been a traditional ESA (European Space Agency) test site used in different international campaigns: SEN2FLEX (SENtinel-2 and Fluorescence Experiment, [47]), SPARC (SPectra bARrax Campaign, [48]), ImagineS (Implementing Multi-scale Agricultural Indicators Exploiting Sentinels, [49]), and DAISEX (Digital Airborne Imaging Spectrometer EXperiment, [50]). This is a patchy agricultural flat area with an average altitude of 700 m above sea level close to Albacete and is frequently used in agriculture water management studies [51,52,53,54,55]. Ground LST measurements concurrent with 13 Landsat 7-ETM+ overpasses in 2015–2016 were carried out in 6 different crop fields: barley, vineyard, corn, garlic, and almond trees, with the aim of covering a variety of vegetation cover conditions. An overview of the study area and the crop fields is shown in Figure 2. Extension of all selected plots is >10 ha, except the vineyard with ~5 ha, which guaranties at least 3 × 3 Landsat pixels within each plot. The temperatures were measured using three hand-held infrared radiometers (IRTs) Apogee MI-210. These radiometers have a broad thermal band (8–14 μm) with an accuracy of ±0.3 K and 22° field of view. In the sparse crops (vineyard and almond trees), special care was taken in the measurements, averaging soil and canopy temperatures so to get representative values of the target LST. “Las Tiesas” site falls within the overlapped area between Landsat 7 tiles WRS 200/33 and WRS 199/33 (path/row), at a distance to their corresponding scene center of 79 km (North West) and 110 km (North East), respectively.

Ground LST data used by [34] in their validation study were also extracted and incorporated to the present work to reinforce the ground-based dataset and further analysis. These authors focused on the Valencia rice fields test site and took ground measurements of LST in 7 different days (see

Table 1. Dates for the Landsat 7-ETM+ scenes used in this work. Information on the test sites is also included: field crop and extension, location coordinates, and altitude. # refers to the case number for further references.

#	Date	Path	Row	Crop	Extension	Latitude	Longitude	Altitude
					(ha)	(°)	(°)	(km)
1	12/08/2004	198	33	Rice	>100	39.250	−0.295	0.007
2	21/07/2005	199	33	Rice	>100	39.265	−0.308	0.006
3	06/08/2005	199	33	Rice	>100	39.265	−0.308	0.006
4	24/07/2006	199	33	Rice	>100	39.265	−0.308	0.006
5	02/08/2006	198	33	Rice	>100	39.265	−0.308	0.006
6	11/07/2007	199	33	Rice	>100	39.265	−0.308	0.006
7	20/07/2007	198	33	Rice	>100	39.265	−0.308	0.006
8	21/05/2015	200	33	Barley	25.4	39.059	−2.099	0.694
9				Vineyard	5.0	39.060	−2.101	0.692
10	22/06/2015	200	33	Barley	25.4	39.059	−2.099	0.694
11				Corn	25.4	39.059	−2.096	0.693
12	01/07/2015	199	33	Barley	25.4	39.059	−2.099	0.694
13				Corn	25.4	39.059	−2.096	0.693
14				Vineyard	5.0	39.060	−2.101	0.692
15	17/07/2015	199	33	Barley	25.4	39.059	−2.099	0.694
16				Corn	25.4	39.059	−2.096	0.693
17	24/07/2015	200	33	Corn	25.4	39.059	−2.096	0.693
18				Vineyard	5.0	39.060	−2.101	0.692
19	02/08/2015	199	33	Corn	25.4	39.059	−2.096	0.693
20				Vineyard	5.0	39.060	−2.101	0.692
21	09/08/2015	200	33	Corn	25.4	39.059	−2.096	0.693
22				Vineyard	5.0	39.060	−2.101	0.692
23	18/08/2015	199	33	Barley	25.4	39.059	−2.099	0.694
24				Corn	25.4	39.059	−2.096	0.693
25				Vineyard	5.0	39.060	−2.101	0.692
26	01/06/2016	199	33	Garlic	25.4	39.059	−2.099	0.694
27				Garlic	25.4	39.059	−2.096	0.693
28	24/06/2016	200	33	Garlic	25.4	39.059	−2.099	0.694
29				Garlic	25.4	39.059	−2.096	0.693
30				Barley	54.4	39.056	−2.080	0.692
31	19/07/2016	199	33	Vineyard	5.0	39.060	−2.101	0.691
32				Garlic	25.4	39.059	−2.099	0.694
33	26/07/2016	200	33	Almonds	11.1	39.043	−2.090	0.692
34				Vineyard	5.0	39.060	−2.101	0.692
35				Barley	11.1	39.057	−2.088	0.695
36				Wet Bare Soil	25.4	39.058	−2.098	0.689

) in July and August 2004–2007, when rice crops have nearly full cover and are well irrigated. It makes the site highly homogeneous in terms of both surface temperature and emissivity, thus easing the radiometric measurement of LST.

Three TIR radiometers were used in the Valencia test site, two CIMEL CE 312-1 with four bands (8–13, 11.5–12.5, 10.5–11.5, and 8.2–9.2 μm) and one CIMEL CE 312-2 with six bands (8–13, 8.1–8.5, 8.5–8.9, 8.9–9.3, 10.3–11.0, and 11.0–11.7 μm). These instruments show absolute accuracies within a range between 0.1 and 0.2 K in all bands. The absolute calibration for those bands are performed using BB as the reference. Field of view of these sensors is 10°. The Valencia site also falls within the overlapped region between two tiles (WRS 199/33 and WRS 198/33), at a distance to their corresponding scene center of 114 km (North West) and 91 km (North East), respectively.

Ground measurements were only performed under cloud-free sky conditions over the sites. The radiometers were placed on the fields and carried back and forth pointing at nadir view, at a height of 1.5–2 m above the ground surface. Temperatures were taken at a rate of 5–10 measurements/min, covering distances of 30–50 m/min and several hectares in total. Ten-min averages centered at the satellite overpass time were calculated. The standard deviation of the ground temperatures was calculated as a measure of the spatial and temporal variabilities of LST in the test site. It varies depending on the crop field, with values ≤2.0 K, except in 4 cases, and <1 K for the rice field dataset (see

Table 2. Average values of the ground-measured temperatures (T_g). Standard deviation (σ) accounts for the spatial and temporal variability in each plot. Total precipitable water (W) obtained from the 5-km resolution MOD05 product [57]. Satellite brightness temperatures (T_b) and emissivity values corresponding to 3 × 3 pixel averaged for cases labeled in Table 1.

#	Tg	±σ	W	ε	Tb	#	Tg	±σ	W	ε	Tb
	(°C)	(°C)	(cm)		(°C)		(°C)	(°C)	(cm)		(°C)
1	28.2	0.6	2.1	0.988	24.9	19	29.5	0.6	1.8	0.988	25.1
2	28.1	0.4	2.0	0.988	24.4	20	46.7	1.9	1.9	0.978	35.2
3	28.1	0.5	1.9	0.988	24.4	21	28.7	0.3	3.1	0.990	23.1
4	28.8	0.4	2.4	0.988	24.9	22	42.3	0.3	3.1	0.967	31.1
5	29.0	0.9	2.9	0.988	23.9	23	45.5	3.4	2.0	0.975	35.6
6	26.9	0.5	2.9	0.988	21.8	24	28.4	0.3	2.0	0.989	24.5
7	28.0	0.4	2.9	0.988	22.8	25	44.1	1.4	2.0	0.980	35.6
8	19.4	0.9	0.9	0.986	17.8	26	36.6	2.6	0.9	0.990	33.2
9	41.3	2.7	0.9	0.970	34.2	27	29.2	1.8	0.9	0.991	25.2
10	39.0	0.9	2.9	0.977	29.9	28	57.1	2.4	2.2	0.965	44.4
11	31.9	1.1	2.9	0.992	25.5	29	33.8	2.1	2.2	0.987	28.5
12	36.6	0.5	2.0	0.977	31.3	30	37.4	1.3	2.2	0.975	32.8
13	29.0	0.3	2.0	0.989	26.0	31	48.3	0.8	2.2	0.973	38.4
14	44.8	1.3	2.0	0.977	37.1	32	36.4	2.3	2.2	0.969	30.1
15	54.9	1.8	2.0	0.972	42.7	33	50.3	1.6	2.3	0.966	40.5
16	31.2	3.6	2.0	0.989	27.9	34	48.2	1.1	2.2	0.973	39.4
17	34.6	0.4	1.7	0.989	29.0	35	49.9	2.3	2.0	0.968	39.5
18	51.6	0.8	1.7	0.968	41.6	36	31.3	0.8	2.2	0.967	28.1

). Radiometric temperatures were corrected from atmospheric and emissivity effects. Downwelling sky radiance was measured with each radiometer and used for the emissivity correction of the ground-measured LST, including the reflection of the sky irradiance. Emissivity data were obtained by applying the TES procedure [24] to in situ thermal radiance measurements using a radiometer CIMEL CE 312-2 [56]. Broadband emissivity values (8–13 μm) were used.

2.4. Satellite Images

The selection of satellite images was constrained to Landsat 7-ETM+ although the methodology could be applied to any image of the Landsat series, including Landsat 8-TIRS and forthcoming Landsat 9. A total of 19 Landsat scenes were used, with a fixed row (33) and three different paths (198, 199, 200). Figure 3 shows the centers of the scenes and the location of the two study sites over the digital elevation model used.

The Landsat 7 images used belong to Landsat collection Tier 1 with two processing levels, top-of-atmosphere reflectance and surface reflectance. The USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov/) provided the Landsat 7 images.

Radiances in the TIR domain were obtained from band 6 low gain of Landsat 7 ETM+. ETM+ TIR data are acquired at 60-m spatial resolution; however, these data are resampled and provided at 30-m resolution. Landsat from ESPA includes a pixel quality image. Based on this quality image, only free pixels of clouds, snow, shadows and water were kept in this study. All scenes were affected by the failure in SLC (SLC-off mode) that occurred in 2003, so part of the data in the scenes were lost. This issue did not affect data in this work corresponding to the test sites. Table 2 also present the brightness temperature (T_b) and the emissivity (ε) obtained following the procedure mentioned before.

3. Results and Discussion

3.1. Local Validation of SBAC

The ACT developed by [32,33] is being successfully used for the atmospheric correction of Landsat thermal data, as mentioned above. To obtain the atmospheric correction parameter as outputs, the user must provide the ACT with date, hour and location. ACT uses atmospheric profiles from NCEP as inputs of the MODTRAN-4 radiative transfer code to calculate the atmospheric correction parameters for the bandpasses of several Landsat series thermal band (Landsat 5/TM, Landsat 7/ETM+ and Landsat 8/TIRS). Several options allow the user to decide the interpolation process of the atmospheric profiles and select the default upper atmospheric profile. The user is also given the possibility of entering surface atmospheric conditions if known, although in this case the first 3 km of the profile are modified to fit with the ground conditions. Coll et al. [34] tested all options and showed the best results when no ground conditions are imposed, and the profile is obtained from the interpolation of the four surrounding profiles. In this work, this option of the ACT tool was selected to compare with SBAC performance.

The potential of the SBAC tool was assessed in two steps: in situ validation with ground LST measurements and distributed application to full scenes.

A comparison with the results obtained through the ACT was considered in this work to reinforce SBAC validity and feasibility. Based on the similarities between SBAC and ACT, results should coincide, or at least be very close, when the exact location of the target of interest is introduced in the ACT-site. However, differences are expected when the latitude/longitude coordinates of the center of a full Landsat scene are used as the location inputs in the ACT-center. Uncertainties introduced in this case were evaluated in this work.

Figure 4 shows the comparison between ground-measured and satellite derived LST, using the three different configurations described above for the atmospheric correction: SBAC, ACT-site and ACT-center.

All differences are within ±5 K except three cases (23–25). The largest deviations are only observed for the ACT-center atmospheric correction, most likely due to the presence of some clouds in the scene, responsible for the heterogeneity in the atmosphere for that date. For these specific three cases (23–25), LST results using ACT-site shows an underestimation above −2.0 K whereas SBAC overestimates in less than 0.9 K. Average results for SBAC show practically no mean bias (−0.1 K) and small standard deviation (±1.6 K), yielding a root mean square error (RMSE) of ±1.6 K. Note thermal ETM+ data are provided with a noise equivalent temperature <0.3 K [58].

Moreover, differences between observed and estimated LST show no correlation with water vapor content (R² = 0.09) nor brightness temperature (R² = 0.3). Regarding the ACT tool, after removing the cases 23–25, similar RMSE values close to ±2 K were observed in both cases, with an underestimation of −1.1 K for ACT-center and −0.7 K for ACT-site. Again, poor correlation with water vapor content and brightness temperature (R² ≤ 0.1 and R² ≤ 0.4,) was observed in both cases.

Deeper analysis implies to compare the atmospheric correction parameters obtained from the different sources. Quantitative results are shown in Figure 5 where transmittance, upwelling radiance and downwelling radiance obtained with both, ACT-site and ACT-center, are plotted against those obtained from SBAC.

Very good agreement, close to 1:1 line, is observed between SBAC and ACT-site outputs, whereas significant deviations are obtained with ACT-center.

The effect of the reflected atmospheric radiance on the LST differences is small in this work due to the high emissivity values considered and the fact that the same emissivity value was used for the three configurations of the atmospheric correction. In order to evaluate the impact of τ and L^↑_i on LST retrieval, Coll et al. [59] defined the atmospheric correction term as the term to be added to ‘at’-sensor brightness to compensate the atmospheric absorption and emission:

$$\Delta T_{atm}=\frac{1-\tau }{\tau }\left(T_b-T_{eff}\right)$$

(8)

where T_eff is the effective atmospheric temperature, which represents the average temperature of the atmosphere in terms of the emitted radiance and is defined from atmospheric transmittance and upwelling radiance as:

$$L^{\uparrow }=\left(1-\tau \right)B\left(T_{eff}\right)$$

(9)

from which T_eff can be inverted using (4). Since the major absorption and emission of the atmosphere is due to water vapor content, T_eff can be related to the averaged air temperature weighted with the water vapor mixing ratio profile, so T_eff is closer to the mean air temperature at low levels.

In (8) we can find two main factors, the transmittance factor (TF = ((1 − τ))/τ) and the difference between brightness temperature and T_eff. Differences between TF obtained by ACT and SBAC are consequence not only of the dependence of the atmospheric transmittance on water vapor content but also on the vertical distribution of water vapor and temperature, due to the temperature dependence on the water vapor continuum absorption [59]. In fact, the difference in atmospheric correction shows no correlation with water vapor content as seen in Figure 6a. However, a clear correlation appears between differences in the atmospheric correction, and their TF differences (Figure 6b). Focusing on the effective atmospheric temperature, ACT yields values higher than SBAC (Figure 6c). The use of SBAC appears to be equivalent in this work to the use of ACT-site, with an overestimation in transmittance around 2% and an underestimation for the upwelling radiance lower than 4%. The comparison with the use of the center of scene coordinates (ACT-center) yields a similar overestimation below 2% in transmittance but a great underestimation in upwelling radiance higher than 13% in this case. These discrepancies increase with the heterogeneity of the atmosphere within.

The scene. For instance, in cases 23–25 the presence of clouds yields differences in transmittance and upwelling radiance of 30% and −85%, respectively, whereas SBAC and ACT-site differ by −2% and 1.2%, respectively.

3.2. Full-Scene Analysis of SBAC Tool Outputs

Beyond the comparison at a specific location, the real potential of the SBAC tool, compared to other atmospheric correction tools, comes when atmospheric correction is conducted at a full-scene scale. LST images obtained from the SBAC tool were compared to those using the ACT applied to the center of the scene coordinates in this case (ACT-center). As an example, two L7 scenes for path 200 and row 33 were selected for two days with different atmospheric conditions, 21 May 2015 and 24 July 2015. Figure 7 shows the difference in LST between SBAC and ACT-center (dT). Images of the total precipitable water content (TPWC), from the MOD05 [57] product, were also included in Figure 7, to illustrate the atmospheric condition distribution for both days. Mean TPWC values of 1.1 ± 0.3 cm and 1.8 ± 0.2 cm, are observed for 21 May and 24 July, respectively.

In order to discard from the analysis pixels that might be affected by surrounding effects of clouds or shadows, or those affected by the SLC failure, a thermal homogeneity filter was applied. A pixel is considered valid if the center-located 3 × 3 array shows a variance <2.5 K and is labeled as clear in the cloud mask. LST differences in the full scene are mainly within the range 0.5–3.0 K, with a small (less than 9% in 21 May and less than 3% in 24 July) percentage out of it. Average LST overestimations are 0.9 K in 21 May and 1.3 K in 24 July, with standard deviation of ±0.5 K and ±0.6 K, respectively.

Figure 8 shows the spatial distribution of the atmospheric correction within the scenes in terms of the transmittance factor obtained by the pixel-by-pixel SBAC results. Note the spatial deviations of this parameter, ranging 0.10–0.30, in comparison to the fixed TF values of 0.11 and 0.19 extracted from the ACT-center, in the 21 May and 24 July scenes, respectively. Lower values of TF and more homogeneity are observed in May due to less atmospheric water vapor content compared to July.

Largest differences are observed in the southeastern sector of the images, with LST differences over 3.0 K in some cases. According to the DEM shown in Figure 3, this area corresponds to some mountains contrasting with the flat plain dominating the scenes. However, sole dependence of the LST deviations on altitude or orography is not so evident. For example, significant LST differences observed in the northeastern borders of the July image denote further influences. For a more in-depth analysis, some linear profiles of dT were extracted from the images in Figure 7 and plotted together with the corresponding DEM profiles versus distance to the image center. Figure 9 shows two representative examples corresponding to the crossing straight diagonal lines in the scenes. Note pixels close to the center of the scene, within a radius of 50 km, show overall lower dT values, set in 0.6 K and 0.8 K in 21 May and 24 July, respectively, with standard deviation of ±0.11 K in both cases.

Plots in Figure 9 indicate there is not a direct effect of the distance to the center. For example, in the east half-side of the May (July) scene average overestimations of 1.8 K (1.1 K), with standard deviations of ±0.8 K (±0.3 K), are obtained, whereas in the west half-side similar results for both scenes are obtained, an overestimation of 0.9 K with a standard deviation of ±0.2 K. Differences appear rather as a consequence of variations in terrain altitude. Lower or higher elevations imply different atmospheric layer, compared to the center of the scene, and then the SBAC tool better captures variations in the corresponding atmospheric conditions.

However, the largest dT values are obtained in the lower-right corner of the scene, as mentioned above, under similar elevation conditions, showing that atmospheric conditions are the dominant factor for these differences. Also, note dT values are generally higher in July than in May, and also variability is stressed in the first case responding to higher TF values.

4. Conclusions

Results in this work give us the first insight on the feasibility and potential of the SBAC tool to provide pixel-by-pixel atmospheric correction parameters and then estimate LST values from single-channel sensors. We focus on L7-ETM+ data although it could be further applied to L8-TIRS or any other satellite thermal sensor.

No significant differences in retrieved LST are observed between SBAC and well-known ACT tools if the atmospheric conditions remain stable with no appreciable variations throughout the scene. SBAC results are in very good agreement with those obtained from ACT when the study site coordinates are specified. Local validation, using ground data from two different study sites covering a variety of surface and atmospheric conditions, shows an average estimation error in LST < 2.0 K using both SBAC and ACT tools, with negligible systematic deviation for SBAC and an underestimation ~1.0 K using ACT. Moreover, the atmospheric correction provided by SBAC and ACT-site are equivalent in terms of transmittance factor. In general, SBAC yields lower values of effective atmospheric temperature than ACT. Differences in atmospheric correction show no significant correlation with water vapor content.

The analysis at a full-scene scale yields no significant dependence of LST differences between SBAC and ACT on the distance of the study site to the coordinates inserted in the ACT tool, but rather on the orography variations, especially for lower altitudes. These differences are also powered by spatial variations in the atmospheric water vapor content. Average LST differences close to 1.0 K were obtained for the two full L7 images used as an example, with punctual differences reaching 3.0 K between SBAC and ACT tools.

In short, the SBAC tool automatically provides a full image pixel-by-pixel atmospheric correction, accounting for the surface elevation, with a precision in very good agreement with that running the ACT tool to every single pixel within the scene. Implementing the SBAC tool in satellite thermal image processing will improve the accuracy in the spatial estimation of land surface temperature. This might boost the exploitation of single-channel satellite thermal sensors in a variety of environmental applications.