Estimation of Land Surface Temperature in an Agricultural Region of Bangladesh from Landsat 8: Intercomparison of Four Algorithms

The presence of two thermal bands in Landsat 8 brings the opportunity to use either one or both of these bands to retrieve Land Surface Temperature (LST). In order to compare the performances of existing algorithms, we used four methods to retrieve LST from Landsat 8 and made an intercomparison among them. Apart from the direct use of the Radiative Transfer Equation (RTE), Single-Channel Algorithm and two Split-Window Algorithms were used taking an agricultural region in Bangladesh as the study area. The LSTs retrieved in the four methods were validated in two ways: first, an indirect validation against reference LST, which was obtained in the Atmospheric and Topographic CORection (ATCOR) software module; second, cross-validation with Terra MODerate Resolution Imaging Spectroradiometer (MODIS) daily LSTs that were obtained from the Application for Extracting and Exploring Analysis Ready Samples (AρρEEARS) online tool. Due to the absence of LST-monitoring radiosounding instruments surrounding the study area, in situ LSTs were not available; hence, validation of satellite retrieved LSTs against in situ LSTs was not performed. The atmospheric parameters necessary for the RTE-based method, as well as for other methods, were calculated from the National Centers for Environmental Prediction (NCEP) database using an online atmospheric correction calculator with MODerate resolution atmospheric TRANsmission (MODTRAN) codes. Root-mean-squared-error (RMSE) against reference LST, as well as mean bias error against both reference and MODIS daily LSTs, was used to interpret the relative accuracy of LST results. All four methods were found to result in acceptable LST products, leaving atmospheric water vapor content (w) as the important determinant for the precision result. Considering a set of several Landsat 8 images of different dates, Jiménez-Muñoz et al.’s (2014) Split-Window algorithm was found to result in the lowest mean RMSE of 1.19 °C. Du et al.’s (2015) Split-Window algorithm resulted in mean RMSE of 1.50 °C. The RTE-based direct method and the Single-Channel algorithm provided the mean RMSE of 2.47 °C and 4.11 °C, respectively. For Du et al.’s algorithm, the w range of 0.0 to 6.3 g cm−2 was considered, whereas for the other three methods, w values as retrieved from the NCEP database were considered for corresponding images. Land surface emissivity was retrieved through the Normalized Difference Vegetation Index (NDVI)-threshold method. This intercomparison study provides an LST retrieval methodology for Landsat 8 that involves four algorithms. It proves that (i) better LST results can be obtained using both thermal bands of Landsat 8; (ii) the NCEP database can be used to determine atmospheric parameters using the online calculator; (iii) MODIS daily LSTs from AρρEEARS can be used efficiently in cross-validation and intercomparison of Landsat 8 LST algorithms; and (iv) when in situ LST data are not available, the ATCOR-derived LSTs can be used for indirect verification and intercomparison of Landsat 8 LST algorithms.


Introduction
Estimation of Land Surface Temperature (LST) and the study of its changes over time is an important topic of research because, these days, global climate is changing fast. Therefore, retrieval of LST with new technologies has become an interesting field to explore in order to better understand the environment all over the world. With the recent advancement in remote sensing earth observation systems, studying LST, as well as land use and land cover (LULC), has become much easier than it was before. Today, many sources of satellite images are available containing optical, as well as thermal, information of earth surfaces.
LST is the thermodynamic skin temperature of land surfaces which can be studied by measuring the infrared radiation coming from the surface [1]. With LST information, urban heat island can be monitored [2,3] and forest fire can be detected [4]. LST information can be useful to estimate the soil moisture [5][6][7]; hence, studies related to many hydrological processes can be explored from LST [8]. It can also help in different climate studies and weather forecast [8][9][10]. Changes in LST over time can be related with changes in LULC types [11]. LST is related with all sorts of processes that control the energy and water fluxes over the interfaces between the Earth's surface and the atmosphere [12,13]. All these applications make the study of LST a crucial parameter to better understand the regional, as well as the global, environment and its changes over time.
There are different algorithms proposed for the retrieval of LST using different sources of Remote Sensing (RS) data [14][15][16][17][18][19][20]. Among those sources, Landsat has the biggest archive of free images and are of great interest among researches for LST study. Studies with Landsat data include LST retrieval with Landsat 5 TM (Thematic Mapper) data over an agricultural region of Spain by Sobrino et al. [21]; LST with Landsat 7 ETM+ (Enhanced Thematic Mapper Plus) data over Maraqeh County of Iran [22]; and, a work by Fu and Weng [23] for consistent and daily LST monitoring using Landsat images from 1984 to 2011 over the Beijing city of China. Mallick et al. [24] used Landsat 7 ETM+ data to perform LULC and LST study over a heterogeneous urban area of India. Sahana et al. [25] studied LULC change and its impact on LST using Landsat 5 TM data and Landsat 8 thermal infrared (TIR) data over the Sundarban biosphere reserve in India.
Study of LST with RS data over different areas in Bangladesh is rather limited compared to other countries of the world. Bangladesh covers an area of 147,570 km 2 [26] with a huge population. Geographically it extends from 20 • 34 N to 26 • 38 N latitude and from 88 • 1 E to 92 • 41 E longitude [26]. It is one of the most densely populated countries in the world [27]. Most of its population live in its capital city. An agricultural region close to the capital city was selected as the study area for this work (see Section 2.2).
Speaking of the limited LST research works in Bangladesh, those found in literature include the capital city and some of its surrounding areas [28][29][30], and at least one work covering the whole country [31]. Changes in LST and land cover over time, as well as future LST simulation, was studied by Ahmed et al. [28] in the capital city of Bangladesh using Landsat 4, 5 TM, and Landsat 7 ETM+ data. Reja [29] retrieved LST using Landsat 4 TM and Landsat 7 ETM+ data. Sultana et al. [31] used NOAA-16 (National Oceanic and Atmospheric Administration -16) and NOAA-17 data to estimate the minimum and maximum LSTs for six different seasons of whole country. Roni [30] studied the relation between LST and Normalized Difference Vegetation Index (NDVI) using Landsat TM data. Ara et al. [32] used Landsat 8, along with Landsat TM and ETM+ data, to study the effects of land use intensity on LST in Chittagong city corporation area. Study of LST with Landsat 8 images in the capital city or its surrounding areas in Bangladesh is not known as of the time of writing this paper.
Retrieval of LST with the precision result depends on data type, environmental conditions, and the particular algorithm used for the calculation. LST retrieval algorithms may depend on the presence of one or multiple thermal bands in the source RS data. Thermal channel in electromagnetic spectrum covers the region of 10 to 12 µm. Previous Landsat missions came with only one thermal channel but Landsat 8 has two. Sensors with more than one thermal band allow the user to extract information from both of these bands with the possibility of better LST retrieval.

Materials and Methods
Landsat 8 images were used as the primary sources of data to retrieve the LST products. Two Split-Window algorithms, a Single-Channel Algorithm, and an RTE-based direct method were used. The Landsat 8 satellite-retrieved LST products obtained in four methods were validated against reference LSTs and MODIS daily LSTs. The idea is to perform an intercomparison, and study the relative performances of four existing LST algorithms from Landsat 8. The reference LST was retrieved with the Atmospheric and Topographic CORection (ATCOR) module. The MODIS daily LST products were extracted using the AρρEEARS online tool. A study area in an agricultural region of Bangladesh was selected for this research.

Dataset
Landsat 8, the primary source of data in this work, travels on the descending (daytime) node from north to south and crosses the equator at 10:11 a.m. ±15 min mean local time. For our study area, primarily one Landsat 8 image was downloaded from the USGS earth explorer website. The image was taken on 21 February 2018-the actual date on which the authors conducted a field work in the study area-with Operational Land Imager (OLI) and Thermal Infrared (TIR) scanners onboard Landsat 8, and has Path 137 and Row 44. The TIR bands were used to retrieve LSTs, along with the red (band 4) and infrared (band 5) channels for estimating NDVI. The shapefiles necessary for our study were downloaded from the Database of Global Administrative Areas (GADM-https://gadm.org/).
In order to test the LST algorithms on other Landsat 8 images, a set of four Landsat 8 scenes of different dates were downloaded. Two of them were images taken before 21 February 2018, while the other two images were taken after 21 February. We did not use the image of 5 February 2018 because this image was found covered with 77.4% cloud. The amount of cloud cover in percentage is available in the Landsat 8 image metadata file as CLOUD_COVER. The approach used for the cloud classification in Landsat 8 images involve multiple algorithms. It is collectively known as the Cloud Cover Assessment (CCA) and includes the Automated Cloud Cover (ACCA), See-5 CCA, Cirrus CCA, AT-ACCA, etc., algorithms [48] to classify clouds. The CCA analysis results are then merged into the final L1 quality band and the cloud cover amount is made available in the image metadata file. Further details regarding the cloud detection algorithms of Landsat 8 can be found in Reference [48].
To perform cross-validation of Landsat 8 LST products with MODIS daily LST data, we downloaded Terra MODIS LST products using the AρρEEARS online application. We extracted the daily MODIS LSTs of different dates by uploading the shapefiles of our study area in the online module. Details of all Landsat 8 and MODIS images used in this study are presented in Table 1.

Study Area
The study area selected for the verification and intercomparison of four LST algorithms is Chandina sub-district under Cumilla district in Bangladesh. Cumilla district is adjacent to the capital city of Bangladesh. Located in South Asia, Bangladesh is virtually surrounded by India and the Bay of Bengal to the south [49]. It is a low-lying country with huge count of rivers.
Cumilla district is situated in the south eastern part of Bangladesh [50]. It has an area of 3085.17 km 2 , and located in between 23 • 2 N to 24 • 47 N latitudes and in between 92 • 39 E to 91 • 22 E longitudes [51]. There are 17 Upazilas (sub-districts) in this district, among which Chandina is one. Figure 1 shows the location of this sub-district as our study area.
Chandina sub-district is an area of about 201 km 2 , with mostly vegetative land surfaces. A visit to the study area was arranged as part of this research, and it was found that different types of vegetables are grown in this area.
The study area observation was conducted in two consecutive days, 21 and 22 February in 2018, under clear sky conditions, starting the visit in the morning at around 8:30 local time till 15:30 in the afternoon. There were no radiosounding instruments available surrounding the study area; hence, in situ LST data were not monitored. It was found that the local experts usually monitor the LST with soil thermometers in these areas, which may not be a good measurement of in situ LSTs compared to radiometer-retrieved in situ LSTs. Therefore, the intercomparison study of four LST algorithms from Landsat 8 was carried out against ATCOR-derived reference LSTs and AρρEEARS-derived MODIS daily LSTs.

Research Methodology
The methodology to perform the intercomparison of different LST algorithms requires the estimation of LSTs from Landsat 8 images using each of these algorithms. A flowchart representing different steps in LST estimation with four algorithms is shown in Figure 2. As shown in this flowchart, processing of the remote sensing data starts with the level-1 product of Landsat 8 images. Then, the first step is to remove the cloud and haze from level-1 DN values, which is a part of the atmospheric correction. In this study, we used the ATCOR module for atmospheric correction, which uses MODerate resolution atmospheric TRANsmission (MODTRAN) codes to classify and remove cloud and haze from a Landsat 8 scene.
Once the remote sensing images are atmospherically corrected, both optical (OLI) and thermal (TIR) bands are necessary to determine the LST. Using optical bands, band 4 (red) and band 5 (infrared) in particular, NDVI is estimated. From the NDVI, first the proportion of vegetation (P v ), then land surface emissivity (LSE) is determined.
Using both thermal bands of Landsat 8 data, top-of-atmosphere radiance (L ToA ) is first determined. Then, top-of-atmosphere temperature (T ToA ) is estimated. After that, using appropriate LST algorithm with its mathematical formula and its coefficient values (see Section 3), LST can be retrieved using emissivity and T ToA as shown in the flowchart.
Once LSTs with different algorithms are estimated, they can be validated against in situ LSTs or reference LSTs. In this study, the algorithm-retrieved LSTs from Landsat 8 were validated comparing them with reference LSTs retrieved using the ATCOR module. Validation of algorithm-retrieved LST against ATCOR-retrieved LST can be called as indirect verification [53]. The indirect verification was chosen because the in situ LST data were not available due to the absence of radiosounding instruments in the study area. In addition to the indirect verification, cross-validation of Landsat 8 LSTs obtained from different algorithms was performed comparing them with MODIS daily LSTs. Intercomparison study of four Landsat 8 LST algorithms was performed against reference LSTs and MODIS daily LSTs. ToA temperature T ToA 10 and T ToA 11 eq. (33) LST retrieval with SW algorithm eq. (14) or eq. (16) emissivity LUT, spectral response curve emissivity: ε s , ε v (band 10 and 11) Table 5 LSE retrieval LSE TIR 10 , LSE TIR 11 eq. (29), (30) LST retrieval with SC algorithm eq. (4) LST retrieval with RTE eq.

Four Methods to Retrieve Land Surface Temperature from Landsat 8
The thermal channels of Landsat 8 images are band 10 and band 11, also known as TIR 1 and TIR 2 channels, respectively. First of these two has thermal window of 10.60 to 11.19 µm, and the latter 11.50 to 12.51 µm [48]. Since objects or land surfaces transmit radiation in different amounts depending on the wavelength of the channel window, reflected radiation recorded from two thermal bands contain different information about the land surface.
Retrieval methods of LST from thermal bands are based on the Radiative Transfer Equation (RTE). This equation can be used directly to estimate LST; or Single-Channel, Split-Window etc. methods can be used. The RTE-based direct method requires the information of atmospheric parameters for the study area. On the other hand, Single-Channel or Split-Window methods do not require the direct input of those parameters but use some coefficient values obtained through simulation calculations.
While LST can be estimated using one thermal band, the calculation of differential absorption in multiple TIR bands minimizes the atmospheric effects and has the potential to provide better results compared to the use of one thermal band [54]. This idea of using more than one thermal bands is called Split-Window technique [10,14]. Split-Window algorithms for Landsat 8 thermal bands were proposed by many researchers including Rozenstein et al. [55], Jiménez-Muñoz et al. [37], Yu et al. [41], and Du et al. [40]. In this study, we used Split-Window algorithms developed by Jiménez-Muñoz et al. [37] and by Du et al. [40] because: (i) these two algorithms were found providing good performances in different w ranges [40], and (ii) first of these two requires direct input of w value, whereas the second comes with algorithm coefficients for several w sub-ranges. Therefore, it would be interesting to see how these two algorithms perform by making use of w in different ways.
When a Single-Channel or Split-Window algorithm is developed, using the MODerate resolution atmospheric TRANsmission (MODTRAN) or other radiative transfer codes, algorithm coefficients are estimated and made available to users so that they can be input directly in the mathematical expression of the algorithm to study LST for different areas of the world. Among various atmospheric profile databases that can be used in MODTRAN codes, examples include the Thermodynamic Initial Guess Retrieval (TIGR) profile constructed by the Laboratoire de Météorologie Dynamique [40], the Global Atmospheric Profiles from Reanalysis Information (GAPRI) database [37,56], the SAFREE database and the Cloudless Land Atmosphere Radiosounding (CLAR) database [56], the National Centers for Environmental Prediction (NCEP) database [57], etc.
The MODTRAN radiative transfer codes take atmospheric parameters and surface parameters as inputs. Spectral parameters needed for the MODTRAN codes are obtained from the spectral response functions of given TIR sensors [37,40,58]. Surface parameters for the MODTRAN codes may include emissivity [57], viewing geometry, etc. The atmospheric parameters retrieved from these codes include the transmittance (τ), up-welling path radiance (L up ), and down-welling path radiance (L down ) for the given thermal channel. In the following, we provide short description of four LST retrieval methods for Landsat 8 that were used for intercomparison in this study.

Radiative Transfer Equation and Atmospheric Parameters to Retrieve Land Surface Temperature
The formula to retrieve land surface temperature using Radiative Transfer Equation (RTE) can be expressed as [41]: where λ i is the effective wavelength of band i; B i (T i ) is the ToA thermal radiance, τ i is the band average atmospheric transmittance, and ε i is the emissivity of the same band; L up and L down are the upwelling and downwelling radiance in the atmosphere obtained in band i; c 1 and c 2 are Planck's first and second radiation constants, respectively, with c 1 = 1.191 04 × 10 8 W µm 4 m −2 sr −1 and c 2 = 1.438 77 × 10 4 µm K. To retrieve LST with RTE using Equation (1) three atmospheric parameters are needed: τ i , L up , and L down ; besides these parameters, land surface emissivity (ε i ) is necessary. The λ i in Equation (1) (or λ eff ) for a given band can be estimated as [41,42]: where f (λ) is calculated as the function of the spectral responsivity of thermal bands; λ 1 and λ 2 are the lower and upper limit of f (λ). As shown in Figure 3, the λ eff for Landsat 8 band 10 is 10.8 µm and for band 11 is 12 µm. For Landsat 8 TIR band 10, Equation (1) can be rewritten as ( The atmospheric parameters needed to retrieve LST using Equation (3) can be estimated with local radiosounding instruments (if available) or from global atmospheric profiles using simulation codes. The online atmospheric correction calculator used in this study to estimate atmospheric parameters is available at https://atmcorr.gsfc.nasa.gov/ [57,60]. This calculator extracts required parameters from NCEP database using MODTRAN codes and spectral response curve of Landsat 8, Landsat 7, or Landsat 6. Coll et al. [61] have validated this web-based tool against ground measurements for Landsat 7 and reported that the atmospheric correction from this tool is comparable with correction from local radiosonde profiles. The tool was also used to validate a newly proposed pixel-by-pixel atmospheric correction method called SBAC for Landsat 7 in Reference [62]. The online atmospheric correction tool requires the user to input some mandatory data including location (longitude and latitude), date, and time for which the atmospheric parameters are to calculate. Optionally some surface conditions can be input, but if left empty, are assumed from the atmospheric database. The calculated parameters (τ i , L up , and L down ) are sent to the user via email.

Single-Channel Algorithm by Jiménez-Muñoz et al. (2014) with Coefficients
In the Single-Channel algorithm, only one thermal band is used to retrieve LST. For Landsat 8, either band 10 or band 11 can be used. The SC algorithm used in this study was developed by Jiménez-Muñoz and Sobrino [36] and was validated for AVHRR channels 4 and 5, ATSR2 channels 1 and 2, and Landsat TM band 6 data. Later this method was adapted for Landsat 8 by Jiménez-Muñoz et al. [37]. According to them, the SC algorithm for Landsat 8 can be expressed as where ε is the land surface emissivity (LSE); ψ 1 , ψ 2 , and ψ 3 are atmospheric functions (AFs). The symbols δ and γ represent two parameters that can be estimated from linear approximation of Planck functions as [21]: The above equation can be rewritten for simplification as where c 1 , c 2 , and λ holds same meaning as described in Section 3.1. The parameter δ in Equation (4) can be calculated according to Reference [21] as The AFs (ψ 1 , ψ 2 , and ψ 3 ) can be approximated with a second-order polynomial fit against atmospheric column water vapor content (w). If atmospheric function Ψ is considered as a function of water vapor content W as Ψ = CW, the matrix notation to determine the values of AFs can be expressed according to Reference [42] as where c ij are coefficients that were determined by Jiménez-Muñoz et al. [37] for Landsat 8 TIR band 10 using the GAPRI_4838 database as Combining Equations (8) and (9), we can write Using Equations (10)- (12), for a given water vapor content (w) the atmospheric functions (AFs) for Single-Channel algorithm can be calculated and then input in the algorithm expressed in Equation (4) to retrieve LST from Landsat 8. Considering different amounts of water vapor content (w), the values for ψ 1 , ψ 2 , and ψ 3 can be calculated as shown in Table 2.
It should be mentioned here that the Equation (4) gives LST products in K (Kelvin) unit when Equation 32 is used to calculate the ToA brightness temperature (T ToA ). In order to get the LST in • C (degree Celsius) unit, the T ToA needs to be calculated in • C unit using Equation 33, which involves the subtraction of 273.15. This T ToA should then be used to compute the γ and δ parameters. Another way to get the LST products in • C unit is to subtract 273.15 directly in Equation (4).

Split-Window Algorithm by Jiménez-Muñoz et al. (2014) with Coefficients
The Split-Window Algorithm by Jiménez-Muñoz et al. (2014) is based on the mathematical structure by Sobrino et al. [38]; it was later modified by Sobrino and Raissouni [63]. The same mathematical expression was used by Jiménez-Muñoz and Sobrino [58] to retrieve LST from different low-resolution thermal data. Then, it was adapted for Landsat 8 by Jiménez-Muñoz et al. [37]. According to this Split-Window algorithm, LST can be estimated with the following expression: where T i and T j are ToA brightness temperature for band i and band j. For Landsat 8, these two bands are band 10 and band 11, respectively; ε m is the mean land surface emissivity of two bands and is calculated with 1 2 (LSE 10 + LSE 11 ); ∆ε is the difference in LSE of two bands, calculated as ∆LSE = LSE 10 − LSE 11 ; c 0 , c 1 , . . . , c 6 are algorithm coefficients derived from atmospheric simulation codes; w is the total atmospheric water vapor content in g cm −2 unit. For Landsat 8 TIR bands 10 and 11, Equation (13) can be rewritten as The coefficients (c 0 , c 1 , . . . , c 6 ) for this algorithm, specifically for Landsat 8 data, are given by Jiménez-Muñoz et al. [37]. They have used the GAPRI database to estimate the coefficients using MODTRAN simulation codes. Further details regarding GAPRI database can be found in [37]; determination procedure for algorithm coefficients from different atmospheric databases can be found in [42,58].
The seven coefficient values for Jiménez-Muñoz et al.'s Split-Window algorithm are presented in Table 3. The atmospheric water vapor content (w) is not included during the determination of coefficients; it has to be put in the mathematical Expression (14) of the algorithm directly.  [40] is based on the generalized Split-Window algorithm by Wan [39], which was developed for MODerate Resolution Imaging Spectrometer (MODIS) data. The Du et al. algorithm was developed specifically for Landsat 8 data; they named this method as practical Split-Window Algorithm. The mathematical expression of the algorithm is [40]: where b 0 , b 1 , . . . b 7 are algorithm coefficients derived from atmospheric profile dataset using simulation codes. All other parameters in Equation (15) are same as in Equation (13), except the w in Equation (13), which is absent in Equation (15). Equation (15) for Landsat 8 data, with more specific expressions, can be rewritten as T ToA 10 + T ToA 11 2 Other than the algorithm coefficients, parameters in this equation are same as in Equation (14). The coefficients (b 0 , b 1 , . . . b 7 ) for this algorithm were estimated by Du et al. [40] under different atmospheric and surface conditions, through numerical simulation, and using atmospheric profiles. As atmospheric profile they have used the TIGR database; and to perform the numerical simulation they used the MODTRAN 5.2 atmospheric transmittance/radiance codes. They used water vapor content (w) divided into five sub-ranges to calculate the coefficients. Under different atmospheric conditions, users can use the coefficients for their needed sub-ranges. Further details regarding the procedure of algorithm coefficients' estimation can be found in [40]. The coefficients for Du et al. algorithm covering different w sub-ranges are presented in Table 4.  [40]. Using the two Split-Window LST retrieval algorithms in Equations (14) and (16), along with their respective coefficients, LST can be determined once the column water vapor (w) of the study area, ToA brightness temperature (T ToA ) for both thermal bands, and LSE information is available. The T ToA is usually estimated from the thermal bands of remote sensing data while the LSE is determined using optical bands as described in the following sections.

Normalized Difference Vegetation Index and the Proportion of Vegetation Cover
In scientific studies of land resources, vegetation index is often used to express the amount of living plant coverage on lands. The term to mathematically express this indication is called Normalized Difference Vegetation Index, or NDVI for short. It is calculated from the visible red and near-infrared (NIR) light reflected by vegetation and can be expressed as where the subscript "DN" stands for "Digital Number"; it implies that the calculation was made using level-1 DN values stored in remote sensing images. As remote sensing images are subject to cloud covers, atmospheric scattering effects, viewing angle problem, etc., the raw RS data need to be converted into surface reflectance values correcting those effects. In practice, the level-1 DN values are first converted into top-of-atmosphere (ToA) spectral reflectance (ρ λ,ToA ). It is a unitless quantity and can be calculated from OLI data of Landsat 8 as [48]: where M ρ is the reflectance multiplicative scaling factor for the given band, available in the image metadata file as REFLECTANCE_MULT_BAND_n, with n being 1 through 9; A ρ is the reflectance additive scaling factor for the given band, available in the same file as REFLECTANCE_ADD_BAND_n, with n being 1 through 9; Q cal is the level-1 pixel value stored as DN values in the image for both OLI and TIR bands; the θ on left hand side means that this reflectance (ρ λ,ToA(θ) ) does not involve the correction for sun elevation angle.
Making the correction for sun elevation angle, we can exclude the θ from reflectance notation and express the reflectance as [48]: where ρ λ,ToA(θ) is the ToA spectral reflectance without sun elevation angle correction as calculated in Equation (18); θ is solar elevation angle, which is the local sun elevation angle at the time of satellite overpass available in the image metadata file of Landsat 8 as SUN_ELEVATION. This value needs to be checked for each Landsat 8 scene. The sun elevation angle can be related with solar zenith angle (θ z ) using cosine function as where θ z is the solar zenith angle with θ z = 90 • − θ. Combining Equations (18) and (19), we can calculate ToA spectral reflectance as To determine NDVI, we need the reflectance values for red and near-infrared bands; this translates to band 4 and band 5 of Landsat 8 image, respectively. The red band in Landsat 8 covers 0.636 to 0.673 µm and NIR band 0.851 to 0.879 µm in the electromagnetic spectrum [48]. For band 4 and band 5, M ρ value is 2.0000 × 10 −5 , and A ρ value is −0.100000 [48], as given in the image metadata file. Once the reflectance is corrected for the sun elevation angle, it needs another correction-the atmospheric correction because the ToA spectral reflectance includes atmospheric scattering effects and therefore does not represent the true reflectance of land surfaces. In order to get the reflectance from the surface (ρ λ,LS ), we need to correct atmospheric effects including the cloud cover and atmospheric gases. This correction can be made using some atmospheric correction modules available in various image processing software. Then, surface reflectance can be determined using Equation (21).
To simplify the notation of surface reflectance, in the following we used ρ in place of ρ λ,LS . According to this notation, the NDVI of land surface (NDVI LS ) can be expressed as For Landsat 8 image data, Equation (22) takes the form: Using NDVI LS values as in Equation (23), the proportion of vegetation cover can be calculated as [64]: where NDVI max = 0.5 indicates the presence of vegetation on lands, and NDVI min = 0.2 represents only bare soil on land surfaces. NDVI value less than 0 indicates the water, and NDVI value greater than 0.5 indicates full vegetation [65]. When NDVI value ranges between 0.2 and 0.5, the surface is considered as a mixture of soil and vegetation, requiring the calculation of P v using Equation (24). Using the P v , the emissivity of the mixed land surface is then estimated as described in the following Section.

Land Surface Emissivity Determination
Quantitatively, emissivity is the ratio of the thermal radiation from a surface to the radiation from an ideal black surface at the same temperature as given by the Stefan-Boltzmann law [1]. To retrieve LST from ToA brightness temperature (see Section 3.7), we must consider emissivity from land surfaces. The term Land Surface Emissivity (LSE) is used to indicate emissivity from land surfaces that are composed of different types of materials (soils, vegetation, water etc.). One way to get LSE is through normalized difference vegetation index since NDVI represents the greenness of land surfaces, giving an idea of types of materials comprising the land surfaces. Thus, different values of NDVI represent different materials of land surfaces. If the NDVI value is less than 0.2, then it represents only bare soil. In this case, the emissivity can be calculated from reflectivity values in red region of the image. If the NDVI value is greater than 0.5, then the land surface is composed of vegetation only. If that is the case in a real study, a constant value of emissivity, typically 0.99, can be used [21]. But in situations when NDVI lies between 0.2 and 0.5, LSE for a given band i can be related with NDVI and proportion of vegetation (P v ) as [41]: where ε v,i is the emissivity of fully vegetated surfaces and ε s,i is emissivity of barren soil, in the band i; P v is the proportion of vegetation as calculated in Equation (24); a i , b i are the coefficients that can be estimated from laboratory spectra of soils using statistical fits, assuming that the emissivity and the reflectivities in red band have a linear relationship [66]. The symbol C i in the above equation denotes the roughness of land surfaces [65]. For plain and homogeneous land surfaces, this C i can be neglected [21], and considered C i = 0 [65]. For rough and heterogeneous surfaces, i.e., soil-vegetation mixed pixels, C i denotes the increment in emissivity resulted from the cavity effect and multiple scattering in the mixed pixels [67].
Taking emissivity values of soil and vegetation into account, and assuming that NDVI values of earth surfaces range from around 0.2 to 0.5, the emissivity of land surfaces (LSE) can be calculated according to the NDVI-threshold method as [41,66]: The C in the above equation is the same as in Equation (25), which can be found also as dε in literature [21]. According to Equation (26), to calculate LSE for Landsat 8 thermal bands, we need the ε s and ε v values estimated for both TIR bands. Yu et al. [41] estimated these values using the MODIS UCSB (University of California, Santa Barbara, CA, USA) emissivity library (https://icess.eri.ucsb.edu/ modis/EMIS/html/em.html), as presented in Table 5. In Equation (26), an approximation of C is given by [41,68]: where F is a shape factor [69]. Sobrino et al. [21] considered this shape factor (F) under different geometrical distributions having a mean value of 0.55.
Taking Equations (26) and (27) into account, the LSE can be calculated as with Based on the emissivity values in Table 5 and the mathematical expression in Equation (28), the LSE for both TIR bands of Landsat 8 data can be calculated. To do so, we first need to calculate the m and n values for these bands. For TIR band 10, we get Thus, using Equation (28), we get the LSE for TIR band 10 as Similarly, for TIR band 11, we get Therefore, LSE for TIR band 11 can be expressed as Using the proportion of vegetation cover as described in Section 3.5 and expressed in Equation (24), LSE for TIR band 10 and band 11 can be calculated using Equations (29) and (30), respectively.

Top-of-Atmosphere (ToA) Brightness Temperature Determination
Determination of ToA brightness temperature (T ToA ) can be described as a two-step process. The first step includes the conversion of level-1 DN values of Landsat 8 thermal infrared data to at-satellite (or at-sensor, or ToA) spectral radiance values. This is the spectral radiance in wavelength of a surface, which is expressed with L λ,ToA and has a unit of watt per meter squared per steradian per micro meter (W m −2 sr −1 µm −1 ). The formula to convert level-1 DN values in RS images to spectral radiance is [48]: where M L is the radiance multiplicative scaling factor for the given band, available in the image metadata file as RADIANCE_MULT_BAND_n, with n being 1 through 11; A L is the radiance additive scaling factor for the given band, available in the same file as RADIANCE_ADD_BAND_n, with n being 1 through 11; Q cal is the level-1 pixel value stored as DN values in the image, available for both OLI and TIRS bands. The M L value for both TIR bands (band 10 and band 11) is 3.3420 × 10 −4 ; and A L value for both of these bands is 0.10000. Spectral radiance can be obtained for both OLI and TIR bands but in order to retrieve LST from Landsat 8 thermal bands, radiance from only thermal bands is necessary. For both TIR bands, L ToA 10 and L ToA 11 can be estimated following Equation (31). The second step to determine T ToA involves the use of L λ,ToA image data from Equation (31). Then, T ToA in K (Kelvin) unit can be calculated by inverting Planck's radiation equation as [48]: where K 1 and K 2 are the thermal conversion constants for the given band, available in the image metadata file as K1_CONSTANT_BAND_n and K2_CONSTANT_BAND_n, respectively, with n being 10 or 11; L λ,ToA is the ToA spectral radiance calculated for band 10 or band 11 with Equation (31 According to Equation (33), the T ToA 10 and T ToA 11 for both TIR bands of Landsat 8 can be estimated in • C unit.

Preparation and Processing of Data for LST Retrieval from Landsat 8
In order to perform atmospheric corrections of Landsat 8 data, as well as for other processing of raster images, different computer tools were used in this study. First, appropriate Landsat 8 images were downloaded from the website of USGS earth explorer. Vector shapefiles were downloaded from the GADM website (http://gadm.org/data.html). Spatial subset of the raster image was created using the shapefile for our study area in QGIS [70] software. Atmospheric corrections including the removal of cloud and haze from the level-1 Landsat 8 data were done using the ATCOR module in PCI Geomatica 2016 software. From the atmospherically corrected subset image data, NDVI was determined using raster calculator of QGIS. Other calculations, for example, proportion of vegetation cover, LSE, LST, etc. were also estimated with the same tool. The LST maps were exported using the print composer function of QGIS. Statistical analyses and all types of plots and histograms were created using the R program [71] with substantial help from the raster [72], rgdal [73], ggplot2 [74], and caret [75] libraries.

Results and Discussion
The LST products retrieved from all four algorithms were validated against reference LST (LST ref ) and cross-validated against MODIS daily LST (LST MOD ). The RMSE (root-mean-squared-error) was used to measure the accuracy of LST retrieved with four methods comparing them against reference LST and MODIS LST. The reference LSTs for all Landsat 8 images were estimated using the ATCOR module available in PCI Geomatica 2016 software. The MODIS daily LST was retrieved using the AρρEEARS online application from the Terra MODIS images. The online atmospheric correction calculator (https://atmcorr.gsfc.nasa.gov/) was used for the estimation of atmospheric parameters, including the water vapor content. For several locations in the study area, a total of 48 calculations of w for each Landsat 8 scene were made using the online tool and the w range was found to be 0.6 to 3.0 g cm −2 considering all five Landsat 8 images. Therefore, we have used several values of w for LST retrieval.
The LST result retrieved with the RTE-based direct method is denoted as LST RTE , LST with Single-Channel algorithm as LST SC

Results from RTE Method using TIR Band 10
The atmospheric parameters required to retrieve LST in direct method using RTE (3) were obtained using the online atmospheric correction calculator for six specific locations of the study area as presented in Table 6. The values in several observations were found to be close to each other.
In order to compute the RMSE, we first calculated ∆LST comparing the RTE-retrieved LST and the reference LST. The validation results are represented with histograms in Figure 4 including the RMSE values. The smallest RMSE for LST RTE was found to be 1.95 • C, with the largest being 2.67 • C.

Results from Single-Channel Algorithm using TIR Band 10
The LST retrieved with Single-Channel algorithm from band 10 of Landsat 8 uses Equation (4) and atmospheric functions calculated for different values of water vapor content as given in Table 2. For our study area, we retrieved LSTs for w range of 0.5 to 3.0 g cm −2 with an interval of 0.5 g cm −2 . This gives LST calculations for six values of w. All LST products obtained in this method were then validated against reference LST as shown in Figure 5 with their RMSEs. The smallest RMSE was found to be 4.16 • C for w = 3.0 g cm −2 , whereas the largest RMSE was 4.17 • C for w = 0.5 g cm −2 .

Results from Jiménez-Muñoz et al.'s Split-Window Algorithm
The coefficient values in Jiménez-Muñoz et al.'s (2014) algorithm (see Table 3) do not include the water vapor content (w); instead, the algorithm in Equation (14) takes the w as a direct input from the user. Like in Single-Channel algorithm, we used six values of w to retrieve LSTs in this method (LST Jim ). To compute RMSEs, all observations of LST Jim were validated against reference LST and presented as histograms in Figure 6. The smallest RMSE for LST Jim was found to be 0.74 • C when w = 0.5 g cm −2 , with the largest RMSE being 0.94 • C when w = 3.0 g cm −2 .

Results from Du et al.'s Split-Window Algorithm
As it is mentioned previously, the coefficient values in Du et al.'s (2015) algorithm were estimated using water vapor contents (w) of various sub-ranges (see Table 4). Therefore, LSTs were estimated with this algorithm (LST Du ) using algorithm coefficients given for w sub-ranges of 0.0 to 2.5 g cm −2 , 2.0 to 3.5 g cm −2 , 3.0 to 4.5 g cm −2 , 4.0 to 5.5 g cm −2 , 5.0 to 6.3 g cm −2 , and 0.0 to 6.3 g cm −2 . All six LSTs were then validated against reference LST and RMSEs were computed (Figure 7). The smallest RMSE (1.13 • C) was found when LST was estimated using w in 0.0 to 6.3 g cm −2 while in w range of 2.0 to 3.5 g cm −2 the LST was found with abnormally the largest RMSE (11.05 • C). This error could be due to the wrong use of w amount (see Section 4.6).

Intercomparison of Four Algorithms to Retrieve LST from Landsat 8 Against Reference LSTs
The LST results as estimated from the Landsat 8 image of 21 February 2018 with all four algorithms (see Sections 4.1-4.4) were compared among them validating the results against reference LST by means of RMSE computation and R 2 (coefficient of determination). In addition, algorithm-retrieved LSTs were also observed against reference LST by means of mean bias error (Bias ref ).
Except for LST RTE , a sub-range of w was considered for the average LST values for three other methods (LST SC , LST Jim , and LST Du ). Considering the water vapor content of 1.68 to 2.43 g cm −2 for the 21 February 2018 Landsat 8 image, as retrieved from the NCEP database with MODTRAN codes, a w range of 1.5 to 3.0 g cm −2 was considered as the possible mean of w for LST SC and LST Jim methods. For the LST Du method, the algorithm coefficients available in w range of 0.0 to 6.3 g cm −2 , that is, the entire range of w was considered.
The minimum and maximum LST, along with the statistical mean and standard deviation of LST retrieved with four methods, and the reference LST are presented in Table 7. The RMSE and R 2 , as well as the mean bias (Bias ref ) of algorithm-retrieved LSTs computed against reference LST, are also presented. As seen in Table 7, the mean RMSE was found 2.20 • C for LST RTE ; 4.17 • C for LST SC ; 0.88 • C for LST Jim ; and 1.13 • C for LST Du . The average correlation of coefficient (R 2 ) for LST Jim and LST Du is 0.86; it is 0.77 for LST RTE , and 0.76 for LST SC method. This indication implies that all four algorithms perform efficiently compared to reference LST. The mean bias against reference LST (bias ref ) is the smallest in LST Jim (0.82 • C) and largest in LST SC (4.14 • C), with LST Du (1.08 • C) and LST RTE (2.30 • C) being in between. Considering all these observations, the best performing LST algorithm in our study is LST Jim , with the other three methods staying in close agreement.
The box plots in Figure 8 represent LST results of four algorithms in terms of RMSE, giving a visual aid for intercomparison. As seen in this figure, variation in RMSE is lowest for LST Jim , with LST RTE and LST SC being next to it and LST Du having the highest variation. The rather different RMSE variation in LST Du is probably due to the wrong use of w, which is described in Section 4.6.
The LST maps created using all four methods and with ATCOR module (LST ref ) are shown in Figure 9. We presented the maps with LST results that were found very close to the reference LST, considering the possible amount of w present in our study area. This makes the LST maps for LST SC and LST Jim computed with w = 2.0 g cm −2 , and LST Du computed with w range of 0.0 to 6.3 g cm −2 .   Considering all five Landsat 8 images used in this intercomparison study, the average RMSE is found to be 2.47 • C for the LST RTE method; 4.11 • C for LST SC ; 1.19 • C for LST Jim ; and 1.50 • C for the LST Du algorithm (not shown in Table). It can be mentioned here that the w values retrieved from the NCEP database were different for Landsat 8 images of different dates (see Table 8). The corresponding w ranges were taken into consideration for the calculation of LSTs with different methods.

Variation in LST Results due to Wrong Amount of Water Vapor Contents
The Single-Channel algorithm and the two Split-Window algorithms used to retrieve LST from Landsat 8 are dependent on atmospheric water vapor content (w). Therefore, it is necessary to study their performance in different amounts of w. To perform this study, we retrieved LSTs for the Single-Channel algorithm and Jiménez-Muñoz et al.'s Split-Window algorithm for w amount up to 4.5 g cm −2 , that is, three more observations than it is seen in Table 7. For Du et al.'s Split-Window algorithm, we used the LSTs retrieved for all w sub-ranges.
The variation in RMSEs of LST results for three algorithms was plotted against varied amount of w as shown in Figure 11. Since the algorithms for LST SC and LST Jim require direct input of w, the two plots (Figure 11a,b) share w values in same intervals. On the other hand, RMSE variation for the LST Du algorithm is shown with several sub-ranges of w (Figure 11c). As seen in Figure 11a, for LST SC the RMSE decreases with increasing w, whereas in Figure 11b, RMSE for LST Jim increases with increasing w. The decrease of RMSE for LST SC in lower w is quite abrupt compared to its change in higher w (Figure 11a). For the LST Jim , the increase in RMSE with increasing w is almost constant (Figure 11b). On the other hand, for LST Du (Figure 11c), LST error is very high when w = 0.0 to 4.5 g cm −2 ; it is lowest when w lies in 4.0 to 6.3 g cm −2 ; it is considerably better when w = 0.0 to 6.3 g cm −2 , that is, algorithm coefficients for the entire w range. It is understandable from these plots that the use of wrong w value may result in unacceptable LST estimation, especially for the LST Du algorithm. Therefore, it is very important to estimate the w of the study area with great precision.

Cross-Validation and Intercomparison of Landsat 8 LSTs Against MODIS Daily LSTs
As mentioned previously, MODIS daily LSTs for the study area were extracted using the AρρEEARS online tool. We considered the nearest available LSTs compared with Landsat 8 images of different dates used in this study.  Figure 12 using box plots.
For cross-validation of Landsat 8 LSTs that were retrieved with four algorithms, as well as with the ATCOR module, the MODIS daily LSTs were converted into • C unit from K unit. Then, LST mean bias was computed for each Landsat 8 LST retrieval method against MODIS mean LST. The cross-validation and intercomparison results between Landsat 8 and MODIS daily LSTs are presented in Table 8. The Landsat 8 mean LST is denoted with LST L8 , MODIS mean LST with LST MOD , and the mean LST bias between them is denoted with Bias MOD-L8 .
As seen in Table 8, the mean bias (Bias MOD-L8 ) between LST L8 and LST MOD is always lower for the LST Jim algorithm among four Landsat 8 LST algorithms. The LST Du algorithm performs very close to the LST Jim method. The LST RTE shows higher LST bias while the LST SC has the highest bias among these four methods. It suggests that the best performing LST retrieval method for Landsat 8 is the LST Jim Split-Window algorithm, with LST Du method being the next, and LST RTE having better performance than the LST SC method. These cross-validation and intercomparison results agree with the intercomparison results obtained for different Landsat 8 LST methods when validated against LST ref (see Section 4.5). It also reveals that the reference LST estimated from ATCOR module (LST ref ) performs efficiently showing a very small mean bias (from −0.58 to −0.29 • C) for most images and the lowest mean bias for the first three Landsat 8 images (

Conclusions
Landsat 8 data are great sources of high resolution remote sensing images with two thermal bands that can be efficiently used to retrieve LST. An intercomparison study among the existing LST algorithms for Landsat 8 was performed against ATCOR-derived reference LSTs and AρρEEARS-derived Terra MODIS daily LSTs. According to the observation of this study, the Single-Channel algorithm can be used in LST retrieval for Landsat 8 images, but use of a good Split-Window algorithm has the potential of ensuring greater accuracy. The challenges with a good, practical, and feasible Split-Window algorithm development can be the precise estimation of coefficient values and determination of atmospheric water vapor content of the study area.
Taking all Landsat 8 images used in this study under consideration against reference LST, the RTE-based method (LST RTE ) gives LST results better than the Single-Channel method (LST SC ) with an average RMSE = 2.47 • C; but it (LST RTE ) performs worse compared to the Split-Window algorithms. Since the RTE-based direct algorithm depends heavily on various atmospheric parameters, precision calculation of those parameters can be the determinant of LST accuracy. The use of NCEP database and MODTRAN codes through the online atmospheric correction calculator seems promising for this method.
The LST results with Single-Channel algorithm provide larger RMSE (average RMSE = 4.11 • C) than the RTE-based method. In contrast, Jiménez-Muñoz et al.'s Split-Window Algorithm (LST Jim ) shows the best performance (average RMSE = 1.19 • C) among other methods. These two LST algorithms (LST SC and LST Jim ) can be chosen when: (a) actual water vapor content is precisely measured, and (b) other atmospheric parameters are not available as they are not necessary in these two methods. Especially for the Single-Channel algorithm, it is not advisable to use this method for images that have more than one thermal band.
The Du et al.'s Split-Window Algorithm (LST Du ) with its coefficients available for w range 0.0 to 6.3 g cm −2 can be applied for an area where the actual amount of water vapor content cannot be determined or very uncertain for precision determination. The reasons behind this recommendation are: (a) this algorithm gives good LST results compared to other methods with RMSE = 1.50 • C as found in this study, and (b) other three LST algorithms require either the direct input of water vapor content (LST SC and LST Jim methods) or several atmospheric parameters (LST RTE method). On the other hand, the potential problem with this method is that the w range of 0.0 to 6.3 g cm −2 is too much of generalization and could give unsatisfactory results in areas with different atmospheric conditions than observed in this study.
The cross-validation and intercomparison results of Landsat 8 LSTs with different algorithms against MODIS daily LSTs were found to agree with the intercomparison results against reference LSTs. The mean bias (Bias MOD-L8 ) here for the LST Jim algorithm was found always lower compared to other Landsat 8 LST algorithms. The ATCOR-derived Landsat 8 LST was found with even lower Bias MOD-L8 for the first three images used in this study revealing that the ATCOR-derived LSTs can be used as references for the indirect verification of Landsat 8 LST algorithms.
The in situ LSTs were not available in this study; therefore, ground validation of Landsat 8 LST algorithms was not performed. Monitoring of in situ LST data with precision radiosounding instruments or radiometers correcting for the effects of emissivity and synchronizing with the actual time of satellite overpass could be taken into consideration to perform the ground validation of LST algorithms.
Although the NCEP atmospheric profile database was found providing with good estimation of atmospheric parameters, other databases (e.g., TIGR, GAPRI, CLAR, etc.) can be used to study their relative performances. Land surface emissivity from MODIS or ASTER remote sensing images can be compared against NDVI-based emissivity in the retrieval of precision LST. Cross-validation of LST from Landsat 8 with MODIS daily LST retrieved using the AρρEEARS online tool was found promising in the intercomparison study of different Landsat 8 LST algorithms. Other sources of remote sensing data can be used in the cross-validation study to further enhance the verification of Landsat 8 LST products.