A Software Tool for Atmospheric Correction and Surface

Preparing to load PDF file. please wait...

0 of 0
A Software Tool for Atmospheric Correction and Surface

Transcript Of A Software Tool for Atmospheric Correction and Surface

Technical Note
A Software Tool for Atmospheric Correction and Surface Temperature Estimation of Landsat Infrared Thermal Data
Benjamin Tardy 1, Vincent Rivalland 1,*, Mireille Huc 1, Olivier Hagolle 1, Sebastien Marcq 2 and Gilles Boulet 1
1 CESBIO, Université de Toulouse, CNES/CNRS/IRD/UPS, 18 av. Edouard Belin, bpi 2801, F-31401 Toulouse Cedex 9, France; [email protected] (B.T.); [email protected] (M.H.); [email protected] (O.H.); [email protected] (G.B.)
2 CNES, 18 av. Edouard Belin, F-31400 Toulouse, France; [email protected] * Correspondence: [email protected]; Tel.: +33-5-61-55-85-06
Academic Editors: Dongdong Wang, Richard Müller and Prasad S. Thenkabail Received: 25 March 2016; Accepted: 2 August 2016; Published: 24 August 2016
Abstract: Land surface temperature (LST) is an important variable involved in the Earth’s surface energy and water budgets and a key component in many aspects of environmental research. The Landsat program, jointly carried out by NASA and the USGS, has been recording thermal infrared data for the past 40 years. Nevertheless, LST data products for Landsat remain unavailable. The atmospheric correction (AC) method commonly used for mono-window Landsat thermal data requires detailed information concerning the vertical structure (temperature, pressure) and the composition (water vapor, ozone) of the atmosphere. For a given coordinate, this information is generally obtained through either radio-sounding or atmospheric model simulations and is passed to the radiative transfer model (RTM) to estimate the local atmospheric correction parameters. Although this approach yields accurate LST data, results are relevant only near this given coordinate. To meet the scientific community’s demand for high-resolution LST maps, we developed a new software tool dedicated to processing Landsat thermal data. The proposed tool improves on the commonly-used AC algorithm by incorporating spatial variations occurring in the Earth’s atmosphere composition. The ERA-Interim dataset (ECMWFmeteorological organization) was used to retrieve vertical atmospheric conditions, which are available at a global scale with a resolution of 0.125 degrees and a temporal resolution of 6 h. A temporal and spatial linear interpolation of meteorological variables was performed to match the acquisition dates and coordinates of the Landsat images. The atmospheric correction parameters were then estimated on the basis of this reconstructed atmospheric grid using the commercial RTMsoftware MODTRAN. The needed surface emissivity was derived from the common vegetation index NDVI, obtained from the red and near-infrared (NIR) bands of the same Landsat image. This permitted an estimation of LST for the entire image without degradation in resolution. The software tool, named LANDARTs, which stands for Landsat automatic retrieval of surface temperatures, is fully automatic and coded in the programming language Python. In the present paper, LANDARTs was used for the local and spatial validation of surface temperature obtained from a Landsat dataset covering two climatically contrasting zones: southwestern France and central Tunisia. Long-term datasets of in situ surface temperature measurements for both locations were compared to corresponding Landsat LST data. This temporal comparison yielded RMSE values ranging from 1.84 ◦C–2.55 ◦C. Landsat surface temperature data obtained with LANDARTs were then spatially compared using the ASTER data products of kinetic surface temperatures (AST08) for both geographical zones. This comparison yielded a satisfactory RMSE of about 2.55 ◦C. Finally, a sensitivity analysis for the effect of spatial validation on the LST correction process showed a variability of up to 2 ◦C for an entire Landsat image, confirming that the proposed spatial approach improved the accuracy of Landsat LST estimations.

Remote Sens. 2016, 8, 696; doi:10.3390/rs8090696


Remote Sens. 2016, 8, 696

2 of 24

Keywords: land surface temperature; Landsat; software tool; atmospheric correction; thermal infrared remote sensing; emissivity

1. Introduction
Land surface temperature (LST) and land surface emissivity (LSE) are two key parameters used in many environmental studies because they are closely connected to the Earth’s surface energy balance. At a regional scale, land and sea surface temperatures (LST and SST, respectively) can be estimated with the help of thermal infrared (TIR) measurements recorded by remote sensors onboard spaceborne platforms. These estimates are used for hydrological and meteorological purposes. In particular, they provide valuable information to monitor evapotranspiration and determine the soil water status on agricultural lands [1–3]. To study the oceans, SST data provide information on sea warming, potential evaporation and water circulation [4]. However, atmospheric effects of absorption and emission affecting the thermal spectrum must first be removed in order to use thermal imagery in absolute temperature research application [5]. At present, very few surface temperature data products obtained through remote sensing observations are available. The MODIS products distributed by the U.S. National Aeronautics and Space Administration (NASA) are the best known. The most commonly-used global and daily land surface temperature and emissivity product is MOD11A1 [6]. MOD11A1 (Collection 5), which uses a grid spatial resolution of 1000 m, is based on the generalized split-window LST algorithm applied to the MODIS thermal multispectral bands [7]. Note that recent works proposed an operational Kalman filter strategy applied to the 3 IR SEVIRI (Spring Enhanced Visible and Infrared Imager) channels from geostationary platform plus ECMWF meteorological analysis (used also here) to directly retrieve surface temperature and the three channels emissivities [8,9]. Spatial resolution of products is around 3 km but earth is fully scanned every 15 min, and the retrieved surface emissivity and temperature accuracy is of ±0.005 and ±0.2 ◦C, respectively. nary platform plus ECMWF meteorological analysis (used also here) to directly retrieve surface temperature and the three channels emissivities [8,9]. Spatial resolution of products is around 3 km but earth is fully scanned every 15 min, and the retrieved surface emissivity and temperature accuracy is of ±0.005 and ±0.2 ◦C, respectively. Proposed by the Ministry of Economy, Trade and Industry (METI) of Japan in collaboration with NASA, the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) kinetic temperature and emissivity products [10] offer a higher spatial resolution of 90 m, but with a low frequency of revisits; indeed data are potentially available only one to eight times per year.
Landsat satellites have been continuously recording multispectral images of Earth since the early 1970s. The U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center provides these data free of charge. In fact, from the first Landsat 1 (launched in 1972) to the current Landsat 8 (launched in 2013), the project has made available more than 40 years of Earth land surface data to study land processes. In total, this represents more than one million scenes, each one offering at least one thermal band. Landsat 5 and 7 provide one thermal channel; Landsat 8 provides two, although it is recommended to use only one band [11,12].
In this sense, the Landsat program appears to be an interesting source of LST records. However, land and sea surface temperatures have not been made available in a free-to-use format. Nevertheless, several methods have been proposed to estimate LST from Landsat’s single TIR channel [13–15] by using atmospheric profiles obtained from databases or through radio sounding. Several studies [16,17] have validated the application of such methods. For example, Barsi [18] proposed a web-based correction tool that used modeled atmospheric global profiles obtained from the National Center for Environmental Prediction (NCEP) in conjunction with the commercially available radiative transfer model MODTRAN [19]. This process allowed the authors to obtain the necessary atmospheric transmission, upwelling, and downwelling radiance values to convert

Remote Sens. 2016, 8, 696

3 of 24

sensor radiance data into surface brightness temperatures (Tb). A subsequent study [20] used the same process to successfully reproduce the results with the North American Regional Reanalysis (NARR) database by using atmospheric profiles as inputs. Note that in both approaches, the computed atmospheric parameters are representative of a region of interest defined by the spatial resolution of the atmospheric profiles. Moreover, in both cases, raw user manipulations are necessary to estimate brightness temperatures, emissivities, and LST from Landsat image radiances. In most studies that used the Landsat thermal channel, LST were corrected for local atmospheric parameter values obtained from [21], or through radio sounding in conjunction with a radiative transfer model such as MODTRAN, or with a single channel method [22]. As a result, the correction process is more adapted to local than to regional applications. For the latter, radio sounding could prove to be inappropriate because it may be either non-coincidental or too far away from the image footprint, leading to errors such as discussed by Mira et al. [23]. The present study expands on Barsi’s approach and proposes an independent and open source software tool that uses the European ERA-Interim global atmospheric reanalysis database in conjunction with the MODTRAN software package. The proposed software tool, named LANDARTs, which stands for Landsat automatic retrieval of surface temperature, is able to automatically estimate surface brightness temperatures and LST for an entire Landsat thermal image with a coherent spatial correction.
Estimating LST through TIR radiance measurements obtained with remote sensors usually requires four steps of correction [5]: (1) converting spectral radiances into at-sensor brightness temperatures, which is typically done by inverting Planck’s law adapted to the spectral sensor window as described by Equation (4); (2) correcting for atmospheric absorption and re-emission, which requires knowledge of the atmospheric profile composition, especially concerning the water content and temperature; (3) correcting for surface emissivity (estimated if unknown); and (4) correcting for surface roughness or topography. This last correction is rarely taken into account at resolutions of 100 m and is neglected at the nadir view angle. Neglecting atmospheric correction results in systematic errors that are proportional to the target temperature for a given atmosphere.
Various approaches have been developed to solve or simplify the atmospheric correction process for TIR data dedicated to single-channel sensors. The available methods can be classified into four categories: the mono-window algorithm, which requires air temperature data only [17]; the single-channel algorithm, which requires total atmospheric water content values only [13]; the look-up tables approach, which supposes an exhaustive learning database [14]; and the physical approach, which requires the calculation of atmospheric parameters with a radiative transfer model and demands detailed knowledge of the atmospheric profile [15,18]. The main goal of our paper is to present and evaluate an automatic Python tool that relies on the MODTRAN model and on ERA-Interim data. The tool was designed to perform atmospheric and surface correction of Landsat (TM and TIRS) thermal data sensors and to retrieve LST and LSE values at a resolution of 30 m (the resolution of radiances provided by the USGS). First, we present the theory of atmospheric correction and the applicable processing algorithm; We then describe the process of comparing the LST reference product from ASTER with in situ temperature measurements to validate the tool’s performance; Finally, we discuss the sensitivity of the results to spatial variations in the atmospheric correction parameters.

2. Materials and Methods

2.1. Atmospheric Correction The spectral distribution of the radiation from a blackbody is described by Planck’s law:

B (T) =








C2 λT



Remote Sens. 2016, 8, 696

4 of 24

where Bλ(T) is the spectral radiance (W· m−2· sr−1 · µm−1) of a blackbody at temperature T (K) and wavelength λ (µm−1); C1 and C2 are two physical constants (C1 = 1.191 × 108 W·µm−4·sr−1·m−2, C2 = 1.439 × 104 µm·K). However, land surfaces measured by remote sensing are gray bodies with
a surface emissivity ( s < 1). Therefore, radiance measured by a sensor of a gray body without
atmospheric perturbation can be expressed as:

Lλ,sen,no−atmo(T) = s,λ Bλ(T)


where Lλ,sen,no−atmo(T) (W· m−2· sr−1 · µm−1) is the spectral radiance of a gray body at temperature T (K) measured by a sensor without atmosphere, s,λ is the surface spectral emissivity at wavelength λ and Bλ(T) (W· m−2· sr−1 · µm−1) is the spectral radiance of the equivalent black-body at temperature T (K) as detailed in Equation (2).
When processing the thermal data of sensor such Landsat (Band 6 on L5 and L7, band 10 on L8; see
Table 1), we need to take into account the atmospheric interaction and use a radiative transfer model
to remove the atmospheric effects. In this goal we need to solve the radiative transfer Equation (3) in
the band wavelength window of the sensor. A radiative transfer model can be used to estimate the
unknowns terms provided the corresponding vertical atmospheric composition is known.

Lλ,sen(T) = τλ s,λ Bλ(T) + Lu,λ,atmo + τλ(1 − s,λ)Ld,λ,atmo


where: Lλ,sen = the spectral at-sensor radiance (top of the atmosphere) (W· m−2· sr−1 · µm−1) Bλ = the radiance of a supposed blackbody surface target at a kinetic temperature T in (K) τλ = the atmospheric transmittance (unitless) at the wavelength λ (µm) Lu,λ,atmo = the upwelling atmospheric radiance in the wavelength window (W· m−2· sr−1 · µm−1) Ld,λ,atmo = the downwelling atmospheric radiance in the wavelength window (W· m−2· sr−1 · µm−1) s,λ = the surface spectral emissivity (unitless).
As a result, the derived spectral atmospheric parameters (τ, Lu and Ld) and the surface emissivity ( s) permit the conversion of the at-sensor radiance into surface radiance exempt from atmospheric
effects. The obtained surface radiance can then be converted into surface kinetic temperatures (LST) by
inverting a simplified Planck’s law from Equation (1) adapted to specific sensor spectral window by
spectral integration according to [24]:




ln (1 + BK(T1 ) )

where: LST = equivalent to T, the land surface temperature (K) B(T) = the surface radiance (W· m−2· sr−1 · µm−1) K1 = the pre-launch calibration constant 1 (W· m−2· sr−1 · µm−1) K2 = the pre-launch calibration constant 2 (K) K1 and K2 are provided by the USGS for L5 and L7 [24], and for L8 [25].
Among these key parameters, the transmittance is known to be inversely proportional to the downwelling and upwelling radiance (see the second and third figures of Section 4.3). The evolution of these atmospheric parameters is primarily controlled by the total water vapor content in the vertical atmospheric column of air [17] and secondary by air temperature profile. The total water vapor content changes from hour to hour, but follows a more or less annual trend depending on latitude. This trend generally attains its maximum in summer (warm and humid atmosphere) and its minimum in winter (cold and dry atmosphere). This atmospheric impact on the surface temperature estimation, that we can associate to the difference between top of atmosphere (TOA) and top of canopy (TOC) brightness temperature (respectively Tb,TOA and Tb,TOC), is proportional to the

Remote Sens. 2016, 8, 696

5 of 24

target temperature. The warmer the surface target temperature is, the higher the atmospheric impact, therefore the correction. Below a temperature limit, around the atmospheric temperature, the mean lower layer atmospheric effect is inversed and Tb,TOA became higher than Tb,TOC. For the two extremes: hot and cold targets, suppose brightness temperature equivalent to LST correspond to under and over estimation of LST, respectively.
The atmospheric transmittance is a very sensitive parameter for LST estimation. Qin et al. [17] shows that an error of 0.05 on transmittance value generated an error comprised between 1–4 ◦C on LST estimation. In addition to daily and annual variations, the existence of spatial variability in the atmospheric composition, which results in spatial variation of the atmospheric correction parameters, can cause large errors if a constant value per image is used.

2.2. Radiative Transfer Modeling with MODTRAN
MODTRAN (moderate resolution atmospheric transmission) is a commercial atmospheric radiative transfer model developed by the U.S. Air Force [19,26]. The model covers a spectrum of 0.2–100 µm that is from the mid-ultraviolet to the far infrared. MODTRAN can be used to estimate the three atmospheric correction parameters in the Landsat thermal spectrum. Because MODTRAN cannot separate downward contributions from upward contributions in a single run, the model was run twice in nadir observation configuration (according to Landsat acquisition mode) using the same atmospheric profile with the following two configurations:
1. The upwelling radiance (i.e., the atmospheric effect between the surface target and the sensor) and the atmospheric transmittance between the surface and the sensor were estimated through simulation by setting the sensor location to 100 km above the surface (considered the sensor altitude) and setting both the surface albedo and emissivity to zero, corresponding to a complete lack of surface reflection for the entire spectrum.
2. The downwelling radiance was estimated through a second run using a configuration in which the sensor was located at 1 m above the surface and the surface albedo was set to 1. For this configuration, we assumed that the downwelling radiance was completely reflected upward toward the sensor.
As the MODTRAN outputs are provided following model spectral resolution (0.1 cm−1), an integration was performed between the bounds defined by the sensor thermal spectrum response .

λi,max var(λ)Rs(λ)dλ

var(λi) = λi,min


λλi,im,minax Rs(λ)dλ

where Rs is the spectral response (unitless) of the sensor at center wavelength λi (µm) of a band with spectral window from λi,min to λi,max. var(i) alternatively denotes the transmittance τ, the upwelling radiance Lu, or the downwelling radiance Ld value extracted from the corresponding column of the MODTRAN output files, between the lower λi,min and upper λi,max wavelengths.

2.3. Atmospheric Profiles Construction
The ECMWF (European Centre for Medium-Range Weather Forecasts) provides full access to several atmospheric reanalysis datasets. In accordance with our requirements, we selected the ERA-Interim global atmospheric reanalysis [27], which provides a variety of vertical and surface atmospheric fields from 1979 to the present (accessible with a 2-month delay) in uniform latitude/longitude grids (0.125◦, 0.25◦, 0.75◦, 1◦, 1.125◦, 1.5◦, 2◦, 2.5◦ and 3◦). The parameters are interpolated from the original N128 reduced Gaussian grid using bilinear methods. LANDARTs integrates the Python API provided by ECMWF to process automatic requests sent to the database [28]. For the Landsat LST data atmospheric correction process, each request sent to the ERA-Interim dataset is defined to cover a footprint that is larger than the image. The returned query

Remote Sens. 2016, 8, 696

6 of 24

results are 2D or 3D matrices of surface and vertical atmospheric fields with a default horizontal grid resolution of 0.125◦, which represents a grid of about 18 × 28 = 504 nodes for each Landsat image. Vertical fields are provided for 37 atmospheric pressure levels from 1–1000 hPa, which corresponds to approximate altitudes of between 100 m and 35 km, depending on the air temperature. In addition to pressure, we were interested in the geopotential, the temperature, and the relative humidity. The altitude of each level can be determined correctly by dividing the geopotential by the gravity constant (g = 9.81 s−2), the air temperature is given in Kelvin, and the relative humidity is given in percent. Analogous to Barsi’s atmospheric correction parameter calculator (Atmocorr, http://atmcorr.gsfc.nasa.gov/), fields from the highest atmospheric levels (for altitudes between 35 km and 100 km) were obtained from the U.S. standard atmosphere of 1976 [29]. For the northern hemisphere, the atmosphere’s summer mid-latitude profile was used for images acquired between March and September and the winter profile was used for all other images; the inverse date ranges were used for the southern hemisphere. At this step, all the vertical profiles for each of the ERA-Interim fields of interest were described using 1000 hPa as the starting value, regardless of the landform or elevation. The surface elevation and air temperature were provided by the surface field data product. However, relative humidity values were not available in the surface field, we estimated it from 2 m dew-point temperature and the surface pressure according to Equation (6) from [30]:

e(t, P)

RH = 100 es(t, P) , with es(td) = e(t)


where RH (%) is the relative humidity, t the air temperature, td the dew-point air temperature, P the atmospheric pressure, es the air saturation vapor pressure and e the air vapor pressure estimated with Magnus formulae [31]. The two datasets were merged to create a consistent meteorological vertical altitude profile extending from the Earth’s surface to the top of the atmosphere. In practice, for surface field pressures higher than or equal to 1000 hPa (first vertical level), surface fields were added at the first level of the vertical profile; for surface field pressures less than 1000 hPa (occurring mainly because of topographic features), surface fields were added to vertical profiles according to the pressure value, and all field profiles below these levels were cut off. For each case, pressure and altitude values were verified to ensure a strictly decreasing order of the vertical profile. However, the values of the vertical profile were not smoothed through interpolation. In fact, failure to respect this decreasing order created a conflict in the MODTRAN execution. Finally, the number of levels in the vertical profile may change depending on the surface elevation. The resulting vertical profile was then assigned to the MODTRAN input file. MODTRAN auto-completes the standard atmospheric profile for standard gases. The ERA-Interim profiles were adapted for each point of the requested grid.

2.4. Emissivity Estimation
Land surface emissivity (LSE) is a key parameter that describes the radiative absorption power of a surface in the longwave radiation spectrum. LSE depends on the target surface top layer composition, such as presence of soil, soil type, vegetation and density, or roughness of the surface [32,33]. Although LSE data are indispensable, they have thus far proven difficult to estimate with remote sensing applications. The sensitivity analysis performed by Olioso [34] showed that an error of ±0.01 for emissivity results in LST values with an error between 0.6 ◦C and 0.9 ◦C (depending on the LST). Qin et al. [17] reported the same range of values. The availability of emissivity data suitable for LST corrections also depends on the sensor spectral windows and on the viewing angles [35]. Landsat spectral windows are given in Table 1. Data were acquired at the nadir, thus avoiding the need for angular corrections. Several solutions exist to retrieve emissivity from satellite data [33]: they may be based (1) on the land cover classification and emissivity coverage tables; (2) on the exploitation of multi-spectral channels, as is the case for temperature and emissivity separation (TES) methods; or (3) on empirical relationships involving vegetation indices, which is the most commonly used method. Li et al. [36] compared six emissivity extraction methods and concluded that all six

Remote Sens. 2016, 8, 696

7 of 24

methods yielded similar results. The first method required a classification of the zone of interest, which is not always easy to obtain; The second method required at least three thermal channels, which are unavailable for Landsat 5 and 7, and the USGS recommends using only thermal band 10 for L8. Therefore, we ultimately decided to rely on the vegetation index method. Most vegetation index methods involve the normalized difference vegetation index (NDVI). The NDVI is defined as:

NDV I = N IR − RED (7) NIR + RED

where NIR represents the reflectance in the near infrared region, and RED is the spectral reflectance in the red region. A commonly used method for calculating mixed pixel emissivity is the estimation proposed by Sobrino et al. [16], adapted from [37]. This empirical relationship is based on an exponential function that is NDVI dependent and defined as follows:

NDV I − NDV Iv k

λ = vλ − ( vλ − sλ ) NDV Is − NDV Iv


where λ represents the spectral band, v is the vegetation emissivity (set to 0.99 in our case), and s is the soil emissivity, which may be adapted by user for more accuracy. NDV Iv is the maximum NDV I for full vegetation (set to 0.99) and NDV Is is the minimum NDV I for bare soil (set to 0.17).
The soil emissivity was not set because its value depends on the soil’s composition and roughness. In our case, the soil emissivity was set to 0.96 in agreement with the mean emissivity value extracted from the JPL library database [38] corresponding to the soil composition of our selected sites. The k parameter was fixed arbitrary to 2, but note that [39] proposed explore values between 2 and 3. This parameter can be adjusted in the future version. Another difficulty in determining soil emissivity is the non-homogeneity of the soil at the pixel scale. However, for acquisitions at a resolution of 100 m, Zhang et al. [40] found that the pixels can be considered homogeneous.

Table 1. Landsat thermal bands specifications.

L5 L7 L8 L8

Wavelength (µm)
10.40–12.50 10.40–12.50 10.60–11.19 11.50–12.51

Native Resolution
120 m 60 m 100 m 100 m

B6 B6 B10 B11

2.5. Ground Measurements
The ground measurements used for the local validation of the LST estimates were retrieved from two meteorological network stations one local near Kairouan in central Tunisia (35◦18 17.14 N, 9◦54 56.62 E) the other located near Toulouse in southwestern France and composed of two plots (43◦54 97 N, 01◦10 61 E and 43◦49 65 N, 01◦23 79 E) identify hereafter as respectively Aurade and Lamasquere. A localization map for both networks is given in Figure 1. The Tunisian site is representative of a semi-arid climate, with high surface temperatures during summer. The French site is representative of a continental climate with moderate summer surface temperatures.
The Tunisian site is part of the SICMED program (http://www.sicmed.net/). The measurements considered in the present study were acquired on non-irrigated olive tree fields composed of trees spaced at regular intervals of 20 m and with crown diameters of about 5 m. The ground was a sandy bare soil. The French site was composed of two cultivated plots located at a distance of 12 km from each other and with different soils and crop rotations. Aurade is a crop rotation of wheat/rapeseed and sunflower each 5 years; Lamasquere is a rotation of wheat and silage maize. The measurement methodology for both sites follows ICOS (integrated carbon observation system) recommendations. A complete description of the agricultural practices and field instrumentation for the two French

Remote Sens. 2016, 8, 696

8 of 24

plots are available in [41,42]. For both zones of interest, two kinds of micro-meteorological data were measured and processed to estimate local LST. First, directional remote temperatures were obtained with an IR-120 thermo-radiometer (Campbell Scientific) located at 2 m height and nadir oriented with a spectral window of 8–14 µm; Second, incoming and outgoing long-wave radiation was measured with a hemispherical CNR1-Pyranometer (Kipp Zonen) with a spectral window of 5–50 µm. Surface brightness temperatures were measured directly with the thermo-radiometer and were then compared with Landsat LST estimations after emissivity correction. For the hemispherical pyranometer measurements, the surface brightness temperature was first estimated with the Stefan-Boltzmann law:


 LW↑  4

Tbλ −λ = 

λ1 −λ2




where Tbλ1−λ2 is the surface brightness temperature, and LWλ↑1−λ2 is the outgoing long-wave radiation for the spectral window λ1 − λ2 equal to 5–50 µm.

Figure 1. Localization map of the both network (green triangles), for France (left) and for Tunisia (right). Background are true-color images of L7 (left) and L8 (right) acquired respectively on 29 May 2003 and 29 March 2013. Images are scaled and oriented. The black boxes on images represent the ASTER image footprint used on each location to validate tool (see Section 4.2).

Then, the surface brightness temperature was converted into surface temperature using the approach proposed in [34]:


+ (1 −

) Tb

(1 − )LW↓




λ 4 fλ(Tbλ)σTbλ3

where λ represents the (λ1 − λ2) band spectral window, σ is the Stefan-Boltzmann constant, LWλ↓ the incoming long-wave radiation, fλ(Tbλ) is a function factor corresponding to the fraction of energy emitted in the λ spectral window by a black body at temperature T relative to the energy emitted
over the full spectrum as defined in [43] and approximate to unity in our configuration, and is the
soil emissivity.

Remote Sens. 2016, 8, 696

9 of 24

The data used in this study are acquired between 2011 and 2013 for France and between 2012 and 2014 for Tunisia. The available data allowed us to compare more than three years of local thermal measurements with our estimates obtained from Landsat cloud-free images.
2.6. Landsat Data
Landsat data are the main inputs for LANDARTs. The present study used data retrieved by three missions: Landsat 5 TM (L5), Landsat 7 ETM+ (L7), and Landsat 8 OLI-TIRS (L8). Both L5 and L7 provide a single-band in the thermal domain (band 6). Landsat 8 provides two thermal bands (10 and 11). However, the USGS recommends against the use of band 11 because considerable uncertainty persists regarding the quality of the acquired data [11,12]. The thermal band characteristics for the three sensors are summarized in Table 1. Image data for the studied area and for the period of interest were downloaded from the USGS website (http://glovis.usgs.gov/). Note that all Landsat thermal data provided by the USGS were re-sampled to a resolution of 30 m (Table 1). In LANDARTs data are read as digital numbers and converted into sensor radiances (Lsen) in accordance with the calibration coefficients given in [24] for Landsat 5 and 7 and in [25] for Landsat 8.
The validation of our atmospheric correction tool considered two main areas corresponding to the ground dataset measurements (Section 2.5). The first area was located in south-western France (near Toulouse) at the intersection of two tiles: 199/030 for the western portion, and 198/030 for the eastern portion, in accordance with USGS WRS-2 nomenclature. The second area was located in Tunisia (near Kairouan) on tile 191/035. The downloaded Landsat images for France were acquired between 1 January 2011 and 31 December 2013. For the Tunisia site the images were acquired between 5 April 2012 and 31 December 2014.
2.7. ASTER Data
ASTER thermal channels constitute a valuable dataset for our investigation because both channel 13 and 14 (Table 2) are similar to the single TIR channel of L5 and L7 (Table 1). The L8 channel 10 also covers a large part of ASTER channel 13. It is therefore possible to compare Landsat and ASTER products for certain ranges of wavelengths. The ASTER database provides several data products. However, our study only had access to one higher-level product: the surface kinetic temperature called AST_08, which is the result of the TES algorithm applied to the five thermal channels at a resolution of 90 m [44].

Table 2. ASTER thermal band specifications.

Band No.
10 11 12 13 14

Wavelength (µm)
8.125–8.475 8.475–8.825 8.925–9.275 10.25–10.95 10.95–11.65

Acquisition Resolution
90 m 90 m 90 m 90 m 90 m

All cloudy images were discarded to prevent problems with automatic cloud detection. A comparison between the ASTER dataset covering southwestern France for the period from 2000–2006 and the USGS database revealed an acquisition time overlap of only a few dates. All dates after 2003 were discarded, to avoid large bands of missing data on the L7 images due to SLC failure and unfortunately corresponding to the ASTER position in the Landsat images. The ASTER dataset for Tunisia consisted of only two images, one of which was cloudy. Incidentally, an L8 image was available for the only cloud-free ASTER image date. The data selected and used in the present study are presented in Section 4.2. It was our aim to implement identical procedures to evaluate the performance of the proposed tool for both sites; we consequently considered only the surface kinetic temperature product (the only product available for the Tunisian site).

Remote Sens. 2016, 8, 696

10 of 24

3.1. Main Algorithm
LANDARTs is a Python tool that works (currently) only with Landsat archives in the new metadata file format (data acquired after 2012) that can be downloaded from the USGS website. Although the tool accepts a custom installation of MODTRAN, we tested the following three versions only: 4.1 on Unix, 4.3 on Windows, and 5.2 on Unix and Windows operating systems. A flowchart for the proposed tool is presented in Figure 2. The implementation process consists of three steps to produce LST maps at a resolution of 30 m using Landsat 5, 7 or 8 images. Step 1 consists of calculating the three atmospheric correction parameters. The details of the parameter calculations are described below:
1. Acquisition dates and times are read automatically according to the metadata file provided by the USGS data. The four corners of the Landsat footprint are also extracted from GeoTiff metadata.
2. The information obtained in (1) are used to create four spatial queries to the ECMWF server: for each two time samples bounding the time input, two queries are done, one for the surface and one for the pressure level dataset. The query zone is defined to be larger than the Landsat footprint to allow interpolations at the edges. Data are in the uniform latitude/longitude grids in one of the resolution proposed by ERA-Interim and selected by user.
3. As explained in the ECMWF Section 2.3, the two time profiles constructed with surface and pressure levels for the height from the Earth’s surface to an elevation of 100 km are linearly interpolated between the two times samples to give acquisition time [18]. This interpolation calculus is performed for each point of the ECMWF grid to obtain a 3D matrix (2D spatial and 1D vertical) of an equivalent atmospheric variables profiles at the acquisition time.
4. In accordance with the MODTRAN documentation [19], the 10 required input cards are written using the atmospheric profiles obtained in (3). The cards are used to create two tape5 files as MODTRAN inputs. This process is repeated for each point of the ECMWF grid.
5. MODTRAN run on each point of the grid defined by atmospheric profiles. When done, the MODTRAN predicted per-wavelength transmission, upwelling radiance and downwelling radiance are integrated over the instrument’s spectral response. According to the metadata, the 2D matrices of the three parameters are finally interpolated using the gdalwarp utility from the open source Geospatial Data Abstraction Library (GDAL, http://www.gdal.org/gdalwarp.html) to create three GeoTiff images oversampled at the resolution of 30 m. These three final parameter’s images can then be used for pixel-by-pixel processing.
Step 2 consists in converting the digital numbers of multispectral data into reflectance using the metadata and the USGS formulas. Next, Red and NIR bands are used to calculate the NDVI. The NDVI map were then converted into an emissivity map at the data resolution of 30 m according to Equation (8).
Finally, Step 3 handles the thermal data: in accordance with the USGS user guide, the digital numbers of the thermal band are converted into top of atmosphere radiance using the metadata file. The emissivity estimates derived from the NDVI in Step 2 are then used to calculate the ground radiance, and LST are converted by inverting radiative transfer Equation (3) with spatialized atmospheric parameters and simplified Planck’s law (4). The LANDARTs tool is available free of charge by contacting the corresponding author. However, note that the user needs to have his or her own official instance of MODTRAN.