Assessment of Methods for Land Surface Temperature Retrieval from Landsat-5 TM Images Applicable to Multiscale Tree-Grass Ecosystem Modeling

Land Surface Temperature (LST) is one of the key inputs for Soil-Vegetation-Atmosphere transfer modeling in terrestrial ecosystems. In the frame of BIOSPEC (Linking spectral information at different spatial scales with biophysical parameters of Mediterranean vegetation in the context of global change) and FLUXPEC (Monitoring changes in water and carbon fluxes from remote and proximal sensing in Mediterranean ―dehesa‖ ecosystem) projects LST retrieved from Landsat data is required to integrate ground-based observations of energy, water, and carbon fluxes with multi-scale remotely-sensed data and assess water and carbon balance in ecologically fragile heterogeneous ecosystem of Mediterranean wooded grassland (dehesa). Thus, three methods based on the Radiative Transfer Equation were used to extract LST from a series of 2009–2011 Landsat-5 TM images to assess the applicability for temperature input OPEN ACCESS Remote Sens. 2014, 6 4346 generation to a Landsat-MODIS LST integration. When compared to surface temperatures simulated using MODerate resolution atmospheric TRANsmission 5 (MODTRAN 5) with atmospheric profiles inputs (LSTref), values from Single-Channel (SC) algorithm are the closest (root-mean-square deviation (RMSD) = 0.50 °C); procedure based on the online Radiative Transfer Equation Atmospheric Correction Parameters Calculator (RTE-ACPC) shows RMSD = 0.85 °C; Mono-Window algorithm (MW) presents the highest RMSD (2.34 °C) with systematical LST underestimation (bias = 1.81 °C). Differences between Landsat-retrieved LST and MODIS LST are in the range of 2 to 4 °C and can be explained mainly by differences in observation geometry, emissivity, and time mismatch between Landsat and MODIS overpasses. There is a seasonal bias in Landsat-MODIS LST differences due to greater variations in surface emissivity and thermal contrasts between landcover components.


Introduction
Land surface temperature (LST) is a state variable that plays a crucial role in many land surface processes [1]. LST is related to the transport of heat between the land surface and the atmospheric boundary layer [1][2][3], and makes possible estimation of sensible heat flux [4] and latent heat flux, or evapotranspiration [5,6]. It is a necessary input for ecosystem modeling [7], which can be performed at local [4], regional, and global scales. While local modeling relies heavily on field data, remote sensing has become the main source for LST estimation at the regional and global scales [8].
Radiance measured at a sensor can be transformed into LST by inverting the Radiative Transfer Equation (RTE) applied to a particular thermal IR band or wavelength: L sensor = τεL Ts + L u + τ(1 − ε)L d (1) where L sensor is the radiance registered by the sensor, also referred to as top of atmosphere radiance, L Ts is the blackbody radiance related to the surface temperature by Planck's law and Ts is the LST, L u and L d are the upwelling and downwelling atmospheric radiances, respectively (all the radiances in W•sr −1 •m −2 •μm −1 ), τ is the atmospheric transmissivity and ε is the land surface emissivity. In the case of dealing with a waveband, all these parameters are integrated according to the spectral response function of this band. The signal coming from the target to the sensor is modified as it passes through the atmosphere, which both emits and absorbs thermal radiation. The latter effect is mainly caused by the presence of water vapor. When atmospheric conditions are known, emission and absorption of radiation in the atmosphere can be quantified and corrected using one of the radiative transfer computer codes, e.g., MODerate resolution atmospheric TRANsmission (MODTRAN) [9]. Atmospheric conditions are typically assessed using in situ atmospheric profile data, which are often not available for the place and time the image was acquired, although on-line atmospheric databases [10,11] or estimations based on empirical models [12] can be used.
At present, there are several satellites providing global data from the thermal region of the spectrum at different scales. Among them are MODIS [13] and Spinning Enhanced Visible and Infrared Imager (SEVIRI) [14] characterized by low spatial and high temporal resolutions, for which LST products are available on a regular basis. At the medium spatial scale Landsat has provided global brightness temperatures since 1984, with Landsat 8 launched at the beginning of 2013 giving continuity to the data record [15]. The assessment of methods for LST estimation from a unique thermal band gains additional importance if we consider problems with data from one of the Landsat 8 thermal bands (band 11) and National Aeronautics and Space Administration (NASA) suggestion not to use band 11 for surface temperature retrieval [16]. The recently published reviews [8,17] mention several single-channel methods based on approximations from the RTE, which can be applied for LST retrieval from Landsat-5 unique thermal band [18][19][20][21]. These methods perform atmospheric correction based on water vapor content [19,20] or both water vapor and near-surface air temperature [18,21]. Apart from the atmospheric correction parameters, the surface emissivity (defined as the ratio between the target emitting capacity and that of a blackbody at the same temperature) is also required. A review of methods for surface emissivity estimation from satellite data is available in Li et al. [22]. Because of the high level of correlation between NDVI and surface emissivity, many methods proposed for estimating emissivity are based on this vegetation index [23][24][25][26][27].
One of the research fields with a great demand of LST data at a local scale is carbon and water fluxes modeling in terrestrial ecosystems. BIOSPEC (Linking spectral information at different spatial scales with biophysical parameters of Mediterranean vegetation in the context of global change) [28] and FLUXPEC (Monitoring changes in water and carbon fluxes from remote and proximal sensing in Mediterranean -dehesa‖ ecosystem) [29] projects carry out the analysis of these processes using information from ground-based measurements of fluxes and vegetation biophysical parameters, and their modeling throughout the integration of spectral data from remote sensors having different spatial, spectral and temporal resolutions (Landsat and MODIS) following the attempts of other scientific teams [30,31]. Landsat can provide LST at a spatial detail much higher than MODIS, but only once in 16 days compared to daily images acquisition by MODIS. Thus, integration of the data from these two satellites would be highly beneficial given the spatial resolution of the former and the temporal resolution of the latter. However, the challenges and persisting uncertainties related to the use of Landsat for LST estimation [32], especially in heterogeneous environments, make it necessary to evaluate the methods and atmospheric information sources looking for those more similar to MODIS. Although there are a number of studies comparing methods for LST retrieval from one thermal band [18,19,33,34], the evaluation is usually based on data from homogeneous environments. On the other hand, this study presents an assessment of the single-channel methods in heterogeneous environments common for most of the land surface.
Our main interest in this study is to compare the performance of the most common methods for LST retrieval from Landsat-5 TM images of the dehesa tree-grass ecosystem [8] and analyze the relationship between LST estimated from Landsat and LST from MODIS product (MOD11_L2), for the use in Landsat-MODIS LST fusion algorithm development to study energy and water exchange between the dehesa landcover and the atmosphere. Three procedures are applied for LST retrieval from a sequence of 13 images of Central Spain, acquired from 2009 to 2011: (1) RTE inversion with atmospheric correction parameters calculated by on-line ACPC tool [10], which is referred to as Radiative Transfer Equation Atmospheric Correction Parameters Calculator (RTE-ACPC) from here on and two methods, which are approximations of the RTE with minimum parameters: (2) single-channel (SC) method by Jimé nez-Muñoz and Sobrino [20], updated in 2009 [19], and (3) mono-window MW method by Qin et al. [21]. The results are compared with LSTs simulated by Radiative Transfer Code MODTRAN 5. We also assess and analyze the relationship existing between Landsat LSTs and those from MODIS LST product (MOD11_L2). In situ grass temperature measurements available for some of the images complete the set of reference data.

Study Area
The study area shown in Figure 1 is located in a dehesa ecosystem near the Las Majadas del Tietar FLUXNET site (geographic coordinates: Lat 39°56′26′′N, Long 5°46′29′′W), which is operated by the Mediterranean Center for Environmental Studies (CEAM). FLUXNET is a network of micrometeorological observation sites established to perform continuous measurement of exchange fluxes in the soil-vegetation-atmosphere system [35].
The dehesa is an open savanna with an integrated agroforestry ecosystem, and has a complex vegetation structure typical of Mediterranean areas. The study site is flat, and is covered by grass (75% of the area) and holm oak trees Quercus ilex ssp. rotundifolia (25% of the area). The zone climate (Csa according to Köppen classification) is characterized by an annual average temperature of 16 °C and approximately 550 mm precipitation, and has a four-month hot dry period from June to September [36].

Landsat-5 TM Images
Landsat-5 TM provides images with six bands in the optical region, and a thermal band with a bandwidth of 10.4-12.5 μm. The LST was retrieved from 13 Landsat-5 TM (path 202, row 32) clear sky images pre-processed by the NLAPS (National Land Archive Production System-USGS) and downloaded from [37] ( Table 1). The images over the study area were acquired at approximately 10:50 a.m. GMT from 2009 to 2011.

MODIS LST Images
The MODIS Terra LST MOD11_L2 product with a 1-km pixel spatial resolution was used for comparison. MOD11_L2 constitutes an output of the split window algorithm [38] applied to MODIS bands 31 (10.780-11.280 µm) and 32 (11.770-12.270 µm). The time difference between Landsat and MODIS passes over the study area is about 20 min (Table 1): MODIS images are acquired approximately 20 min earlier. FLUXNET tower data corresponding to the same dates show an average air temperature increase of about 0.5 °C for the same time period, while in situ grass surface temperature measurements available for three summer dates in 2011 (Table 2) demonstrate an average increase of 1.5 °C . Following the procedure applied by other researchers [39,40] to account for different spatial resolution of the sensors, MODIS temperature value corresponding to a pixel centered in the study area was compared with the mean value of the Landsat-5 TM pixels within that MODIS pixel. Moreover, to minimize the effects of the differences in the observation geometry only the images with the best quality MODIS pixel of the study area (MODIS product quality flag 0) were used for the comparison. According to the MOD11_L2 product description quality flag 0 is assigned to the cloud-free pixels with LST error less than 1 °C and the emissivity errors in channels 31 and 32 involved in LST estimation less than 0.01.

Atmospheric Correction Parameters Sources
We obtained and compared data on the atmospheric water vapor content from three online sources: Aerosol Robotic Network (AERONET) database, National Center for Environmental Prediction (NCEP) Reanalysis (hereafter called REANALYSIS) database and from MODIS MOD05 product. AERONET is part of the NOAA Observing System Architecture, which includes more than 500 sites distributed worldwide. Precipitable water content values (g•cm −2 ) were downloaded from an online database [41] for Cá ceres; the observation site is located approximately 50 km from the study area. The National Center for Environmental Prediction (NCEP) and the National Center of Atmospheric Research Reanalysis Project (NCAR) maintain a free access online database of gridded and continuously updated meteorological data at 2.5° × 2.5° spatial and 6 h temporal resolution extending back to 1948 [42]. Precipitable water values (kg•m −2 ) for 2009-2011 were downloaded from [43]. The noon values, approximately 1 h later than the Landsat-5 TM overpass, were extracted for the study area location and used in the water vapor sources comparison. Atmospheric profiles containing information on vertical distribution of pressure, geopotential height, temperature and relative humidity for simulation of the reference LSTs were generated by ACPC tool based on the interpolation of the NCEP profiles resampled to 1° × 1° spatial resolution [11]. Interpolated profiles were completed with the data from the standard atmospheres for the altitude range from 30 km to 100 km and user-supplied information for the lowest level, resulting in the 31 levels in each profile. Precipitable water from MODIS MOD05 product at 1-km spatial resolution close in time to Landsat overpass was obtained from MODIS web archive [44].
FLUXNET tower was used as the source of ACPC tool meteorological inputs. Due to the limited extension of the study site, meteorological data provided by the tower were considered characteristic for all the analyzed area.

In Situ Grass Temperature Measurements
To put the obtained results in site context and take into account the difference in LST between the overpass times of Landsat and MODIS on board of Terra (from Latin -land‖) satellite, we used the data on grass temperature obtained from an infrared sensor Campbell IR120 installed on a tower at a height of 8 m ( Table 2). The sensor registers data every 10 min with an accuracy of ±0.2 °C . The data are available for a part of 2011 beginning 3 March 2011. The device offers a non-contact means of measuring the surface temperature of an object by sensing the infrared radiation in the wavelength range of 8 to 14 μm in the field of view of 20°. The in situ LSTs coincident with the Landsat image acquisition (10:50 a.m. GMT) were only used to assess the significance of time mismatch between Landsat and MODIS TERRA overpasses because the data are available only for one of the landcover components (grass) and for less than 25% of the images.

Land Surface Temperature (LST) Estimation
Prior to LST retrieval optical bands of Landsat images used in emissivity estimation were corrected for atmospheric effects using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) algorithm implemented in the ENVI (software package, a geospatial imagery analysis and processing application marketed by Exelis Visual Information Solutions) [45]. The LST was retrieved from the thermal band; the digital numbers were first converted into radiance using the header files parameters and then to the at-sensor brightness temperature, which was then transformed to LST. Three procedures used to transform the at-sensor brightness temperature into LST are: (1) RTE inversion using atmospheric correction parameters from on-line ACPC tool [10] available at [46]; and two algorithms based on the approximations of RTE: (2) single-channel SC method [19,20]; and (3) mono-window MW method [21]. The most recent SC modification [18] is highly sensitive to water vapor changes and was not considered, because in situ measurements of water vapor content were not available. Since LST estimation methods require clear sky, only cloud-free images were used for processing.

Radiative Transfer Equation (RTE)
As mentioned in Section 1, LST can be obtained from RTE (Equation (1)) and Planck's law inversion once parameters for the atmospheric corrections (L u , L d and τ) are estimated and the surface emissivity is known. The first tested procedure used the atmospheric correction parameters from the Atmospheric Correction Parameter Calculator (ACPC). It is an on-line tool developed for atmospheric correction of the Landsat 5 and 7 thermal data using MODTRAN 4 radiative transfer code [10,11]. The tool receives as input user-provided information on geographical coordinates, site elevation, date and time of the image acquisition and calculates site-specific atmospheric transmission, upwelling, and downwelling atmospheric radiances to be used in LST estimation through RTE inversion. Henceforth, the LST values obtained in the study by this procedure are referred to as RTE-ACPC. NCEP atmospheric databases are used to interpolate the profile for the specified place, date, and time; the profiles resulting from time interpolation can be provided for the closest lat/long grid corner or interpolated for the user-specified location. The latter option was used in this study. The tool processes data corresponding to one set of conditions (one Landsat image) at a time; the results are forwarded to the user's e-mail address. The set of parameters generated by the tool for the images analyzed in this study is presented in Table 3. According to developers, the tool provides parameters allowing LST estimation through RTE (Equation (1)) inversion within ±2 °C [11]. In the MW algorithm [21] the LST is determined through decomposition of Planck's radiance function using a Taylor's expansion and calculation of two empirical coefficients a and b. Three a priori known parameters are required for the algorithm: transmissivity (τ)/water vapor content, effective mean atmospheric temperature (T a ) and emissivity (ε). All the temperatures are in K. LST (Ts) is calculated from the equation (2): where a = −67.355351 and b = 0.458606 are constants, Tsensor is the at-sensor brightness temperature, C and D are calculated using Equation (2a,2b) respectively: The suggested method for calculation of T a is based on the relationship between T a and the vertical water vapor distribution in the atmosphere [47]. Simulations performed using LOW resolution TRANsmission 7 (LOWTRAN 7) [21] indicate that, while water vapor content differs significantly depending on the atmospheric conditions, the distribution of the ratio of water vapor content at a particular altitude to the total is very similar for all atmospheric profiles. This enabled formulation of the Equation (3a-3c) for calculation of T a from the total water vapor content and the near surface local air temperature (T 0 ), according to the atmospheric conditions [21]: The most important parameter of the algorithm τ is estimated using the expressions obtained from simulations using LOWTRAN 7 [21] for two air temperature profiles: Equation (4a,4b) for high (35 °C) and Equation (4c,4d) for low (18 °C) [21]: The algorithm performs well for atmospheric conditions where the water vapor content is 0.5-2.5 g•cm −2 [18,19,21].

Single-Channel (SC) Method
SC method [19,20] is also an approximation of RTE and requires only atmospheric water vapor content for atmospheric correction. In this method LST is obtained from the following Equation (5) where w is total atmospheric water vapor content in g· cm −2 . Similar to the MW, the optimal performance of the SC algorithm is observed for the atmospheres with water vapor content in the range of 0.5-2.5 g•cm −2 [18,19,21].

Reference Land Surface Temperature (LST)
Because of the incompleteness of the in situ data, LSTs simulated by the latest version of the radiative transfer code MODTRAN 5 are used as a reference set. As suggested in previous studies [8,17,48,49], LSTs simulated using radiative transfer code can be an alternative for validation when field measurements at a required spatial scale are not available. The method was earlier applied for Landsat [34] and MODIS [48,49] LST assessment. Among the most important improvements in MODTRAN 5 compared to MODTRAN 4 is the incorporation of band model parameters based on HITRAN2008, with 2009 updates [9]. MODTRAN 5 performs calculations based on the information about observation geometry and atmospheric profiles at the moment of observation. The best results are achieved when data come from in situ radiosoundings synchronized in time with image acquisition. Unfortunately, they were not available in this study. When discussing the difficulty of obtaining local radiosounding data, multiple studies [17,50,51] suggest the use of the atmospheric profiles from the reanalysis products as a viable solution. Thus, we use NCEP atmospheric profiles interpolated for the exact location and time of Landsat overpass, the choice validated by previous research [17,50,51]. The NCEP atmospheric profiles interpolated for the study area and conditions by ACPC tool are complemented with on-site meteorological data for the lowest atmospheric layer, which together with the newer MODTRAN version (5 vs. 4) marks the difference with the RTE-ACPC procedure. To simulate the reference LSTs, the profiles are inserted into MODTRAN input file. Then the first MODTRAN run is performed with 0% surface albedo; atmospheric transmissivity (τ) and upwelling radiance (L u ) are extracted from the MODTRAN output files and integrated over the Landsat-5 TM thermal band using the sensor filter function. To calculate downwelling radiance (L d ) MODTRAN 5 is run for the second time with 100% surface albedo. Next, the obtained atmospheric correction parameters τ, L u and L d together with previously estimated emissivity ε are substituted into RTE (Equation (1)) to calculate the radiance from the target (L Ts ). The final step consists in transformation of the calculated target radiance into LST (LST ref ) by inversion of the Planck's law.

Emissivity Estimation
Most of the emissivity retrieval methods from remotely sensed data, such as TES [52] or TISI [53] cannot be used with Landsat images because there is only one thermal band. The possible solution is to apply one of the methods based on the normalized difference vegetation index (NDVI) [22]. Among the advantages of these methods is that they rely on the information from the image used for the LST retrieval [22]. The NDVI thresholds method (NDVI THM ) [25,54] based on the findings of Valor and Caselles [26] was applied to estimate surface emissivity in this study. The emissivity of the pixel is determined based on its NDVI. Different functions are applied to calculate emissivity depending on the NDVI range (Table 4).
In case of the mixed pixels category the NDVI values (thresholds) selection is based on an analysis of the images histograms. The soil emissivity εs value of 0.984 is based on in situ field measurements using box method [24] with an estimated error of 0.003 [24], and is similar to the values reported by previous research [34,55]. The vegetation emissivity v  is assigned the value of 0.990 [34]; ε d = 0.01 is the term accounting for surface roughness different from zero for heterogeneous covers [3]; The validation of NDVI THM method performed by Sobrino et al. [34] gets the error of less than 0.01, which in terms of LST would mean the error below 0.5 °C [26]. Of the three dehesa landcover components, soil emissivities show the greatest variation in the thermal region of the spectrum [24,34]. As the soil emissivity measured in situ is high in present study, the related error should be smaller.

Results and Discussion
We present and discuss below the results of LST estimation in heterogeneous Mediterranean tree-grass (dehesa) ecosystem with the RTE-ACPC, MW and SC procedures described in Section 3.1. Emissivity ε is calculated using the NDVI Thresholds method presented in Section 3.2. Section 4.1 compares three sources of the atmospheric water vapor (w) and explains the choice of the NCEP REANALYSIS for this study. Section 4.2 analyses the differences between the LST ref and LST generated by the tested procedures. Next, Section 4.3 discusses the relationship between Landsat LST and MODIS LST product. Both LST comparisons (LST ref and MODIS) include the use of the in situ values of grass temperature measured in 2011 to assess the implications of time mismatch on the LST differences.

Atmospheric Water Vapor Content
Atmospheric conditions on the images acquisition dates are shown in Table 5. The registered mean w values were relatively low (1.292 g•cm −2 , 1.515 g•cm −2 , and 1.600 g•cm −2 for REANALYSIS, AERONET, and MODIS, respectively), and the maximum values were close to 2.5 g•cm −2 . Therefore, the data were considered adequate as inputs to the MW and SC methods. The average difference between w sources was around 0.3 g•cm −2 .
A detailed case-by-case analysis revealed important differences among databases on some dates. For example, the difference between MODIS and other sources was greater than 0.7 g•cm −2 for 4 August 2011, while AERONET exceeded w values from REANALYSIS in more than half a gram per square centimeter on 30 June 2010, and 17 October 2009. Although a clear pattern of differences among the data sources was not observed, the REANALYSIS water vapor values were lower than those of the other two databases; only once the w value from this source was marginally greater than the value from MODIS (11 April 2010) and in two cases the w levels were greater than those of the AERONET database (1 August 2010, and 4 August 2011). The comparison of three different atmospheric water vapor (w) sources did not reveal statistically significant differences between them (F-Test = 1.16; p-value > 0.05). Hence, the REANALYSIS w values were used in atmospheric correction since this database is the result of modeling which assimilates data from multiple sources and is continuously updated. We did not use the MODIS product as a w source, because one of the objectives of the study is the comparison of the Landsat-retrieved LSTs with those from MODIS LST product, which employs MOD05 w values in the algorithm. There were no statistically significant differences between the values obtained using tested procedures (F-Test = 0.111; p-value > 0.05) and the degree of correlation between the values obtained by different methods is very high (R 2 > 0.986). It is not strange considering that all the four algorithms are based on successive versions of the same radiative transfer code: LOWTRAN 7 (Mono-Window (MW)), MODTRAN 4 (RTE-ACPC and SC) and MODTRAN 5 (LST ref ), developed in 1988 [58], 1999 [59] and 2011 [9] respectively. Moreover, all of them employ the fewest (although different) possible number of parameters for atmospheric correction (w for SC; T a and w for MW; profiles of RH, T a and atmospheric pressure for RTE-ACPC and LST ref ) and the same emissivity.  (Table 7) [19], which show that MODTRAN 4 generates greater w values (around 1 g· cm −2 for high w values) resulting in higher LSTs. At the same time, the SC and RTE-ACPC (methods based on MODTRAN4) are closer to the in situ data: (averages of 5.46, 3.19, and 3.31 °C for LST in_situ -LST MW , LST in_situ -LST SC , and LST in_situ -LST RTE-ACPC , respectively), although this comparison is not fully accurate since LST in_situ corresponds only to grass component of the landcover. Even though MW systematically underestimates LST, the size of the differences varies from 0.11 to 4.66 °C depending on the date (Table 8); the range of variations for SC and RTE-ACPC is much smaller (below 1 °C and 3 °C for SC and RTE-ACPC, respectively). Considering that both procedures use the same emissivity, explanation of the anomalies lies in different sensitivity of the algorithms to atmospheric variables. Good correlation of the differences between LST ref and MW with w and air temperature (R = 0.8) can be appreciated in Figures 2 and 3; high atmospheric water vapor concentration and high temperatures in summer time explaining the biggest LST deviations. The same graphics reveal that there is no relationship between atmospheric parameters and the differences between SC and LST ref (R < 0.2). Bigger errors in hot and wet conditions have already been detected in other studies [19,50]. Modeling [60] shows that a typical w error of 10% [61] may lead to LST error of 0.4 K and 0.2 K for SC and MW algorithms respectively for summer atmosphere [60]. For MW it is also necessary to consider the 0.2 °C error due to the air temperature [21]. Because in MW algorithm coefficients are developed only for two air temperature values and a reduced number of standard atmospheres, the algorithm fails to represent real atmospheric conditions in the study area correctly, especially on in summer. However, SC incorporates atmospheric functions based on extensive atmospheric profile databases allowing more precise representation of atmospheric conditions over the study site at the moment of satellite pass [19,50].
Based on statistical analysis we can conclude that SC and RTE-ACPC procedures are capable of retrieving LSTs in the study area of Mediterranean tree-grass ecosystem with an error below 1 °C, which is similar to the results of the previous studies conducted in the homogeneous areas [34,62]. Thus, Sobrino et al. [34] compared LSTs from MW and SC methods applied to Landsat images with LSTs simulated using radiative transfer code and in situ emissivity in agricultural area obtaining the errors of around 0.9 °C for SC and around 2 °C for MW procedures; similar errors were reported by Copertino et al. [33] who applied the same methods for estimating LST over different landcover types in Southern Italy, in this case retrieved LSTs were compared to the soil temperatures. Limin et al. [60] compared LST estimated from HJ-1B satellite by MW and SC with MODTRAN 4 simulations of LST registering errors below 1 °C in summer for nadiral view of the sensor.

Landsat LST vs. MODIS Land Surface Temperature (LST)
We now present the comparison of Landsat LSTs and LSTs from MODIS LST product. Before the comparison some adjustment was performed to account for differences in data format and spatial resolution between Landsat and MODIS. MODIS LST images (MOD11_L2 product) were reprojected to match spatial reference of Landsat. Since the study site is in the middle of the much more extensive tree-grass ecosystem area with similar LST variability at the MODIS scale, the average LST value of the Landsat pixels inside the MODIS pixel covering the center of the study area is calculated for each date and method and is used for the comparison.
The results of the comparison with MODIS product LST and the intercomparison of the LST values retrieved by the tested methods (Table 7) show that SC and RTE-ACPC are more similar to each other than to the LSTs from MODIS product (RMSD of 4.16 and 4.29 °C for RTE-ACPC and SC respectively). On the contrary, the MW-estimated LST values are much closer to MODIS LSTs (RMSD of 2.27 °C ).
Compared to Landsat-estimated values MODIS product underestimates LST, the bias is 1.5 °C for MW and 3.5 °C for SC procedures. This is in agreement with the results reported in previous studies [40,63,64], which mention that LST values from MODIS product are lower than those obtained from other sensors or in situ measurements. Thus, Trigo et al. [65] observed a negative bias of 2.6 °C in MODIS LST compared with ground values, especially at night. The underestimation also occurs when comparing MODIS with other sensors, such as SEVIRI [65] and AATSR [40]. In case of AATSR sensor, which is the most accurate infrared radiometer currently being flown in space according to [40], the biases of −0.5 and −1.2 °C were observed both during day and night respectively. So, it is evident that there is a problem related to spatial scale differences, which makes complicated the comparison of satellite and in situ data [17,48,49,66], although the differences are also affected by other factors. One of the most important is the impact of the observation angles on the measurements: while Landsat angles of observation are almost nadiral, MODIS views the study area at an angle of 60°, i.e., the sensor observes the surface from the west, detecting higher fraction of shadow and vegetation surfaces considerably decreasing LST. Previous studies show that the differences in the LST measured in nadir and off-nadir observations can be as large as 5 °C depending on the angle and cover type [17,67].
On the other hand, the 21 min time mismatch in the study area overpass between the sensors (ranging from 7 min to 38 min, see Table 1) also operates in the same direction. The analysis of the time differences between Landsat and MODIS is performed using data from thermal sensor installed in the study area. The average temperature increment between 10:30 and 10:50 GMT for the three dates in 2011, all of them in summer, is around 1.5 °C (Table 2). These coincide with [66] who indicate that LST difference between the LSTs at the moments of Landsat and MODIS Terra overpasses can range from 0.8 to 2 K, depending on the vegetation cover. If this time mismatch and the corresponding surface temperature increase were taken into account the gap between Landsat and MODIS would be reduced. The LSTs were not adjusted because only grass temperatures are available, not so the temperatures of tree canopies and shadows. However, even though tree canopies cover only about 20% of the area, we would expect significant decrease of the LST due to their presence within the MODIS pixel, since some studies [68] indicate that the difference between the grass and tree canopy temperatures in summer can be around 6-15 °C depending on the species and time of the day.
Although MW apparently generates LST values, which are closer to those from MODIS, they may not be more accurate than LSTs estimated by other procedures. The similarity between MW and MODIS LSTs results from two trends acting in the same direction: one is the underestimation of the LST by MW algorithm due to the use of the older radiative transfer code version (LOWTRAN) and another is the underestimation of the LSTs by MODIS due to the differences in time and observation angles between MODIS and Landsat and implications of these differences on the emissivity.
When SC results (the closest to the LST ref ) are compared to MODIS LST, a seasonal bias is observed: the greatest variances (above 6 °C ) occur in summer (Table 9) and the lowest (0.00-0.38 °C ) in winter and autumn. This fact was already mentioned in other studies [51]. Trigo et al. [65] observed that greater LST dispersion in summer can be related to the great thermal contrasts between landcover components (bareground, grass, tree canopy) taking place during this season. Because of higher spatial resolution and higher variability in emissivity, Landsat is more sensitive to this dispersion. Greater thermal range of around 8 °C on summer dates can be appreciated in Figure 4 showing Landsat LST variability within MODIS pixel. We should also consider that MODIS surface emissivity estimation is based on landcover types from the map updated annually [69], while NDVI Thresholds emissivity algorithm used in this study is based on NDVI (see Section 2.4).  Another explanation for the magnitude of LST Landsat -LST MODIS is the greater spatial and temporal variability of emissivity values estimated from Landsat-5 TM NDVI. This wider range is caused by the higher spatial resolution of the Landsat-5 TM, different algorithms used for emissivity estimation for two sensors and differences in viewing angles between Landsat and MODIS Terra (almost nadiral for Landsat vs. around 60° viewing angles for MODIS Terra) resulting in greater sensitivity of Landsat to an increase in the soil component and greater temperature contrasts between areas with and without vegetation, characteristic to summer as a consequence of grass senescence.

Conclusions
The study demonstrates that LST of dehesa ecosystem can be estimated from Landsat-5 TM thermal band using SC and RTE-ACPC procedures with RMSDs lower than 1 °C and the RMSD of 2.3 °C using MW algorithm, with expected uncertainties in energy fluxes modeling of around 10-30 W•m 2 for SC and RTE-ACPC [17]. The differences with the reference LSTs (LST ref ) are due to the fact that the tested methods are based on the different versions of the radiative transfer code: LOWTRAN 7 for MW and MODTRAN 4 for SC and RTE-ACPC. Moreover, there is a seasonal bias in the MW results, as evident from the correlations between MW-LST ref and near-surface air temperature and atmospheric water vapor w (R = 0.8), explained by the worse fit of MW coefficients to real atmospheric conditions in the study area compared to other procedures. This dependence is not evident in the LSTs obtained by the SC and RTE-ACPC procedures.
On the other hand, the existing LST mismatch between Landsat and MODIS is due mainly to (1) the time differences in the satellites overpasses and (2) the differences in the viewing angles which make Landsat much more sensitive to changes in the proportion of different landcover components with high thermal contrasts (soil and vegetation) and decrease of emissivity, especially during hot summer months.
Considering the generally-accepted error at the level of 1-2 K [70,71], the three tested procedures (SC, RTE-ACPC, and MW) can be used for LST estimation from Landsat-5 TM thermal data. RMSDs obtained for SC and RTE-ACPC procedures are below 1 °C , with the best results for SC (RMSD = 0.5 °C). This algorithm, which does not require radiosounding data, is considered the most adequate for integration with LST from MODIS MOD11_L2 product. However, the between-sensors differences due to time mismatch and observation angles should be taken into account. It was not possible to estimate the precise magnitude of Landsat-MODIS LST differences due to the lack of information on the contribution of each of the landcover components to ensemble radiance from heterogeneous and non-isothermal pixel characteristic for dehesa ecosystem.
Group AIRE of Physics Department, University of Extremadura, for establishing and maintaining the Cá ceres AERONET site database used in this research.

Author Contributions
This research was conducted by Lidia Vlassova and Fernando Perez-Cabello. Lidia Vlassova performed data processing and modelling. Hector Nieto and Pilar Martí n contributed to interpretation of the results and supervision of the methods employed. David Riaño and Juan de la Riva contributed to the organization of the paper. All authors helped in editing and revision of the manuscript, and responding to reviewers comments.