Land Surface Temperature Retrieval from Landsat 5, 7, and 8 over Rural Areas: Assessment of Different Retrieval Algorithms and Emissivity Models and Toolbox Implementation

Abstract

Land Surface Temperature (LST) is an important parameter for many scientific disciplines since it affects the interaction between the land and the atmosphere. Many LST retrieval algorithms based on remotely sensed images have been introduced so far, where the Land Surface Emissivity (LSE) is one of the main factors affecting the accuracy of the LST estimation. The aim of this study is to evaluate the performance of LST retrieval methods using different LSE models and data of old and current Landsat missions. Mono Window Algorithm (MWA), Radiative Transfer Equation (RTE) method, Single Channel Algorithm (SCA) and Split Window Algorithm (SWA) were assessed as LST retrieval methods processing data of Landsat missions (Landsat 5, 7 and 8) over rural pixels. Considering the LSE models introduced in the literature, different Normalized Difference Vegetation Index (NDVI)-based LSE models were investigated in this study. Specifically, three LSE models were considered for the LST estimation from Landsat 5 Thematic Mapper (TM) and seven Enhanced Thematic Mapper Plus (ETM+), and six for Landsat 8. For the accurate evaluation of the estimated LST, in-situ LST data were obtained from the Surface Radiation Budget Network (SURFRAD) stations. In total, forty-five daytime Landsat images; fifteen images for each Landsat mission, acquired in the Spring-Summer-Autumn period in the mid-latitude region in the Northern Hemisphere were acquired over five SURFRAD rural sites. After determining the best LSE model for the study case, firstly, the LST retrieval accuracy was evaluated considering the sensor type: when using Landsat 5 TM, 7 ETM+, and 8 Operational Land Imager (OLI), and Thermal Infrared Sensor (TIRS) data separately, RTE, MWA, and MWA presented the best results, respectively. Then, the performance was evaluated independently of the sensor types. In this case, all LST methods provided satisfying results, with MWA having a slightly better accuracy with a Root Mean Square Error (RMSE) equals to 2.39 K and a lower bias error. In addition, the spatio-temporal and seasonal analyses indicated that RTE and SCA presented similar results regardless of the season, while MWA differed from RTE and SCA for all seasons, especially in summer. To efficiently perform this work, an ArcGIS toolbox, including all the methods and models analyzed here, was implemented and provided as a user facility for the LST retrieval from Landsat data.

1. Introduction

Remote sensing technology is an important source of Earth observation from different platforms and sensors, and it offers work on a large scale with cheap, accurate (depending on the research design), and faster results compared to the conventional methods. Thermal remote sensing is one of the branches of remote sensing that deals with the acquisition, processing, and interpretation of data acquired primarily in the Thermal Infrared (TIR) region of the Electromagnetic (EM) spectrum [1,2,3]. Thermal remote sensing captures the radiation emitted from the ground primarily to estimate the surface temperature. In addition to surface temperature, surface emissivity, soil moisture, and evapotranspiration are the other crucial biophysical parameters estimated from TIR observations. Since these parameters govern the land-atmosphere interactions and the energy fluxes, their accurate evaluation is required to understand the behavior of the Earth.

Land Surface Temperature (LST) represents the temperature of the Earth’s surface, and it is one of the key parameters that affect surface energy balance, regional climates, heat fluxes, and energy exchanges [4,5,6,7,8,9,10,11,12,13,14,15]. Many researchers have investigated the importance and effects of LST on various topics, including urban climate and Surface Heat Island (SHI) studies [16,17,18,19,20], evapotranspiration [21], forest fire monitoring [22], geological, and geothermal studies [23,24,25,26,27]. Besides, LST has been approved as one of the high-priority parameters for the International Geosphere and Biosphere Program (IGBP) [7,28]. LST can be estimated from radiance measurements by meteorological stations. However, this method does not generally allow a large scale monitoring since it is a point-based measurement [29,30]. Remotely sensed TIR data allow temporal and spatial LST analysis on a large scale, even globally [31].

Accurate LST retrieval from TIR data depends on atmospheric effects, sensor parameters, i.e., spectral range and viewing angle, and surface parameters such as emissivity and geometry [32,33,34,35,36,37,38]. Since emissivity and atmospheric effects are two fundamental factors to derive LST from thermal data, many researchers have proposed different approaches for LST retrieval considering these factors [39,40,41,42,43,44,45,46]. These algorithms are named considering the number of TIR bands used. For instance, single-channel or mono-window algorithms use one TIR band. However, split window or multi-channel methods include more than one TIR band.

Accuracy assessment of space-based LST retrievals is one of the most important challenging procedure for the remote sensing community. In general, there are three methods utilized to validate LST values obtained from space, namely, the Temperature-based method (T-based), the Radiance-based method (R-based), and cross-validation [7]. The R-based method considers the satellite-derived LST and in-situ atmospheric profiles and LSEs as initial input parameters to simulate the TOA radiance using radiative transfer simulations at the moment of the satellite overpass [47]. The difference between the adjusted LST and the initial satellite-derived LST represents the accuracy of the retrieved LST [7]. The cross-validation method considers a well-validated LST product as a reference and compares the satellite-derived LST with the referenced (well-validated) LST derived from other satellites [7]. The T-based method, used by many researchers and also considered in this study, directly compares the satellite-based LST with ground-based LST measurements at the satellite overpass [46,48,49,50,51,52,53,54,55]. The main advantage of the T-based method is that it enables evaluating the radiometric quality of the satellite sensor and the performance of LST retrieval methods depending on atmospheric and emissivity parameters. However, the effectiveness of the T-based assessments relies largely on the accuracy of the ground-based LSTs and how well they represent the LST at the satellite pixel scale [7]. In addition, another issue that affects the correctness of the T-based validation method is the accuracy of calibration in TIR bands [56].

In this study, LST retrieval algorithms, namely, Radiative Transfer Equation (RTE) method [39,42], Single Channel Algorithm (SCA) [44] and Mono Window Algorithm (MWA) [43] were evaluated using Landsat 5 Thematic Mapper (TM), 7 Enhanced Thematic Mapper Plus (ETM+) and 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) data. Additionally, Split Window Algorithm (SWA) [45,46] was assessed for Landsat 8 OLI/TIRS data. Since LSE is one of the most important factors influencing the LST estimation reliability, the effects of different Normalized Difference Vegetation Index (NDVI)-based LSE models on LST accuracy were investigated in this study. In previous studies, many researchers have already examined the validation of different LST retrieval methods using Landsat data and in-situ LST measurements; however, they just considered one LSE model in the validation. Meng et al. [52] used the National Oceanic and Atmospheric Administration (NOAA) Joint Polar Satellite System (JPSS) Enterprise algorithm and a hybrid LSE model to retrieve LST from Landsat-8 data. Yu et al. [46] compared RTE, SWA, and SCA methods using Landsat 8 data and their NDVI-based LSE model reported in Section 4.3. Zhang et al. [57] utilized Sobrino et al.’s LSE model [58] and the SCA method for LST retrieval from Landsat 8 data. Zhang et al. [53] also investigated the accuracy of SCA using Landsat 8 imagery and Surface Radiation Budget Network (SURFRAD) measurements. Wang et al. [54] proposed a Practical Single-Channel Algorithm (PSCA) using Sobrino et al.’s LSE model. Sekertekin [59] used Skoković et al.’s LSE model [60] and compared RTE-based LST from Landsat 8 with SURFRAD measurements. As pointed out above, researchers generally focused on the validation of Landsat 8 derived LST images with in-situ measurements. However, the LST validation results of Landsat 5 TM and 7 ETM+ still remain insufficiently explored. Therefore, this study provides the LST validation results of Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS data examining different LST retrieval methods and LSE models.

This study aims to evaluate the performance of LST retrieval methods using different NDVI-based LSE models and data of old and current Landsat missions. The U.S. Geological Survey (USGS) has been producing and publishing LST products of Landsat missions considering LSE from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Emissivity Database (GED) data covering the US, Africa, Arabian Peninsula, Australia, Europe, and China; however, these LST products are geographically limited within the boundary of the North American Regional Reanalysis (NARR) grid, which is the climate data set used in the atmospheric correction algorithm [61,62]. Thus, it is still important to analyze Landsat data with different LST retrieval methods and over other climatological regions than North America. The ground-based measurement of LST requires accurate upwelling and downwelling thermal radiation measurements, and there are few stations in the world that measure these parameters to obtain in-situ LST. SURFRAD stations, established by the NOAA Office of United States (US) and located at different climatological regions of the US, are unique sources of information about in-situ LST over rural areas [63]. In this work, a total of forty-five Landsat scenes, fifteen images for each Landsat mission, acquired in the Spring-Summer-Autumn period over rural areas in the mid-latitude region in the Northern Hemisphere were obtained over five SURFRAD stations in the period of 2000–2019. Simultaneous in-situ LST data with satellite acquisitions, obtained from the correspondent SURFRAD station, were utilized for accuracy analyses. For the aim of the work, we developed an enhanced toolbox for automated LST retrieval from Landsat data by RTE, SCA, MWA, and SWA algorithms using different LSE models (Supplementary Materials). This toolbox is a first step to fill the gaps in the availability of different LST retrieval methods/LSE models in packaged Remote Sensing (RS) or Geographic Information System (GIS) software.

2. Data Set

2.1. Landsat Imagery

Landsat series of satellites have provided space-based moderate-resolution remote sensing data continuously for more than four decades. From 23 July 1972, in total, eight series of Landsat satellites were launched for Earth Observation (EO) purposes. Landsat 6 was the only satellite that failed to achieve orbit. The rest of the satellites have provided a unique resource for global change research and applications in agriculture, cartography, geology, forestry, regional planning, surveillance, and education over the last four decades. In this study, fifteen images for each Landsat series of 5, 7, and 8 were utilized for LST retrieval. The acquisition years of Landsat data range from 2000 to 2019, and only clear-sky images were considered. The selected dates ensure the presence of the in-situ data and the same number of images for the three Landsat missions. Landsat data can be downloaded through the USGS ’Earth Explorer’ website free of charge.

Landsat 5 TM and Landsat 7 ETM+ have six reflective bands (visible, near-infrared, and short-wavelength infrared, 30-m spatial resolution) and one band in the TIR region (Band 6). The thermal band has a native spatial resolution of 120-m and 60-m for TM and ETM+, respectively, but it is delivered by USGS at 30-m after cubic convolution resampling. The Landsat 8 OLI sensor has nine reflective bands with 30-m spatial resolution, and Landsat 8 TIRS sensor has two bands in the TIR region (Band 10 and Band 11). These thermal bands have a 100-m native spatial resolution but resampled and published at 30 m by USGS.

Appendix A reports the information of Landsat data utilized in the study (forty-five images from 2000 to 2019), as well as meteorological data (near-surface temperature and relative humidity) and NDVI value information corresponding to the acquisition times of Landsat data.

2.2. SURFRAD Data and Validation Sites

SURFRAD network was established by NOAA Office in 1993 in order to support climate-related researches over the US. SURFRAD stations have been measuring accurate, continuous, and long-term in-situ surface radiation budget [63]. The system became operational in 1995 with four stations, and, currently, eight SURFRAD stations are operating in different climatological regions of the US. Upwelling and downwelling components of solar and infrared radiation are the primary measurements. Besides, ancillary observations include direct and diffuse solar radiation, ultraviolet-B radiation, and meteorological parameters. Since SURFRAD stations provide unique in-situ LST information over rural sites, many researchers have used these data to validate satellite-based LST retrievals [15,46,52,53,56,64,65,66,67]. In this study, five SURFRAD stations were considered as ground-based stations, and

Table 1. Information about the Surface Radiation Budget Network (SURFRAD) experimental sites used in the study.

Site Name	Site Code	Latitude	Longitude	Elevation	Land Cover Type
Bondville, Illinois	BND	40.05° N	88.37° W	230 m	Cropland
Desert Rock, Nevada	DRA	36.62° N	116.02° W	1007 m	Open Shrub-lands
Fort Peck, Montana	FPK	48.31° N	105.10° W	634 m	Grassland
Goodwin Creek, Mississippi	GWN	34.26° N	89.87° W	98 m	Cropland/Natural Vegetation Mosaic
Penn. State Univ., Pennsylvania	PSU	40.72° N	77.93° W	376 m	Cropland

reports information regarding the SURFRAD experimental sites. We also analyzed the LST retrieval at the SURFRAD station of Table Mountain, Boulder, Colorado (TBL). However, the LST differences between the satellite and the ground were high due to its elevation and the heterogeneity of the land cover as also assessed in other studies [15,59,68,69]. Thus, TBL station was not considered in the analyses. SXF (Sioux Falls, South Dakota) and SGP (ARM Southern Great Plains Facility, Oklahoma) station data were not processed in this study.

The native spatial resolution of the Landsat thermal channels spans from 60 to 120 m, even if pixels are resampled at 30 m by USGS. The SURFRAD pyrgeometer used to measure the upwelling radiation is deployed at a 10-m high tower, producing an effective diameter of the field-of-view of about 40 m at the surface, i.e., roughly of the same order of the Landsat pixel size. Therefore, the Landsat pixel covering the SURFRAD instrument was selected for the comparison test.

3. LST Retrieval Methods

As previously pointed out, the following four LST retrieval methods will be considered: Mono Window Algorithm (MWA) [43], Single Channel Algorithm (SCA) [44], Radiative Transfer Equation (RTE) method and Split Window Algorithm (SWA) [45,46]. While the first three methods can be applied to Landsat 5 TM, 7 ETM+ and 8 OLI/TIRS data, the SWA is applicable only to Landsat 8 OLI/TIRS data, since it requires at least two TIR bands. The essential differences between these methods are in the mathematical formulation and the input parameters [70]. In addition to the emissivity and the atmospheric transmissivity common to all methods, MWA needs near-surface air temperature for the effective mean atmospheric temperature computation unlike other methods. Conversely, RTE and SCA require the upwelling and downwelling atmospheric radiances for LST retrieval. The sensitivity of the input parameters on LST retrieval methods is reported in Appendix D.

3.1. Mono Window Algorithm

Mono Window Algorithm (MWA) was developed by Qin et al. [43] for Landsat TM data. The method requires three main parameters, i.e., emissivity, atmospheric transmittance, and effective mean atmospheric temperature. LST values from MWA can be estimated as:

$$\begin{array}{c}{\mathrm{T}}_{\mathrm{s}}=\left\{\mathrm{a}·\left(1-\mathrm{C}-\mathrm{D}\right)+\left[\mathrm{b}·\left(1-\mathrm{C}-\mathrm{D}\right)+\mathrm{C}+\mathrm{D}\right]·\mathrm{T}-\mathrm{D}·{\mathrm{T}}_{\mathrm{a}}\right\}÷\mathrm{C}\\ \mathrm{a}=-67.355351, \mathrm{b}=0.458606, \mathrm{C}=\mathsf{\epsilon }\times \mathsf{\tau }, \mathrm{D}=(1-\mathsf{\tau })[1+(1-\mathsf{\epsilon })\times \mathsf{\tau }]\end{array}$$

(1)

where T_s is the LST in Kelvin, T is the at-sensor brightness temperature in Kelvin, T_a is the effective mean atmospheric temperature in Kelvin, τ is the atmospheric transmittance, ε represents LSE, a and b are the algorithm constants, C and D are the algorithm parameters calculated using LSE and transmittance. A detailed description of the computations of the T, T_a and τ parameters adopted in this work, are reported in Appendix B. The different LSE models tested in this work will be described in Section 4.

3.2. Single-Channel Algorithm

Jiménez-Muñoz et al. [44] introduced a revision of the Single-Channel Algorithm (SCA) to retrieve LST from Landsat TIR data. Considering SCA, LST (T_s) can be computed using the following general equation:

$${\mathrm{T}}_{\mathrm{s}}=\mathsf{\gamma }\left[\frac{1}{\mathsf{\epsilon }}\left({\mathsf{\psi }}_1{\mathrm{L}}_{\mathrm{sen}}{+\mathsf{\psi }}_2\right){+\mathsf{\psi }}_3\right]+\mathsf{\delta }$$

(2)

where ε is the LSE, L_sen is the at-sensor radiance of thermal band, ψ₁, ψ₂, and ψ₃ are atmospheric functions, and γ, δ are two parameters given by:

$$\mathsf{\gamma }\approx \frac{{\mathrm{T}}^2}{{\mathrm{b}}_{\mathsf{\gamma }}{\mathrm{L}}_{\mathrm{sen}}}$$

(3)

$$\mathsf{\delta }\approx \mathrm{T}-\frac{{\mathrm{T}}^2}{{\mathrm{b}}_{\mathsf{\gamma }}}$$

(4)

where b_γ = c₂/λ_i with c₂ = 14,387.7 µm∙K and λ_i is the effective band wavelength for band i, which is defined as:

$${\mathsf{\lambda }}_{\mathrm{i}}=\frac{{\int }_{{\mathsf{\lambda }}_{1,\mathrm{i}}}^{{\mathsf{\lambda }}_{2,\mathrm{i}}}{\mathsf{\lambda }\mathrm{f}}_{\mathrm{i}}\left(\mathsf{\lambda }\right)\mathrm{d}\mathsf{\lambda }}{{\int }_{{\mathsf{\lambda }}_{1,\mathrm{i}}}^{{\mathsf{\lambda }}_{2,\mathrm{i}}}{\mathrm{f}}_{\mathrm{i}}\left(\mathsf{\lambda }\right)\mathrm{d}\mathsf{\lambda }}$$

(5)

where f_i(λ) is the spectral response function for the corresponding band. λ_{1, i} and λ_{2, i} are the lower and upper boundary of f_i(λ), respectively. The value of b_γ is equal to 1256 K and 1277 K for Band 6 of Landsat 5 and Landsat 7, respectively; for Band 10 and Band 11 of Landsat 8, it is equal to 1320 K and 1199 K, respectively.

Atmospheric functions ψ₁, ψ₂, and ψ₃ are defined as:

$${\mathsf{\psi }}_1=\frac{1}{\mathsf{\tau }}{ ; \mathsf{\psi }}_2{=-\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }-\frac{L_{\lambda }^{\uparrow }}{\mathsf{\tau }}{;\mathsf{\psi }}_3{=\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }$$

(6)

where ${\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$ (W∙m⁻²∙sr⁻¹∙μm⁻¹) is upwelling or atmospheric path radiance, ${\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }$ (W∙m⁻²∙sr⁻¹∙μm⁻¹) is downwelling or sky radiance. In this study, the atmospheric parameters τ, ${\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$ and ${\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }$ used for the ψ₁, ψ₂, and ψ₃ computation are reported in Appendix B.

3.3. Radiative Transfer Equation Method

A straightforward method to retrieve LST from a single TIR band is the inversion of the radiative transfer equation (RTE) according to the following expressions:

$${\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}=\left[{\mathsf{\epsilon }\mathrm{B}}_{\mathsf{\lambda }}\left({\mathrm{T}}_{\mathrm{s}}\right)+\left(1-\mathsf{\epsilon }\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }\right]{\mathsf{\tau }+\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$$

(7)

where ${\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}$ (W∙m⁻²∙sr⁻¹∙μm⁻¹) is at-sensor registered radiance of the related thermal band, B_λ (W∙m⁻²∙sr⁻¹∙μm⁻¹) is the blackbody radiance. Blackbody radiance (B_λ) at a temperature of T_s can be obtained by inverting the Equation (7):

$${\mathrm{B}}_{\mathsf{\lambda }}\left({\mathrm{T}}_{\mathrm{s}}\right)=\frac{{\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}{-\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }-\mathsf{\tau }\left(1-\mathsf{\epsilon }\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }}{\mathsf{\tau }\mathsf{\epsilon }}$$

(8)

and, finally, T_s can be obtained by inverting Planck’s law as:

$${\mathrm{T}}_{\mathrm{s}}=\frac{{\mathrm{K}}_2}{\ln \left(\frac{{\mathrm{K}}_1}{\frac{{\mathrm{L}}_{\mathsf{\lambda }}^{\mathrm{sen}}{-\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }-\mathsf{\tau }\left(1-\mathsf{\epsilon }\right){\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }}{\mathsf{\tau }\mathsf{\epsilon }}}+1\right)}$$

(9)

where K₁ and K₂ are calibration constants for Landsat data reported in Appendix B.

3.4. Split-Window Algorithm

In previous studies, various Split Window Algorithms (SWAs) have been introduced for different sensors [4,71,72,73,74] and detailed information of SWAs is reported in [7]. Among the different SWAs in the literature, in this study we considered SWA developed by Mao et al. [45] with coefficients re-parameterized by Yu et al. [46], corresponding to the Landsat 8 TIRS’ spectral response curve. The USGS recommended not to use Band 11 of Landsat 8 for LST retrieval due to the large calibration uncertainty [75]. However, some researchers claimed that they obtained satisfactory results via SWA [46,76]. Thus, we also analyze and present SWA results in this study. According to SWA, LST (T_s) can be calculated using the following equations:

$${\mathrm{T}}_{\mathrm{s}}{=\mathrm{T}}_{10}{+\mathrm{B}}_1\left({\mathrm{T}}_{10}{-\mathrm{T}}_{11}\right){+\mathrm{B}}_0$$

(10)

$${\mathrm{B}}_0=\frac{{\mathrm{C}}_{11}\left({1-\mathrm{A}}_{10}{-\mathrm{C}}_{10}\right){\mathrm{L}}_{10}{-\mathrm{C}}_{10}\left({1-\mathrm{A}}_{11}{-\mathrm{C}}_{11}\right){\mathrm{L}}_{11}}{{\mathrm{C}}_{11}{\mathrm{A}}_{10}{-\mathrm{C}}_{10}{\mathrm{A}}_{11}}$$

(11)

$${\mathrm{B}}_1=\frac{{\mathrm{C}}_{10}}{{\mathrm{C}}_{11}{\mathrm{A}}_{10}{-\mathrm{C}}_{10}{\mathrm{A}}_{11}}$$

(12)

$${\mathrm{A}}_{10}{=\mathsf{\epsilon }}_{10}{\mathsf{\tau }}_{10}$$

(13)

$${\mathrm{A}}_{11}{=\mathsf{\epsilon }}_{11}{\mathsf{\tau }}_{11}$$

(14)

$${\mathrm{C}}_{10}=\left({1-\mathsf{\tau }}_{10}\right)\left(1+\left({1-\mathsf{\epsilon }}_{10}\right){\mathsf{\tau }}_{10}\right)$$

(15)

$${\mathrm{C}}_{11}=\left({1-\mathsf{\tau }}_{11}\right)\left(1+\left({1-\mathsf{\epsilon }}_{11}\right){\mathsf{\tau }}_{11}\right)$$

(16)

where ${\mathsf{\epsilon }}_{10}$ and ${\mathsf{\epsilon }}_{11}$ represent LSE for Band 10 and 11, respectively, ${\mathsf{\tau }}_{10}$ and ${\mathsf{\tau }}_{11}$ the atmospheric transmittance for Band 10 and 11, respectively, calculated as reported in Appendix B. L₁₀ and L₁₁ can be computed from

Table 2. Linear regression coefficients for the parameters L₁₀ and L₁₁.

TIR Bands	Range	a	B (K)
Band 10	−10–20 °C	0.4087	−55.58
	20–50 °C	0.4464	−66.61
Band 11	−10–20 °C	0.4442	−59.85
	20–50 °C	0.4831	−71.23

within a specific brightness temperature range for Band 10 (T₁₀) and Band 11 (T₁₁), respectively. In Table 2, “a” is the slope and B(K) is the intercept of linear regression. For example, if the brightness temperature of B₁₀ ranges between 20 and 50 °C, L₁₀ can be calculated by 0.4464 * T₁₀ − 66.61.

4. Land Surface Emissivity (LSE) Models

Surface emissivity stands for the surface ability that transforms heat energy into radiant energy [36]. LSE (ε) is one of the key parameters to retrieve accurate LST from remotely sensed imagery. Semi-Empirical Methods (SEMs), Physically-Based Methods (PBMs), and multi-channel Temperature/Emissivity Separation (TES) methods are three distinctive categories for LSE retrieval from space [7]. PBMs and multi-channel TES methods are not operational for Landsat data to obtain LSE due to the limitations presented in many studies, such as the requirement of more than two TIR bands or nighttime images [7,46,76,77,78,79]. SEMs contain the Classification Based Emissivity Method (CBEM) [74,80] and the NDVI Based Emissivity Method (NBEM) [81,82], which are suitable for LSE estimation from Landsat data. The CBEM generates an LSE image from a classified image by applying an emissivity value for each class. However, CBEM is not practical since it requires a good knowledge of the study area and emissivity measurements on the surfaces representative of the different classes [70]. NDVI-based methods are operative and the most commonly utilized LSE retrieval methods since they are easy to apply and presenting satisfying results [36,58,83]. Li et al. [7] presented a detailed study showing the advantages, disadvantages, and limitations of different LSE models for LST retrieval from satellite data. Considering the study of Li et al. [7] and other researches, a state-of-the-art table showing different LSE categories and models, as well as the correspondent satellite data used is reported (

Table 3. The state-of-art table showing different LSE categories, LSE models, and the corresponding satellite missions used.

Category	Surface Emissivity Determination Methods	References	Platform
Semi-Empirical Methods (SEMs)	Classification-based emissivity method (CBEM)	[80]	MSG1/SEVIRI
		[84]	MODIS
	NDVI-based emissivity method (NBEM)	[83]	NOAA/AVHRR Landsat TM
		[82]	NOAA/AVHRR Landsat TM
		[81]	NOAA/AVHRR
		[85]	TERRA/MODIS
		[58]	ENVISAT/AATSR MSG1/SEVIRI Landsat TM
		[60]	Landsat 8
		[46]	Landsat 8
		[86]	MODIS
		[87]	TERRA/MODIS
		[76]	Landsat 8
Multi-channel TES methods	The two-temperature method (TTM)	[88]	TIMS
		[89]	MSG/SEVIRI
		[90]	MSG/SEVIRI
	Grey-body emissivity (GBE) method	[91]	TIMS
	The iterative spectrally smooth temperature emissivity separation (ISSTES) method	[92,93]	Hyperspectral infrared data
	The emissivity bounds method (EBM)	[94]	TIMS
	Reference channel method (RCM)	[95]	multispectral aircraft scanner data
	TES method	[40]	ASTER
		[58]	ASTER Airborne Hyperspectral Scanner (AHS)
	Temperature-independent spectral indices (TISI) based methods	[33]	NOAA/AVHRR
		[96]	NOAA/AVHRR
		[97]	TERRA/MODIS
		[98]	MSG-SEVIRI
Physically-based methods (PBMs)	Physics-based day/night (D/N) method	[99]	TERRA/MODIS
	Two-step physical retrieval method (TSRM)	[100,101]	TERRA/MODIS
		[102]	AQUA/AIRS

).

As presented in Table 3, there are six NDVI-based models introduced for Landsat data, specifically three for Landsat TM and three for Landsat 8 OLI/TIRS. Therefore, we investigated the effect of these six LSE models on the accuracy of LST retrieval methods. Details about LSE models are presented in the following sub-sections. The sensitivity of the LSE on LST retrieval methods is reported in Appendix D.

4.1. LSE Model of Van de Griend and Owe

This model was applied to LST retrieval methods of all Landsat series (Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS). Van de Griend and Owe [83] proposed a logarithmic approach for an LSE retrieval model based on NDVI ranging from 0.157 to 0.727. NDVI is obtained using the Near-Infrared (NIR) and Red (R) bands—the calculation steps of NDVI for Landsat 5, 7, and 8, are presented in Appendix C. The proposed model is given by:

$$\mathsf{\epsilon }=1.0094+0.047\ln \left(\mathrm{NDVI}\right)$$

(17)

4.2. LSE Model of Valor and Caselles

This model was applied to LST retrieval methods of all Landsat series (Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS). Valor and Caselles [82] proposed a theoretical model that relates the emissivity to the NDVI of a given surface by:

$${\mathsf{\epsilon }=\mathsf{\epsilon }}_{\mathrm{v}}{\mathrm{P}}_{\mathrm{v}}{+\mathsf{\epsilon }}_{\mathrm{s}}\left({1-\mathrm{P}}_{\mathrm{v}}\right)+4\langle \mathrm{d}\mathsf{\epsilon }\rangle {\mathrm{P}}_{\mathrm{v}}\left({1-\mathrm{P}}_{\mathrm{v}}\right)$$

(18)

${\mathsf{\epsilon }}_{\mathrm{v}}$ and ${\mathsf{\epsilon }}_{\mathrm{s}}$ represent the emissivity of vegetation and soil, respectively. $\langle \mathrm{d}\mathsf{\epsilon }\rangle$ is a term accounting for the cavity effect, which depends on the surface geometry. Pv (also referred to as fractional vegetation cover, FVC) is the proportion of vegetation calculated as [103]:

$${\mathrm{P}}_{\mathrm{v}}= {\left[\frac{{\mathrm{NDVI} - \mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }{ - \mathrm{NDVI}}_{\min }}\right]}^2$$

(19)

where NDVI_max = 0.5 and NDVI_min = 0.2 in a global situation [70]. As Valor and Caselles [82] suggested, ${\mathsf{\epsilon }}_{\mathrm{v}}$ and ${\mathsf{\epsilon }}_{\mathrm{s}}$ as 0.985 and 0.960, respectively, for unknown emissivity and vegetation structures, we also regarded these emissivity values in the calculation. Besides, they calculated the mean value for $\langle \mathrm{d}\mathsf{\epsilon }\rangle$ term as 0.015, and we utilized this value in LSE retrieval with this model. The final version of the LSE model can be given by:

$$\mathsf{\epsilon }=0.985{\mathrm{P}}_{\mathrm{v}}+0.960\left({1-\mathrm{P}}_{\mathrm{v}}\right)+0.06{\mathrm{P}}_{\mathrm{v}}\left({1-\mathrm{P}}_{\mathrm{v}}\right)$$

(20)

4.3. NDVI Threshold (NDVI^THM)-Based LSE Models

Sobrino et al. [58], Skoković et al. [60], Yu et al. [46], and Li and Jiang [76] estimated LSE from NDVI threshold (NDVI^THM) values considering three different cases as presented in Equation (21). In the first case (NDVI < 0.2), the pixel is considered as bare soil, and the emissivity is obtained from the reflectance values in the red region. In the second case (0.2 ≤ NDVI ≤ 0.5), the pixel is composed of a mixture of bare soil and vegetation, and in the third case (NDVI > 0.5), the pixels with NDVI values higher than 0.5 are considered as fully vegetated areas.

$$\mathsf{\epsilon }=\{\begin{array}{cc}{\mathrm{a}}_{\mathrm{i}}{\mathsf{\rho }}_{\mathrm{R}}+{\mathrm{b}}_{\mathrm{i}}\hfill & \hfill \mathrm{NDVI}<0.2\\ {\mathsf{\epsilon }}_{\mathrm{v}}+{\mathsf{\epsilon }}_{\mathrm{s}}(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon },\mathrm{d}\mathsf{\epsilon }=(1-{\mathsf{\epsilon }}_{\mathrm{s}})(1-{\mathrm{P}}_{\mathrm{v}})\mathrm{F}{\mathsf{\epsilon }}_{\mathrm{V}}\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ {\mathsf{\epsilon }}_{\mathrm{V}}+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$$

(21)

In Equation (21), ${\mathsf{\rho }}_{\mathrm{R}}$ is the reflectance value of the red band, ${\mathrm{a}}_{\mathrm{i}}$ and ${\mathrm{b}}_{\mathrm{i}}$ are estimated from an empirical relationship between the red band reflectance and Moderate Resolution Imaging Spectroradiometer (MODIS) emissivity library. ${\mathsf{\epsilon }}_{\mathrm{v}}$ and ${\mathsf{\epsilon }}_{\mathrm{s}}$ are the soil and vegetation emissivity, respectively. $\mathrm{d}\mathsf{\epsilon }$ is the cavity effect due to surface roughness as in the previous model ($\mathrm{d}\mathsf{\epsilon }$ = 0 for flat surfaces). F is a geometrical shape factor assumed as the mean value of 0.55 [70].

Table 4. The expressions of Normalized Difference Vegetation Index (NDVI) threshold models used in this study.

Sensor	LSE Equations	Reference
Landsat 5 TM and 7 ETM+ (Band 6)	$\mathsf{\epsilon }=\{\begin{array}{cc}0.979-0.035{\mathsf{\rho }}_{\mathrm{R}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.004{\mathrm{P}}_{\mathrm{v}}+0.986\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.99\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Sobrino et al. [58]
Landsat 8 TIR1 (Band 10)	$\mathsf{\epsilon }=\{\begin{array}{cc}0.979-0.046{\mathsf{\rho }}_{\mathrm{R}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.987{\mathrm{P}}_{\mathrm{v}}+0.971(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.987+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Skoković et al. [60]
Landsat 8 TIR1 (Band 11)	$\mathsf{\epsilon }=\{\begin{array}{cc}0.982-0.027{\mathsf{\rho }}_{\mathrm{R}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.989{\mathrm{P}}_{\mathrm{v}}+0.977(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.989+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Skoković et al. [60]
Landsat 8 TIR1 (Band 10)	$\mathsf{\epsilon }=\{\begin{array}{cc}0.973-0.047{\mathsf{\rho }}_{\mathrm{R}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.9863{\mathrm{P}}_{\mathrm{v}}+0.9668(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.9863+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Yu et al. [46]
Landsat 8 TIR1 (Band 11)	$\mathsf{\epsilon }=\{\begin{array}{cc}0.984-0.0026{\mathsf{\rho }}_{\mathrm{R}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.9896{\mathrm{P}}_{\mathrm{v}}+0.9747(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.9896+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Yu et al. [46]
Landsat 8 TIR1 (Band 10)	$\mathsf{\epsilon }=\{\begin{array}{cc}{\mathrm{a}}_{\mathrm{l}\mathrm{i}}+{\displaystyle \sum _{\mathrm{j}=2}^7}{\mathrm{a}}_{\mathrm{j}\mathrm{i}}{\mathsf{\rho }}_{\mathrm{j}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.982{\mathrm{P}}_{\mathrm{v}}+0.971(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.982+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Li and Jiang [76]
Landsat 8 TIR1 (Band 11)	$\mathsf{\epsilon }=\{\begin{array}{cc}{\mathrm{a}}_{\mathrm{l}\mathrm{i}}+{\displaystyle \sum _{\mathrm{j}=2}^7}{\mathrm{a}}_{\mathrm{j}\mathrm{i}}{\mathsf{\rho }}_{\mathrm{j}}\hfill & \hfill \mathrm{NDVI}<0.2\\ 0.984{\mathrm{P}}_{\mathrm{v}}+0.976(1-{\mathrm{P}}_{\mathrm{v}})+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill 0.2\le \mathrm{NDVI}\le 0.5\\ 0.984+\mathrm{d}\mathsf{\epsilon }\hfill & \hfill \mathrm{NDVI}>0.5\end{array}$	Li and Jiang [76]

presents the expressions of NDVI^THM for all models mentioned above.

In Li and Jiang’s LSE model, ${\mathsf{\rho }}_{\mathrm{j}}$ is the apparent reflectance in the OLI band j and ${\mathrm{a}}_{1\mathrm{i}}$–${\mathrm{a}}_{7\mathrm{i}}$ are coefficients obtained from [76]. The LSE models of Band 11 were just utilized in the SWA method. It is important to point out again that the USGS announced caution in the use of Band 11 of Landsat 8 due to the calibration uncertainties [75]. However, some researchers published satisfactory results by using SWA [46,76].

5. LST Computation Using Ground-Based SURFRAD Data

As stated in Section 2.2, image-based LST results were validated using the data of five ground-based SURFRAD stations. Since these stations do not provide LST measurements directly, LST is calculated from the upwelling and downwelling longwave radiation measurements using the following Equation (22) with regard to the Stefan–Boltzmann law:

$$\mathrm{LST}={\left[\frac{{\mathrm{F}}_{\mathsf{\lambda }}^{\uparrow }-\left({1-\mathsf{\epsilon }}_{\mathrm{b}}\right)·{\mathrm{F}}_{\mathsf{\lambda }}^{\downarrow }}{{\mathsf{\epsilon }}_{\mathrm{b}}·\mathsf{\sigma }}\right]}^{1/4}$$

(22)

where ${\mathrm{F}}_{\mathsf{\lambda }}^{\uparrow }$ and ${\mathrm{F}}_{\mathsf{\lambda }}^{\downarrow }$ represent upwelling and dowelling thermal infrared (3–50 μm) irradiance in W/m², respectively, measured during satellite passages. $\mathsf{\sigma }$ is the Stefan–Boltzmann constant (5.670367 × 10⁻⁸ W∙m⁻²∙K⁻⁴), and ${\mathsf{\epsilon }}_{\mathrm{b}}$ represents the broadband longwave surface emissivity, which is not measured by the station instruments. In previous studies on SURFRAD stations [53,56,65], the broadband emissivity was computed as reported in [104,105] by regression from narrowband emissivity of MODIS thermal bands, which are available through the MODIS monthly emissivity data set. The results in [104,105] proved that the longwave broadband emissivity for the SURFRAD sites could be considered 0.97, as also assumed in [64] and [53].

Therefore, in this study, the broadband emissivity was assumed 0.97. This assumption impacts only the SURFRAD LST estimation, not the satellite-derived estimation. Heidinger et al. [66] investigated the impact of changing the assumed broadband emissivity from 0.97 to 0.98 on the SURFRAD LST observation. The results indicated that a 0.01 error in broadband emissivity produces a SURFRAD LST error that rarely exceeds 0.25 K. In addition, Wang and Liang [105] proved that the sensitivity of the SURFRAD LST to broadband emissivity ranged from 0.1K/0.01 to 0.35K/0.01, which means the accuracy of LST varies between 0.1 K to 0.4 K when the broadband emissivity error is about ±0.01. While this error is not negligible, it does not appear to be a dominant source of uncertainty in the SURFRAD-based performance metrics considering the magnitude of the other uncertainties [66].

Satellite-based LST and SURFRAD-based LST were compared using statistical criteria, namely, the Root Mean Square Error (RMSE). RMSE, in Equation (23), is a widely used statistical metric evaluating the efficiency of the models.

$$\mathrm{RMSE}=\sqrt{\frac{\sum {\left[{\mathrm{T}}_{\mathrm{SAT}} {- \mathrm{T}}_{\mathrm{SURF}}\right]}^2}{\mathrm{n}}}$$

(23)

where ${\mathrm{T}}_{\mathrm{SAT}}$ and ${\mathrm{T}}_{\mathrm{SURF}}$ are the satellite-based LST and SURFRAD-based LST, respectively, and n represents the pixel count.

6. Results

In order to compare the results of LST retrieval methods, fifteen images of each Landsat mission (Landsat 5, 7, and 8), a total of forty-five images, were utilized in this study. MWA, RTE, and SCA were performed for all satellite data, whilst SWA was only utilized with Landsat 8 data due to the requirement of two TIR bands. The values of the atmospheric and model parameters used in LST retrieval methods for the 45 images are reported in Appendix E. Furthermore, the effects of different LSE models on the accuracy of the LST retrieval methods were investigated.

6.1. Results of LST Algorithms and LSE Models Derived from Landsat 5 TM

In

Table 5. Validation results of the Land Surface Temperature (LST) retrieval methods for Landsat 5 Thematic Mapper (TM) data based on different Land Surface Emissivity (LSE) models. The best result is in bold.

Landsat Mission	Emissivity Method	LST Retrieval Method	RMSE (K)
Landsat 5 TM	Van De Griend & Owe (1993)	MWA	4.89
		RTE	4.96
		SCA	5.22
	Valor & Caselles (1996)	MWA	2.93
		RTE	3.25
		SCA	3.46
	Sobrino et al. (2008)	MWA	2.41
		RTE	2.35
		SCA	2.47

, the accuracy of LST retrieval methods for Landsat 5 TM data is presented based on the different LSE models. Considering the LSE models, Sobrino et al.’s model provided the best results for all three LST retrieval methods. Valor and Caselles’ model was the second LSE model presenting satisfying results for LST estimation, whilst Van De Griend and Owe’s LSE model provided very high RMSE. The SURFRAD test assessed whether the RTE method was a bit better than MWA and SCA, with a lower RMSE value of 2.35 K.

6.2. Results of LST Algorithms and LSE Models Derived from Landsat 7 ETM+

In

Table 6. Validation results of LST retrieval methods for Landsat 7 ETM+ data based on different LSE models. The best result is in bold.

Landsat Mission	Emissivity Method	LST Retrieval Method	RMSE (K)
Landsat 7 ETM+	Van De Griend & Owe (1993)	MWA	9.10
		RTE	8.18
		SCA	9.51
	Valor & Caselles (1996)	MWA	4.64
		RTE	4.95
		SCA	5.25
	Sobrino et al. (2008)	MWA	2.24
		RTE	2.48
		SCA	2.77

, validation results for Landsat 7 ETM+ data are reported. Again, Sobrino et al.’s model provided the best results whilst both Valor & Caselles’ model and Van De Griend & Owe’s LSE models presented much higher RMSE values. Although all LST retrieval methods with Sobrino et al.’s LSE model presented good results when using Landsat 7 ETM+ data, the results revealed that MWA provided slightly better results than RTE and SCA (RMSE value equals to 2.24 K).

6.3. Results of LST Algorithms and LSE Models Derived from Landsat 8 OLI/TIRS Data

In

Table 7. Validation results of LST retrieval methods for Landsat 8 OLI/TIRS data based on different LSE models. The best result is in bold.

Landsat Mission	Emissivity Method	LST Retrieval Method	RMSE (K)
Landsat 8 OLI/TIRS	VanDeGriend & Owe (1993)	MWA	4.24
		RTE	4.28
		SCA	4.53
	Valor & Caselles (1996)	MWA	5.16
		RTE	4.21
		SCA	5.11
	Sobrino et al. (2008)	MWA	2.52
		RTE	2.85
		SCA	2.94
	Skoković et al. (2014)	MWA	2.73
		RTE	3.01
		SCA	3.11
		SWA	2.79
	Yu et al. (2014)	MWA	2.79
		RTE	3.07
		SCA	3.18
		SWA	3.02
	Li & Jiang (2018)	MWA	2.85
		RTE	3.11
		SCA	3.22
		SWA	2.94

, the accuracy of LST retrieval methods for Landsat 8 OLI/TIRS data is assessed considering six LSE models. Again, Sobrino et al.’s model presented the best results as for Landsat 5 TM and 7 ETM+. The three LSE models Skoković et al., Yu et al., and Li & Jiang, specifically proposed in the literature for Landsat 8 data, also showed good results for all LST retrieval methods (the highest RMSE is 3.22 K). The LSE models of Valor & Caselles and Van De Griend & Owe presented the worst RMSE values again. This test suggests that MWA with Sobrino et al.’s LSE model provides better results than RTE and SCA, with a lower RMSE value (2.52 K). Considering the SWA method, requiring two emissivity images corresponding to the two TIR bands of Landsat 8 (Band 10 and 11), the three LSE models (Skoković et al., Yu et al., and Li & Jiang) provided satisfactory results. However, Skoković et al.’s LSE model demonstrated a slightly better RMSE value.

6.4. Comparison of LST Retrieval Algorithms Considering All Landsat Missions

In Section 6.1, Section 6.2 and Section 6.3, LST results were analyzed by dividing the sensor types of Landsat missions. In this section, the accuracies of the LST retrieval algorithms with respect to the ground-based LST data were evaluated considering the best LSE Model (Sobrino et al.’s) and all Landsat missions. This comparison can be significant for users who will conduct time-series analyses of LST over rural areas using the data of all Landsat missions. Since SWA was only used with Landsat 8 data, but it is not the best retrieval method, it was not considered in this section for the comparison purposes. In Figure 1a–c, MWA-based, RTE-based, and SCA-based LST results derived from Landsat 5, 7, and 8 data were compared with SURFRAD LST, respectively. The RMSE was 2.39 K for MWA, 2.57 K for RTE, and 2.73 K for SCA. The average biases (ground LST-satellite LST) for MWA, RTE, and SCA are −0.72 K, −1.63 K, and −1.81 K, respectively. Moreover, the error standard deviations (STD) are 2.28 K, 1.98 K, and 2.05 K for MWA, RTE, and SCA, respectively. Even though MWA has a slightly greater error STD, overall, the RMSE for RTE and SCA is higher due to the greater biases. This in-situ test over different rural areas showed that MWA, RTE, and SCA can provide good results with Sobrino et al.’s LSE model, and MWA presented slightly better performance. Figure 1 and biases show that although all methods tend to overestimate LST slightly, the MWA overestimation is lower. The different results, especially the bias, can be ascribed to the accuracy of the input parameters: as reported previously, RTA and SCA have the same input parameters, whilst MWA uses T_a instead of ${\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$ and ${\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }$. We must also consider the different formulation of the methods: since SCA is derived from a mathematical approximation of RTE [41], it is expected they provide similar results.

6.5. Analysis of Spatio-Temporal and Seasonal LST Variations Between LST Retrieval Methods

As stated in previous sections, Sobrino et al.’s LSE model provided the best performance on all LST retrieval methods. In this section, we focus on investigating the spatio-temporal and seasonal relationship between the LST retrieval methods using the Sobrino et al.’s LSE model. This analysis does not represent an evaluation of the best LST retrieval method, which was assessed in previous sections by a test with SURFRAD data, but it suggests us if there is similarity or not between the three methods (MWA, RTE, and SCA) with the same LSE model under changing emissivity values due to seasonal variations. In this analysis, we consider the LST retrieval methods in a rural area of 6 km × 6 km (40,000 pixels) centered on the SURFRAD station of each Landsat image. Therefore, a total of 45 × 3 LST sub-images (45 Landsat images × 3 LST methods) were investigated. Figure 2 shows an example of Landsat 8 LST image over the 6 × 6 km² rural area, acquired on 27 April 2018 and covering the BND station, for the three methods used in this analysis.

To analyze the spatio-temporal and seasonal LST variations among the LST retrieval methods, the 45 Landsat images were categorized into three seasons: spring (18 images), summer (13), and autumn (14), and Root Mean Square (RMS) differences were calculated for each image. Winter images were not available due to cloudy conditions (see

Table A1. Information about Landsat data used in the study (45 images from 2000 to 2019). The last column shows the SURFRAD station included within the satellite scene. SURFRAD codes can be found in: https://www.esrl.noaa.gov/gmd/grad/surfrad/sitepage.html or in Table A1.

Sensor	Scene ID	Scene Acquisition Date and Time (UTC)	Path-Row	To (°C)	RH (%)	NDVI Value	SURFRAD Station
LANDSAT 5 TM	LT50230322007167PAC01	16/06/2007–16:29	23–32	30.8	35.3	0.348	BND
	LT50220322008243GNC01	30/08/2008–16:15	22–32	25.5	51.2	0.672
	LT50230322010255PAC01	12/09/2010–16:26	23–32	24.8	32.2	0.438
	LT50400352006267PAC01	24/09/2006–18:15	40–35	21.8	14.5	0.074	DRA
	LT50400352007142PAC01	22/05/2007–18:16	40–35	20.7	9.6	0.075
	LT50400352011281PAC01	08/10/2011–18:10	40–35	16.8	31.3	0.088
	LT50350262006136PAC01	16/05/2006–17:39	35–26	23.8	26.5	0.228	FPK
	LT50350262008238PAC01	25/08/2008–17:32	35–26	30.2	35.1	0.170
	LT50360262011253PAC01	10/09/2011–17:43	36–26	27.6	28.6	0.330
	LT50230362002249LGS01	06/09/2002–16:11	23–36	29.4	55.7	0.566	GWN
	LT50230362008218PAC01	05/08/2008–16:23	23–36	30.6	60.2	0.578
	LT50230362011242PAC01	30/08/2011–16:26	23–36	31.2	33.6	0.552
	LT50160322003267GNC02	24/09/2003–15:30	16–32	15.9	57.5	0.698	PSU
	LT50160322008233GNC01	20/08/2008–15:39	16–32	19.5	50.2	0.674
	LT50160322009139GNC01	19/05/2009–15:40	16–32	15.9	29.9	0.390
LANDSAT 7 ETM+	LE70230322000284EDC00	10/10/2000–16:26	23–32	13.0	35.5	0.515	BND
	LE70230322001254EDC00	11/09/2001–16:24	23–32	24.1	40.0	0.535
	LE70220322002186EDC00	05/07/2002–16:18	22–32	30.7	51.4	0.097
	LE70400352001165EDC00	14/06/2001–18:12	40–35	23.3	12.8	-0.025	DRA
	LE70400352001213EDC00	01/08/2001–18:11	40–35	32.5	10.9	-0.030
	LE70400352002168EDC00	17/06/2002–18:10	40–35	30.9	6.9	-0.027
	LE70350262000112EDC00	21/04/2000–17:39	35–26	21.2	22.6	-0.068	FPK
	LE70360262001217EDC00	05/08/2001–17:43	36–26	28.2	32.8	-0.068
	LE70350262002181EDC00	30/06/2002–17:36	35–26	23.3	32.5	-0.073
	LE70220362000117EDC00	26/04/2000–16:23	22–36	18.9	43.4	0.300	GWN
	LE70230362000220EDC00	07/08/2000–16:28	23–36	32.5	55.1	0.258
	LE70220362001167EDC00	16/06/2001–16:21	22–36	27.3	44.9	0.389
	LE70160322000091EDC00	31/03/2000–15:45	16–32	7.8	37.7	-0.029	PSU
	LE70160322002192EDC00	11/07/2002–15:41	16–32	17.9	44.4	0.375
	LE70160322002256EDC00	13/09/2002–15:40	16–32	21.5	37.3	0.250
LANDSAT 8 OLI/TIRS	LC80230322013247LGN01	04/09/2013–16:38	23–32	23.9	57.2	0.621	BND
	LC80230322018101LGN00	11/04/2018–16:35	23–32	12.8	57.2	0.421
	LC80230322018117LGN00	27/04/2018–16:35	23–32	15.2	32.4	0.622
	LC80400352017121LGN00	01/05/2017–18:22	40–35	23.5	14.0	0.110	DRA
	LC80400352018124LGN00	04/05/2018–18:21	40–35	26.4	14.6	0.108
	LC80400352018236LGN00	24/08/2018–18:22	40–35	32.8	8.7	0.088
	LC80350262017198LGN00	17/07/2017–17:48	35–26	24.7	22.1	0.213	FPK
	LC80360262018160LGN00	09/06/2018–17:53	36–26	27.5	51.2	0.370
	LC80350262018249LGN00	06/09/2018–17:47	35–26	21.8	38.4	0.232
	LC80220362016281LGN01	07/10/2016–16:32	22–36	27.5	44.0	0.410	GWN
	LC80220362017251LGN00	08/09/2017–16:32	22–36	22.5	44.8	0.626
	LC80220362018094LGN00	04/04/2018–16:31	22–36	8.6	43.1	0.397
	LC80160322015124LGN01	04/05/2015–15:51	16–32	24.3	24.1	0.365	PSU
	LC80160322016111LGN01	20/04/2016–15:51	16–32	15.2	15.2	0.512
	LC80160322019263LGN00	20/09/2019–15:53	16–32	20.6	53.8	0.637

). This is not an accuracy test as the one performed by the SURFRAD LST ground measurements, but it is an image-based analysis to highlight the differences among the three retrieval methods. Therefore, the RMS difference is used instead of the RMSE. In the 6 × 6 km² selected areas, the minimum LST values from satellite data for spring, summer, and autumn are 280.29 K, 281.70 K, and 284.67 K, respectively; the maximum are 321.88 K, 330.67 K, and 317.95 K, respectively. Figure 3 shows the box-plot graph presenting the seasonal RMS differences between the LST retrieval methods. The box-plot is used to display distributional characteristics of data [106]. The box-plot information, reported in Figure 3 by numbers, is the minimum (1) (the lowest data point excluding any outliers), first quartile (2), median (3), third quartile (4) and maximum (5) (the largest data point excluding any outliers). The cross “x” in the boxes refers to the mean value of the data set, and the points outside the minimum, and maximum values are assumed as outliers. Concerning Figure 3, blue box-plots represent the RMS differences between MWA and RTE-based LST values across the seasons. Red and orange box-plots refer to the RMS differences between RTE and SCA-based LST values, and SCA and MWA-based LST values, respectively.

Figure 3 highlights that RTE and SCA have a high level of agreement with each other regardless of the season. MWA is slightly different from RTE and SCA, and in the summer this difference is more evident due to the higher LST dynamics. Since SCA is derived based on RTE [41], and they have the same input parameters, their similarities can also be seen in Figure 3. MWA, on the other hand, uses T_a instead of ${\mathrm{L}}_{\mathsf{\lambda }}^{\uparrow }$ and ${\mathrm{L}}_{\mathsf{\lambda }}^{\downarrow }$ considered by RTE and SCA. Although the median values of all box-plots and seasons are close to zero (0.11–0.35 K for spring, 0.18–0.94 K for summer, and 0.12–0.43 K for autumn), MWA provides clearly different LST values than RTE and SCA in some summer images. In addition, the mean RMS differences (the cross “x” in boxes) (0.11–0.59 K for spring, 0.24–1.91 K for summer, and 0.12–0.64 K for autumn) reveals the higher variations between MWA and the other two methods in summer. Besides, there are two evident outliers over the maximum value in the summer and autumn for MWA-RTE and SCA-MWA.

6.6. Automated LST Extraction Toolbox for Landsat Missions

Different LST retrieval methods and LSE models are not available in packaged RS or GIS software. To overcome these difficulties, some researchers have developed plugins for different software such as ENVI [107], QGIS [108], ArcGIS [109], C++ [110], Python [111], Visual Basic [112] and ERDAS Imagine [113] to extract LST automatically. ERDAS Imagine and ArcGIS software present a visual programming interface which is of vital importance for users without specific knowledge of classical textual programming languages. The ModelBuilder is the visual programming language of ArcGIS software that enables connecting different steps of algorithms to automate the whole process.

GIS models and many remote sensing algorithms require a series of serial tasks called geo-processing workflows. Thanks to the geospatial models, all steps of an algorithm can be connected to each other to automate the whole process. Although ESRI’s ArcGIS is a commercial GIS software, many people and governmental institutions around the world utilize this software since it presents a substantial context for GIS users. Besides, it offers a powerful geospatial model builder (the ModelBuilder) that allows automation and documentation of an algorithm or method. Please see Appendix F for more information about the LST toolbox created based on this study (Supplementary Materials).

7. Discussion

As stated in Section 4, CBEM and NDVI-based LSE models are two appropriate models for Landsat data; CBEM can be used if emissivity measurements are obtained from in-situ campaigns, but it is not practical if the land cover information is not known accurately. Thus, we evaluated the impact of different NDVI-based LSE models, introduced in previous works, on LST retrieval methods over rural areas. It should be noted that the parameters and coefficients in NDVI-based LSE models were used as proposed in the original articles, and they can also be adapted to any study area with field campaigns. Valor and Caselles [82] stated that the error of their methodology ranges from 0.5% (due to the experimental limitations of the field methods) to 2% (considering the case in which there is no information about the study area). Van de Griend and Owe [83] did not conduct any sensitivity/error analysis but indicated that the correlation coefficient between NDVI and thermal emissivity was 0.941 at a 0.01 level of significance. One of the limitations of NDVI-based LSE models is its ineffectiveness in estimating LSE values for water and urban environments [114,115]. In this study, we just investigated the validation of daytime LST images. For nighttime LST evaluation, the LSE image can be estimated using CBEM derived from in-situ measurements or daytime NDVI-based LSE acquired on close dates. Overall, the proposed LST retrieval methods and LSE models can be implemented for regions other than the US, as well as for nighttime, non-rural, and winter data in clear sky conditions.

Considering the LST validation, error sources come from both satellite-based LST and ground-based LST. Satellite-based LST retrieval is still a challenging process due to the great variability of Earth surfaces and the necessary a priori knowledge about several parameters such as the atmosphere, the LSE, the meteorological conditions and the sensor specifications (spectral responses, signal to noise ratio, spectral resolution, spatial resolution, and viewing angle) [7,32,48,116,117,118]. Moreover, LST retrieval methods for satellite data are generally proposed considering different conditions and assumptions. Therefore, there is no universal method that always provides accurate LSTs from all satellite TIR data, and it is not easy to say which algorithm is superior to others [7]. The accuracy of the radiometric measurements and emissivity is the primary uncertainty for ground-based LST retrievals [119,120,121,122,123]. Sobrino and Skoković [119,122] presented an example of an error budget for The Global Change Unit (GCU) sites at the University of Valencia, and

Table 8. A literature example showing the impact of the parameter uncertainty ranges on ground-based LST, from [119,122].

Quantity	Uncertainty	Estimated Impact on Ground-Based LST
Radiometric Calibration	± 0.2 to 0.5 K	0.2 K
Emissivity	± 1%	0.3 K
Downwelling atmospheric radiance	± 10%	0.1 K

indicates the impact of the parameter uncertainty ranges on ground-based LST.

It is interesting to discuss our results in comparison with those of other studies that utilized SURFRAD LST measurements and Landsat data for LST retrieval. Meng et al. [52] estimated LST from Landsat-8 data using the NOAA Joint Polar Satellite System (JPSS) Enterprise algorithm and a hybrid LSE model [82,124]. At the SURFRAD sites, the LST RMSE by the Enterprise algorithm was 3.22 K. Considering our analyses, SWA presented close results to the above analysis under three different LSE models ranging from 2.79 K to 3.02 K. Yu et al. [46] compared RTE, SWA, and SCA methods using Landsat 8 data with their LSE models reported in Section 4.3. They obtained satisfying RMSE values, i.e., 0.9 K, 1.39 K, and 1.03 K for RTE, SCA, and SWA, respectively. However, in this study, using Yu et al.’s LSE model, we obtained the RMSEs as 3.07 K, 3.18 K, and 3.02 K for RTE, SCA, and SWA, respectively. Zhang et al. [57] used Sobrino et al.’s LSE model and SCA method for LST retrieval from Landsat 8 data and compared four Landsat 8 LST images with SURFRAD measurements. Their results showed 1.11 K RMSE with reference to four LST images, whilst we computed 2.94 K RMSE using 15 LST images of Landsat 8 for the same LSE model and retrieval method. In addition to the previous study, Zhang et al. [53] also investigated the accuracy of SCA using Landsat 8 imagery and SURFRAD measurements using 40 Landsat 8 scenes acquired in different seasons and different years, and they obtained 1.96 K RMSE. Wang et al. [54] reported that Practical Single-Channel Algorithm (PSCA) and generalized SCA provided 1.77 K and 2.24 K RMSE, respectively, in line with our SCA results based on Landsat 8 and Sobrino et al.’s LSE model. Sekertekin [59] computed 3.12 K RMSE from 20 Landsat 8 images, close to the RTE-based Landsat 8 LST accuracy found in this test (2.73 K RMSE).

The limitations of this study and future investigations can be reported as follows: (1) Thermal bands have a native spatial resolution of 120 m, 60 m and 100 m for Landsat 5 TM, 7 ETM+, and 8 TIRS, respectively, but they are delivered by USGS at 30-m after cubic convolution resampling. Therefore, we also considered 30 m resampled TIR images which may lead to an error source when conducting pixel-scale validation. Different downscaling methods for TIR or LST data can be utilized as future work to examine the LST accuracy. (2) Although previous work proved that the accuracy of SURFRAD LST varies between 0.1 K to 0.4 K when the broadband emissivity error is about ±0.01, the use of fixed broadband emissivity (0.97) in this study may influence the ground-based LST calculation, but not dominantly. Broadband emissivity estimation models can be implemented in the future to show the optimal model for ground-based LST retrieval. (3) The LST toolbox, presented as a user facility in this study, does not include error analysis. Thus, users should carry out accuracy analysis after obtaining LST images. Besides, it can be improved and generated in an open-source environment as future work.

8. Conclusions

In this study, three LST retrieval algorithms (RTE, SCA, and MWA) were evaluated using Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS data, and additionally, SWA were assessed for Landsat 8 OLI/TIRS data. Since LSE is one of the most important factors affecting the accuracy of LST retrieval methods, the effects of different NDVI-based LSE models on satellite-based LST accuracy were also investigated. Forty-five images acquired in the Spring-Summer-Autumn period over rural areas in the mid-latitude region in the Northern Hemisphere were obtained over five SURFRAD stations and simultaneous in-situ LST data were utilized for accuracy analyses. To get rid of step-by-step calculation for all LST methods and LSE models as well as for time consuming in processing of the images, an enhanced toolbox was generated for automated LST extraction. This toolbox can be utilized by all GIS users to obtain LST in an easy and practical way.

Three NDVI-based LSE models, namely, Sobrino et al.’s, Valor & Caselles’ and Van De Griend & Owe’s LSE model, were considered for Landsat 5 TM and 7 ETM+ data to investigate their effects on LST methods. In addition to these three LSE models, three more NDVI-based LSE models (Skoković et al.’s, Yu et al.’s, and Li and Jiang’s LSE models) were added to the analyses of Landsat 8 based LST methods. That is, the effects of six LSE models on the performance of LST methods from Landsat 8 data were investigated. To sum up, this study only considered NDVI-based LSE models for the evaluation of LST retrieval methods. Two different approaches were considered: 1) Sensor types of Landsat missions (Landsat 5, 7, and 8) were evaluated separately. 2) LST retrieval methods were compared with each other independently of sensor type, i.e., considering all Landsat missions together. In the toolbox, users can decide which LST method and LSE model they can utilize if they are dealing with the use of Landsat data. Furthermore, if they have their own LSE image, the toolbox makes it possible to use any external LSE image.

The obtained results showed that Sobrino et al.’s LSE model provided the best performance to extract LST for all Landsat missions and LST methods. Although all LST retrieval methods with Sobrino et al.’s LSE model presented satisfying and close results when using Landsat 5 TM data, RTE offered the best accuracy (2.35 K RMSE). The same would apply to Landsat 7 ETM+ data, even if MWA presented the best results (2.24 K RMSE). Again, Sobrino et al.’s LSE model provided higher accuracy for all LST retrieval methods from Landsat 8 data, with MWA as the best method (2.52 K RMSE). Considering all Landsat missions, MWA offered slightly better accuracy than RTE and SCA. Concerning the analyses above, it is hard to say which method is globally the best one, since the accuracy of the input parameters largely affects the performance of the methods. In addition, the spatio-temporal and seasonal comparison among LST retrieval methods revealed that RTE and SCA have a high level of agreement with each other regardless of the season. Instead, MWA presented different results than RTE and SCA, especially in summer.

The results indicated that LSE models have a great impact on the accuracy of satellite-based LST retrieval methods. This study also revealed that Sobrino et al.’s LSE model was the most appropriate model for all Landsat missions and LST retrieval methods over rural areas. Moreover, validation of LST retrieval methods with different LSE models over urban areas is a further challenge to be faced that deserves future investigations.