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 +

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.


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.

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): 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: 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 Land surface temperature, T, can be finally obtained by inverting Planck's function as: where k 1 = 666.09Wm −2 sr −1 µm −1 and k 2 = 1282.7K 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): 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 = (ρ N IR − ρ R )/(ρ N IR + ρ R )), being ρ N IR and ρ R reflectance of near infrared and red bands, as follows: 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, , where ρ N IR,j and ρ 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.

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 = 0, 50, 100, 150, 200, 300, 500, 750, 1000, 1500, 2000, 3000 and 5000 m).The z 0 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 0 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 0 being ignored.(b) When the prescribed altitude z 0 is below the lowest level of the NCEP profile, we take the NCEP altitude as the lowest z 0 in the 1 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 z0 altitude levels for the two times closest to the sensor overpass.For the 4 selected 1° × 1° nodes and the two closest z0 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 z0 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: where Qk is the value for node k, and dk is the horizontal distance between the pixel location and node k [46].Finally, linear interpolation between the two closest z0 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.

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 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 0 altitude levels for the two times closest to the sensor overpass.For the 4 selected 1 • × 1 • nodes and the two closest z 0 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 0 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: 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 0 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.

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) 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.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) 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).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.
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.

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.
This issue did not affect data in this work corresponding to the test sites.Table 2 also present the brightness temperature (Tb) and the emissivity (ε) obtained following the procedure mentioned before.

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.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.

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.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 19 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)(24)(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)(24)(25), LST results using ACT-site shows an underestimation above −2.0K 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].  1. Dashed line marks the 1:1 line.Moreover, differences between observed and estimated LST show no correlation with water vapor content (R 2 = 0.09) nor brightness temperature (R 2 = 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 2 ≤0.1 and R 2 ≤ 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.  1. Dashed line marks the 1:1 line.
All differences are within ±5 K except three cases (23)(24)(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)(24)(25), LST results using ACT-site shows an underestimation above −2.0K 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 2 = 0.09) nor brightness temperature (R 2 = 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 2 ≤ 0.1 and R 2 ≤ 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.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: where Teff 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: from which Teff can be inverted using (4).Since the major absorption and emission of the atmosphere is due to water vapor content, Teff can be related to the averaged air temperature weighted with the water vapor mixing ratio profile, so Teff 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 Teff.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 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: 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: 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.

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.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.

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.0K, 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.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.0K, 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.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.11K 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.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.11K 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.

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 singlechannel 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.

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.

Figure 1 .
Figure 1.Scheme of the tridimensional grid cell defined by the SBAC method.A-D are 1 • × 1 • nodes at altitude z 0 and E-H are 1 • × 1 • nodes at altitude z 0 + ∆z.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 19 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 Figure2.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.

Figure 2 .
Figure 2. Location and overview of "Las Tiesas" experimental farm and the Valencia validation site, together with some pictures of the different crop fields.

Figure 2 .
Figure 2. Location and overview of "Las Tiesas" experimental farm and the Valencia validation site, together with some pictures of the different crop fields.

Figure 3 .
Figure 3. Digital elevation map of the southeastern region of Spain."Las Tiesas" and "Valencia" site locations are spotted.The geographic center of the different Landsat 7-ETM+ path and row used in the validation are also highlighted.

Figure 3 .
Figure 3. Digital elevation map of the southeastern region of Spain."Las Tiesas" and "Valencia" site locations are spotted.The geographic center of the different Landsat 7-ETM+ path and row used in the validation are also highlighted.

Figure 4 .
Figure 4. Ground-measured LST versus satellite LST derived from ACT (site and center) and SBAC, for all cases displayed in Table1.Dashed line marks the 1:1 line.

Figure 4 .
Figure 4. Ground-measured LST versus satellite LST derived from ACT (site and center) and SBAC, for all cases displayed in Table1.Dashed line marks the 1:1 line.

Figure 6 .
Figure 6.(a) Differences between LST obtained from SBAC and ACT (-site and -center) versus water vapor content (W).(b) Differences between LST obtained from SBAC and ACT (-site and -center) versus the transmittance factor (TF) difference between the same models.(c) Atmospheric effective temperature (Teff) obtained from ACT (-site and -center) versus SBAC.Dashed lines show the 1:1 line.

Figure 6 .
Figure 6.(a) Differences between LST obtained from SBAC and ACT (-site and -center) versus water vapor content (W).(b) Differences between LST obtained from SBAC and ACT (-site and -center) versus the transmittance factor (TF) difference between the same models.(c) Atmospheric effective temperature (T eff ) obtained from ACT (-site and -center) versus SBAC.Dashed lines show the 1:1 line.

Figure 7 .
Figure 7. Maps of LST differences between SBAC and ACT-center (dT) for two scenes (WRS 200/33) corresponding to dates 21 May 2015 (left) and 24 July 2015 (right).White pin and black triangle locate the scene center and study site, respectively.Histograms of the LST differences are shown on top.Lower images show corresponding TPWC maps (cm).

Figure 7 .
Figure 7. Maps of LST differences between SBAC and ACT-center (dT) for two scenes (WRS 200/33) corresponding to dates 21 May 2015 (left) and 24 July 2015 (right).White pin and black triangle locate the scene center and study site, respectively.Histograms of the LST differences are shown on top.Lower images show corresponding TPWC maps (cm).
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.Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 19 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.

Figure 8 .
Figure 8. SBAC Transmittance Factor (TF) for the same scenes as in Figure 7: 21 May 2015 (left) and 24 July 2015 (right).TF for the ACT (center) is marked in the color ramp for both cases.

Figure 8 .
Figure 8. SBAC Transmittance Factor (TF) for the same scenes as in Figure 7: 21 May 2015 (left) and 24 July 2015 (right).TF for the ACT (center) is marked in the color ramp for both cases.

Figure 9 .
Figure 9. Profiles of dT (LST SBAC-LST ACT-center) versus distance to center of the image.Data were extracted from crossing straight diagonal lines: upper-left to lower-right corners (up), and lower-left to upper-right corners (down).Elevation of each pixel, referred to the value in the center of the image (0.941 km), is also plotted.

Figure 9 .
Figure 9. Profiles of dT (LST SBAC-LST ACT-center) versus distance to center of the image.Data were extracted from crossing straight diagonal lines: upper-left to lower-right corners (up), and lower-left to upper-right corners (down).Elevation of each pixel, referred to the value in the center of the image (0.941 km), is also plotted.
• × 1 • node and do not consider any z 0 value below it.

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.