Spatial and Multi-Temporal Analysis of Land Surface Temperature through Landsat 8 Images: Comparison of Algorithms in a Highly Polluted City (Granada)

: Over the past decade, satellite imaging has become a habitual way to determine the land surface temperature (LST). One means entails the use of Landsat 8 images, for which mono window (MW), single channel (SC) and split window (SW) algorithms are needed. Knowing the precision and seasonal variability of the LST can improve urban climate alteration studies, which ultimately help make sustainable decisions in terms of the greater resilience of cities. In this study we determine the LST of a mid-sized city, Granada (Spain), applying six Landsat 8 algorithms that are validated using ambient temperatures. In addition to having a unique geographical location, this city has high pollution and high daily temperature variations, so that it is a very appropriate site for study. Altogether, 11 images with very low cloudiness were taken into account, distributed between November 2019 and October 2020. After data validation by means of R 2 statistical analysis, the root mean square error (RMSE), mean bias error (MBE) and standard deviation (SD) were determined to obtain the coefficients of correlation. Panel data analysis is presented as a novel element with respect to the methods usually used. Results reveal that the SC algorithms prove more effective and reliable in determining the LST of the city studied here.


Introduction
The land surface temperature (LST) is an essential factor when studying physical processes on earth [1]. It is used for reference worldwide, in research dealing with local or regional space and in the fields of hydrology, meteorology, and superficial energy balance [2][3][4], as well as studies surrounding climatic change [5] and recent work involving the phenomenon known as urban heat island (UHI) [6][7][8][9][10]. It can be of vital importance to obtain a reliable LST, and this value could bear an impact upon the quality of life of populations. It has therefore become an essential factor, used to evaluate the interchange of water and superficial energy with the atmosphere [11].
Due to the important economic cost of measurements in situ and the difficulties involved in running them on vast extensions of territory, since the 1990s retrieval of remote sensing satellite images with sensors of thermal infrareds (TIRS) has become a more common approach. Training for the implementation of these systems is an important line of research [12] evidenced by the ample literature on the subject [12][13][14][15][16][17][18][19][20]. Still, the technological improvements incorporated by the infrared sensors of today´s satellites with respect to their predecessors urge researchers to keep on conducting research aimed at further advances in the field.
Numerous LST products have been developed based on satellite data. One of them, the Landsat, was the eighth satellite put into orbit in the year 2013. It features thermal infrared sensing (TIRS) bands 10 and 11, providing LST estimates of 16 days, and a 100-m resolution. When put into operation, the bands gave radiometric errors that surpassed the maximum value permitted (2%)-absolute values of 4%-5% were obtained for the measurements of band 10, and 8%-9% for band 11 [13]. To correct these deficiencies, authors Gerace and Montanaro (2017) [21] proposed an optical flare correction algorithm named the Stray Light Correction Algorithm (SLCA). It came into use in February 2017 for the images of the Terrestrial System of the US Geological Service (USGS) and is used for training in the field through a process of intervening validation [9,10,[13][14][15].
Since the onset of its use, Landsat 8 has proved useful for determining the LST in large cities as well as prosperous, growing mid-size cities such as Lyon (France) [22], Krakow (Poland) [23] and Vienna (Austria) [24]. A variety of validated algorithms may be adopted to recover the LST through Landsat 8 satellite images [16][17][18][19][20]25,26]. They are based on the principles of fluid dynamics and atmospheric physics [27] and comprise three major groups: the mono window (MW) developed by Weng et al. (2004) and Qing et al. (2001) [17] and later improved by   [28]; single channel (SC), developed by Jiménez and Sobrino et al. (2014) [16]; and, finally, the split window (SW) algorithms developed by Jimenez and Sobrino et al. (2014) [16], Du et al. (2015) [20] and Mao et al. (2005) [18]. The two first groups used a single TIRS channel for calculation, typically band 10; the latter group relied on two, bands 10 and 11. MW algorithms are based on individual atmospheric parameters: ambient temperature, humidity, and atmospheric mean temperature [12]. In turn, the SC algorithm [16] accounts for ground emissivity (ε) and steam-driven content of water (w) in the atmosphere (Yu et al. 2014). The great advantage of SW algorithms is that they do not require atmospheric data of great precision [29]. Thus, authors Jiménez et al. (2014) [16] established that they work well in global conditions and provide better results than SC algorithms. Xiangchen et al. (2019) [14] hold that SW algorithms are best suited for Landsat 8 estimates, though further investigation is necessary to validate the precision of these algorithms.
Currently, the use of algorithms to retrieve the LST presents two major challenges. 1. The spatial resolution of the pixels of LST systems can vary from 10 −2 to 10 Km 2 depending on the satellite used. This circumstance can give rise to pixels that represent various types of surfaces, which makes it difficult to obtain the LST [30]. 2. The satellite LST requires a verification process to validate the performance of the algorithms used. Common LST algorithm validation techniques are: in situ measurements, the radiance method, or cross reference and comparison with ambient temperatures near the ground. The first method is difficult and expensive, requiring radiometers located on terrains of homogeneous coverage (lakes, grasslands, dunes and regions with vegetation) and a dimension equal to the pixel of the satellite image. This method can give differences of 2 to 5 Kelvin degrees (K) with respect to actual measurements, according to previous studies [31]. For this reason, only a few validation studies using this method are found in the literature [32][33][34][35]. The main drawback of the second method is that it calls for knowing the emissivity of the soil and atmospheric profiles in situ. These can be simulated using MODTRAN 4.0 software, whose use is common in the literature [13,36,37] but requires verification that the simulated profiles are adapted to the actual atmosphere of the point studied [38]. The third method involves the validation of LST values retrieved against alternative LST values, either well documented or validated by other satellites [19,39]. It has the advantage of not requiring in situ or simulated measurements, but the precision is very sensitive to spatial and temporal mismatches of the measurements [30], and similar spatial resolutions are needed to proceed with the validation. In recent decades, the method of comparison with ambient temperatures near the ground (1-2 m) has gained importance as a validation system for LST algorithms [40][41][42][43][44]. It consists of comparing the recovered LST with the ambient temperature obtained from meteorological stations or temperature probes located near the ground. It is a fast, simple and low-cost method, but its main drawback lies in finding pure pixels of vegetation within cities that do not comprise several surfaces [30,37]. This inconvenience is minimized by using Landsat 8, since the resolution of its thermal bands is 100 m-an adequate distance to locate green or bare areas within or adjacent to cities. The precision of this method has reported differences ranging from 0.2 to 7.8 K between the recovered LST and the measured ambient temperatures [40][41][42][43][44][45].
The city under study presents important problems related to the local climate, conditioned by its proximity to the Sierra Nevada mountain range (30 km away), with average heights of 2045 m above sea level, and to the Mediterranean coast (60 km away). Granada is located in a former basin of the Mediterranean Sea and is considered highly vulnerable to climate change, given its high rate of increasing temperature as compared to the rest of the planet [46,47]. In addition to these geographic and climatological elements, Granada faces environmental challenges, having become the third most polluted city in Spain [48]. To the best of our knowledge, ours is the first study dedicated to LST estimates using satellite images of this area, and one of the first conducted on a Spanish city, after research on Barcelona [49,50] and Madrid [51].
Numerous studies have come to positively link environmental pollution with an increase in LST in urban areas [30,[52][53][54]. Paradoxically, although many LST studies involving Landsat 8 images have been validated by the method of comparison with ambient temperatures [40][41][42][43][44][45], the possible effect of contamination on the environment has not been assessed in terms of the precision and subsequent validation of the LST. In sum, the geographic, climatic and high pollution conditions of the city of Granada make it an ideal place for this research.
The main statistical approaches and methods commonly used in studies of LST recovered through Landsat 8 images are: descriptive statistics of linear correlations and regressions [35,[55][56][57], multiple regression analysis (MRA) [58,59], bivariate correlation [60] and ANOVA [61][62][63]. The panel data method used in this research allows for the incorporation of a greater number of data and variables, admitting the inclusion of individual effects of a certain variable to derive global results. It also allows inclusion of the spatial residual values of the results and eliminates the usual problem of collinearity between variables. Such considerations are often neglected by traditional methods, impeding the arrival at more accurate and complete results.
Therefore, the questions that we aim to answer through this research, in view of the particular characteristics of the city of Granada, are the following: (1) How does one choose the right LST algorithm for Landsat 8? (2) Is the validation method of comparison with ambient temperatures adequate? (3) Does panel data statistical analysis lead to different data than traditional analysis methods?
The foremost objective of this research is to determine the precision and maturity of six LST algorithms on Landsat 8 images validated using the technique of comparing ambient temperatures over the city of Granada (Spain), with its unique geographic, environmental and pollution characteristics. Specific goals include: 1. Determining if high contamination can affect the validation process of the chosen LST. 2. Providing a new lowcost, fast and accurate approach that allows the variations of the LST of cities to be monitored in an urgent and precise way, taking into account physical, environmental and climatic variables. In the long term, such efforts may contribute to decision-making by urban planners and public administrations about the distribution of the LST inside cities.

Area of Study and Data Source
The area studied is the city of Granada, in southern Spain. The geographic coordinates (UTM) of the metropolis are: latitude 37.111741º N and longitude 03.362401º W, and its altitude is 680 m above sea level ( Figure 1). The local weather is strongly conditioned by its location at the feet of the Sierra Nevada, a major Andalusian mountain range with an average height of 2045 m a.s.l. After the Alps, it is the highest range in western Europe, The peak of Mulhacén is 3482 m a.s.l.. According to the Koppen Geiger climatic classification, it features transitional weather between the Mediterranean climate (Bsk) and a cold semiarid climate of moist, mild winters and warm, dry summers [64]. The average yearly temperature fluctuates between 279.65 K in January and 298.45 K in the month of July, with a winter average of 270.15 K and summertime extremes of 316.15 K. The number of hours of sunshine per year is 2,917 hours, giving a mean of 7.99 hours of sun a day. The population is 232,462, although it reaches 540,000 with the inclusion of Granada´s metropolitan area, which occupies a total surface of 880,000 km 2 .

Landsat 8 Images
The satellite Landsat 8 was put into orbit in February 2013 and features two novel sweeping instruments: thermal infrared sensing (TIRS) and the Operational Land Imager (OLI). Images consist of eight spectral bands with 30 m resolution (bands 1 to 7 and 9), one band with 15 m resolution (band 8) and two thermal bands that allow for determining the LST (bands 10 and 11) with a resolution of 100 m [10,6]. These two TIRS bands work with a wavelength interval from 10.6 to 11.19 µm (band 10) and from 11.50 to 12.51 µm (band 11).
The city of Granada lies under the path of satellite Landsat 8, road 20, route 34. The images selected for this study are temporally located between November 2018 and October 2019. In the course of this span, 11 images with an index of cloudiness below 15 were used to enhance accuracy in obtaining the LST object. The images, indexed in Table 1, were acquired from the Geological Service of the United States (USGS). Meteorological conditions (solar radiation, wind, precipitation, etc.) can alter the measurements of the satellite spectrometer. The days selected for obtaining the images did not present precipitation in the previous 96 hours, and the wind speed remained at constant values between 0.83 m/s and 2.22 m/s. The atmospheric parameters (L ↑ λ, L ↓ λ and wind speed) were obtained from Spain´s National Meteorological Agency (AEMET). The solar zenith angle θse was obtained from Landsat metadata files and the θsz values are equal to θse−90º.
After downloading the images, the spatial resolutions of bands 10 and 11 were reclassified to a resolution of 30 m, equal to the resolution of the rest of the Landsat 8 bands. This process was carried out using the raster calculator option of QGIS software, version 3.10.5, providing homogeneous and compatible pixels for the LST and emissivity. The processed images were georeferenced using the projection system ETRS89/UTM Zone 30N. Performing atmospheric correction in bands OLI entailed use of the algorithm Dark Object Subtraction (DOS) [12] by means of the Classification Plugin (SCP) with open source software QGIS [10,[63][64][65][66][67][68].

Measurements of Ambient Temperatures
Spain´s National Meteorological Agency (AEMET) has two meteorology stations in Granada with instrumentation for measuring environmental parameters, calibrated according to international standards. The first is located in a rural area 12 km from the city (UTM: 37.11.45 N; 3.47.22 W), while the second is within the city (UTM: 37.11.23 N; 3.35.44 W). In order to minimize the impact that the urban environment and pollution conditions can have on the stations, we decided not to take into account the data of the second. The station used for reference, then, is found in an area of bare soil having a sufficient surface area for the Landsat 8-pixel (100 × 100 m) to allow for uniform and homogeneous results.
To ensure a greater number of environmental temperature comparison points, five certified temperature and humidity monitoring elements (Testo brand, model 184H) with calibration certificates were placed in a well-distributed manner throughout the city. In choosing the location of the measurement points, two criteria were established: (1) Aspire to an equitable distribution of probes in urban, suburban, green and rural areas. (2) Have consent from Granada´s City Council when the probes were situated on public urban furniture. Figure 1 and Table 2 indicate the positions selected for these five monitoring trials. A sixth measurement device was placed next to the meteorological station so as to check for possible differences or errors in readings that could invalidate the results obtained. Before commissioning, all the probes were set to measure the environmental conditions of a single point for two days, giving no differences in the measurements obtained.
The dataloggers were programmed so that they reliably registered the temperature and humidity values every 15 minutes during the entire study period. A total of 175,200 data for humidity and another 175,200 for temperature were thus obtained. They presented temperatures that oscillated between 253.15 K and 343.15 K, with an accuracy of ±0.5 K (273.15-343.15 K) and of ±0.8 K (253.15-273.15 K) with ±0.1 K. With regard to humidity, accuracy was ±2% to 298.15 K (20% to 80% DH) and ±3% to 298.15 K (20% and 80% DH). The datalogger confirmed calibration.

Methodology
The methodology applied in this research is depicted in Figure 2. After downloading the Landsat 8 images identified in Table 1, we processed and georeferenced them. Next, we derived the values of normalized vegetation index (NDVI), fractional vegetation cover (PV), top of atmosphere (TOA), at sensor temperature brightness (Tb) and soil emissivity required by the different algorithms (MW, CS, SW) to determine the LST. In order to validate the recovered LST values, we compared them with the ambient temperatures near the ground (1-2 m) obtained using temperature and humidity probes, and the correlation coefficients were calculated. Subsequently, and using the panel data technique, the relationships between the LST variable and the rest of the variables included in each algorithm were obtained in order to identify the variability of the data and their reliability. The use of this method in our research allowed us to reflect possible variations in the conditions of each variable contemplated in the final results pertaining to each algorithm, which makes it a powerful approach.

Determination of Soil Emissivity
LST recovery algorithms call for knowledge of the Earth's surface emissivity (LSE). This can vary significantly depending on the vegetation, humidity, roughness and other soil parameters [12,30]. The method usually used in the literature to determine the LST together with the LSE is through the NDVI THM calculation threshold [11,12,25,42,49,69]. Its use is considered suitable for the spectral ranges of 10 to 12 µm, among which lie the bands 10 and 11 of Landsat 8 [25]. This method was therefore selected for use in our study.

The Threshold of Calculation NDVI (NDVI THM )
The normalized vegetation index (NDVI) is calculated by means of the spectral bands of the near-IR spectrum (band 5) and the red band (band 4). This index gives an estimate of the presence of vegetation in a zone. The NDVI values range between −1 and 1, the former indicating clear, depopulated grounds, the latter, dense vegetation. The calculation formula is as follows: where 4 and 5 are the values of band 4 of bands 4 and 5, respectively. With the results obtained from Equation 1, the fractional vegetation cover (PV) can be calculated. This index establishes the proportion of a covered area for the vegetation or soil type [70,71]. Equation 2, for the calculation of the PV, is derived from the NDVI [11]: where NDVI is the index of vegetation normalized by means of Formula (1); and NDVI max and NDVI min are the maximum and minimal values of the NDVI interval. With the data obtained through Formulas (1) and (2), the emissivity of the soil can be calculated using Formula (3) for band 10 and Formula (4) for band 11 [11,15] (Table 3). Table 3. Estimation of ground emissivity normalized vegetation index (NDVI THM ).

Bands
Ground Emissivity NDVI 0.973 -0.0744 4 Here, εv and εs are emissivities of vegetation and ground, Ci is the coefficient of irregularity, and 4 is the value of band 4. The subindexes εv and εs indicated above correspond to Landsat 8 bands 10 and 11, respectively [11]. The Ci factor employed was Ci = 0 for flat surfaces [12].

Top of Atmosphere (TOA) and at Sensor Temperature Brightness (Tb)
The digital numbers (DNs) of the TIRS bands were converted into spectral radiance using the following equation [10,72]: where is the spectral radiance at the top of the atmosphere (TOA) expressed in W/(m²•sr•µm); is the specific multiplying factor of the band located in the metadata index of the Landsat 8 images. For the TIRS bands, this factor is 3.342 × 10 −4 W/(m²•sr• µm); is the digital variant value (DN) of the bands in the interval from 0 to 255, and is the additive factor of specific re-scaling of the TIRS bands, likewise included in the archives of meta-bits of information from the images.
From here on, the spectral radiance (Lλ) is converted and expressed as temperature brightness (Tb) in Kelvin. To this end, the conversion constants used are thermal K1 and K2 of the TIRS bands that appear in the file of metadata (band 10: K1 774.8553 and K2 = 1321.0789; Band 11: K1 = 480.8883 and K2 = 120.1442). The formula used is [10,72]:

Description of the Algorithms of Temperature Recovery
The first method used in the MW group to calculate the LST entails the validated algorithm [17]. The equation for its calculation is as follows: where LST is the land surface temperature, Tb is the temperature brightness of Landsat 8 TIRS 10, is the wavelength of the emitted radiation (10.8 µm for band 10) and ε is ground emissivity. In Equation 8, C2 = 1.4388 × 10 −2 m K, h is Planck's constant of 6.626 × 10 −34 Js, c is the speed of light 2.998 × 10 8 m/s, and s is the Boltzmann constant 1.38 × 10 −23 J/K [10,17].
The second method used involves the validated algorithm [19,28] calculated by: ( where a10 and b10 are the counters used to approximate the derivative of Planck´s radiance function for TIRS band 10 with respective values of -62.7182 and 0.4339, T10 is the Landsat 8 TIRS 10 temperature of brightness, Ꞇ10 is the atmospheric transmittance of Landsat 8 TIRS band 10 for the periods of summer and winter, ε10 is the ground emissivity for band 10, To is the air temperature close to the surface obtained with the datalogger, Ta is the effective atmospheric mean temperature in summer and winter, and is the total vapor content of the atmosphere [28,37].

Single Channel (SC) Algorithm
The method applied in the SC group to calculate the LST is the validated algorithm [16]. It is based on a polynomial approximation of the second order that depends on the content total of atmospheric vapor ( ) expressed in a matrix notation by means of Equation 19. LST is estimated as: where ε is the emissivity of the lot for band 10, LTOA 10 is the spectral radiance of band 10, bϒ is a constant value of 1324 for band 10, Tb10 is the temperature of brightness of Landsat 8 TIRS band 10, and is the atmospheric vapor.

Split Window (SW) Algorithms
This type of TIR band algorithm refers to two channels, 10 and 11, and is based on the attenuation of radiance generated by the atmospheric absorption, established in two wavelengths between 10 and 12 µm. The first SW algorithm is validated [16] through the expression: where Tb10 and Tb11 respectively represent the Landsat 8 TIRS temperature of brightness of bands 10 and 11, is the atmospheric vapor, ε and Δε are the values of the average emissivity and difference of emissivity between bands 10 and 11, respectively. The C coefficients used were: C0 = −0.268, C1 = 1.378, C2 = 0.183, C3 = 54.30, C4 = −2.238, C5 = −129.20 and C6 = 16.40 [16] according to the steam-driven status of water in the atmosphere marked by the probes.
The second SW algorithm used in this study, validated after the previous one, was based on [20] and is expressed as: where Tb10 and Tb11 respectively represent the Landsat 8 TIRS temperature of brightness of bands 10 and 11, ε and Δε are the values of the average emissivity and difference of emissivity between bands 10 and 11, respectively. The coefficients used were: b0 = −2.78009, b1 = 1.01408, b2 = 0.15833, b3 = −0.34991, b4 = 4.04487, b5 = 3.55414, b6 = −8.88394 and b7 = 0.9152, included in the status w of 0 and 2.5 gr cm −2 [20]. The last SW algorithm used to calculate the LST was validated with new coefficients [18]. The expression for the calculation of the LST is as follows: where Tb10 and Tb11 represent the temperature of brightness of Landsat 8 TIRS bands 10 and 11; ε10 and ε11 are the emissivity of the lot for bands 10 and 11. For the calculation of L10 and L11, the following expressions were used [15] for a temperature status between 293.15 and 323.15 K: = 0.4831 − 71.23.
Calculating Ꞇ10 and Ꞇ11 (atmospheric transmittance for bands 10 and 11) entailed the simulation of local weather conditions, mainly with steam-driven water content [18,19]. We adopted the North American standard of 1976 and the MODTRAN summer atmospheric profile of mid-latitude, with a steam-driven status of water in the atmosphere between 0 and 3 gr cm −2 [11]. Its expressions are: = −0.0164 − 0.04203 + 0.9715, where is the total vapor content of the atmosphere.

Identification of the Parameters for Each Algorithm
When calculating the LST, the algorithms (MW) must account for the parameters NDVI, VP, emissivity and temperature of brightness [17] and (SW) [16]. We added humidity to calculate the algorithms of the previous parameters (MW) [19,28], (SC) [16] and (SW) [18,20].
The reference serving to calculate air quality was the air pollution index (API), which is commonly used to measure air pollution [73][74][75][76]. Taking the form of a non-dimensional number in the range of 0 to 500, it provides a generalized degree of air pollution and the accompanying health hazards. Its statistical characteristics are described in Table 4.
The API is an integrated measure that reflects the levels of three fundamental air pollutants: SO2, NO2 and PM10. Each is calculated individually, and the overall API of the city is derived from the highest API of these three air pollutants. The interpolation to calculate the API of the city of Granada was performed using the following equation [74,76]: where API is the index for pollutant i (the most substantial of the three pollutants SO2, NO2, and PM10), Ci is the observed concentration of pollutant i, and Cu and Cl are the upper and lower limits of the range (shown in Table 5), within which the Ci, APIu and APIL are the upper and lower limits of the corresponding API interval (shown in Table 5). Finally, APILi is the value of the lower range when moving to a higher range [74,76].

Strategy of Analysis
Data analysis was carried out using specialized statistical software, STATA, version 15. The first stage of analysis involved simple regression between the LST values obtained by means of the validated algorithms and the LST obtained in situ. Secondly, intervening analysis of the data set was performed using the panel data technique to appraise the statistical relation between the parameters of each algorithm employed.
The panel data technique for intervening statistical analysis is widely used in the literature and includes multi-changing models of relation [77][78][79]. This method allows for a greater number of data points in the analysis, which increases the degree of freedom while reducing the inconvenience of collinearity between the variables [80]. By taking into account individual effects, the final function obtained for the set of individuals is totally different from the one that would have been obtained using other statistical techniques [81]. It is mainly indicated for dealing with time series involving multiple individuals and quantitative variables, where the possibility exists that they change the explanatory variables in the relations between individuals [82]. Hence, this system of analysis was very well suited to the experimental data of our study.
Introducing this method of statistical analysis in our model entailed two phases [77]. In the first place, by means of the Hausman proof [80,81], the effects of the analysis are determined to be either fixed or random. Then, the model is estimated based on the results obtained in the tests of Wooldridge and Wald [80,81]. There are three methods of calculation: method of ordinary squares (MOS), generalized least squares (GLS) and method of intragroup estimators (MIE) [81].
The first of the three, while widely used for years, does not enable the effects of every individual to be analyzed over the course of time, which can give rise to biased estimators.
The second one is considered to be a more efficient extension than the first; it is assumed that the individual effects are not reflected in the explanatory variables of the model. In this case, the individual effects contribute to the error term, the expression of calculation being: where represents the individual effects, is the error of the model, β is coefficient and X are explanatory variables, = individual and = time.
The third method cited above assumes that the individual effects are in line with the explanatory variables, so that the individual effect is separated after error, the expression of calculation being as follows: where, again, are the individual effects, is the error of the model, β is coefficient and X are explanatory variables, = individual and = time.

Results
Before comparing the LST obtained by means of the different algorithms, we proceeded to compare the temperature and humidity data obtained from the AEMET meteorological station and the attached probe. The mean differences detected for temperature and relative humidity over the entire period were, respectively, ±0.1 K and ±1%. The readings of the temperature and humidity probes used in the city can thus be considered adequate for later use as a validation method supporting the results of this research.

General Comparison of LST Algorithms
The statistics obtained using the algorithms to calculate the LST are reflected in Figure 3. Overall, the algorithms presented higher annual mean values than the ones obtained with dataloggers, except for those reported by authors Wang et al. (2015) [28], which were lower but very close values.  [18], which presented a greater variability of data among the top 25%.
As seen in Table 6, a larger difference in means was obtained by Wang et al. (2015) [28], with a value of -6.90 K. The smallest difference corresponded to the SWA study by Jiménez et al. (2014) [16] with 0.28 K. As for the groups of algorithms, MW presented a difference of −5.40 K, algorithm SC an average difference of 0.28 K, and SW a difference of 3.29 K. In terms of the standard deviation of the algorithms, LST algorithms based on MW and SC presented smaller standard deviation (10.61 K and 11.82 K, respectively) than the SW model (12.41 K). The method applied by authors Du et al. (2015) arrived at a much lower standard deviation using SW algorithms and a value similar to ours for the MW group.     Figure 5 presents the differences in mean temperatures between the LST calculated by means of the algorithms and the one obtained with dataloggers, distinguishing the winter period, with high humidity (October-April,) and the dry summer period (May-September). Figures 4 and 5 show that over the course of the winter, the LST calculated by means of algorithms was half the value obtained in situ by means of the datalogger. On the contrary, and during the summer period, the mean LST values based on the SC and MW algorithms were higher than those recorded in situ, although those of SW continued to be lower.

Distribution of LST in the City of Granada
As can be seen in Figure 6, the city of Granada has a longitudinal distribution, with more development in the east-west direction, where the city center and the neighborhoods of higher density are found. The highest temperatures are indicated in reddish and appear at the built-up areas of greater density (points 2 and 5 of Figure 6). They have a high percentage of waterproof grounds: asphalt, sidewalks, cobblestone pavement, and tiles. Contrariwise, the lower temperatures are indicated in blue, green and orange tones and are found in the suburban neighborhoods of average or low density (points 1 and 3) with more limited waterproof surfaces. Other light-shaded areas designate parks and "green areas" inside the city (point 4 in Figure 6) that act like cold poles or islands of urban cooling [82][83][84]. Table 7 offers the calculated statistical data of the LST for each algorithm and for each measurement point of the datalogger. The average LST of the urban areas was higher (299.73 K) than the LST of the suburban areas (299.10 K) and the LST of the large interior green area (298.39 K). The maximum thermal difference detected between the urban versus suburban area was 12.8 K, as compared to the difference between the urban area and the central green area, 13.15 K. Also detected was LST seasonal variability between urban areas and suburban ones. Accordingly, the average LST difference between the two zones in winter was 0.77 K, while this value reached 1.61 K in the summer. For the delimitation of urban, suburban and green areas, the different structures, coverage and densities that the city presents were taken into account. The steps carried out for the selection of zones were: (1) Compilation of the city´s metadata based on the real coverage through high resolution images from Google Street View and Sentinel 2 satellites; (2) Compilation of the metadata regarding building density in the city based on information from the National Institute of Statistics; (3) Selection and delimitation of zones.
According to the data obtained in points 1, 2 and 3, the zones were cataloged taking into account the greatest coincidence, not necessarily exact, of the sites with the zone classes. Table 7. Average values and temperature differences by points and algorithms.

Statistical Analysis of Data
The R 2 coefficients obtained for each algorithm are indicated in Figure 7. Overall, the algorithms used presented adequate R 2 coefficients of linear adjustment, with values over 0.95. This suggests good agreement between the analyzed values, and they are considered statistically significant (over 95%). The SW group presented a higher mean value with a very similar R 2 = 0.9764 when compared to the value of the SC algorithm, with R 2 = 0.9733. The MW group presented R 2 = 0.9684, the lowest value. Individually, the method used by Jiménez et al. (2014) [16] within group SW showed the highest value (R 2 = 0.9782), while the method applied by Qin et al. (2001) and  [28] gave the lowest value (R 2 = 0.9621).
The root mean square error (RMSE) and the mean bias error (MBE) values obtained are represented in Figure 8. The RMSE measures the error existing between two sets of data. The highest mean value corresponded to MW algorithms (RMSE = 2.20 K), an intermediate mean value to the SC (RMSE = 1.27 K), and the lowest mean value to SW algorithms (RMSE = 0.796 K). For the MBE the same classification was obtained, the highest mean value corresponding to MW algorithms (MBE = −0.077 K), an intermediate mean value to the SC (MBE = −0.034 K), and the lowest mean value to SW algorithms (MBE = −0.031 K). The group of algorithms with a lower R 2 value therefore showed a higher RMSE and MBE, and vice versa. To carry out statistical analysis using the panel data method, it was necessary to determine if calculation should be done by means of fixed or random effects. Hausman's proof [79,80] [16].
To develop the panel data, the generalized least squares method was applied through Equation 34; the results are shown in Table 8. For each algorithm of calculation, positive and negative correlations were found between the independent variables and the dependent variable. Independent variables did not correlate with the dependent variable in the same way under any method of calculating LST: there were positive and negative series of alternating marks and spaces. Although the independent variable temperature of brightness (Tb) correlated positively in all the methods, the rest presented alternations of positive correlations or negative ones, depending on the algorithm used. This variable correlated positively under all the methods and, moreover, at a level of significance over 95%, since in all cases p < 0.05. Table 8. Results of the statistical analysis using the panel data.

Group
Author  [25] gave the largest number of negative correlations of the independent variables with the dependent one, a total of 4: humidity, NDVI, VP and constant value. Nevertheless, it is important to underline that they did not all do so in the same way. Humidity and the constant value presented a greater correlation, with p = 0.000 and statistically significant to a level over 99%. The independent variable NDVI also correlated negatively and was considered statistically significant, but to a level of significance over 95%, as 0.05 < p < 0.01. In turn, the variable VP correlated negatively, but it was not statistically significant (p = 0.832).
The algorithm used by Mao et al. (2005) [18] was the one that presented the fewest negative correlations, just one: emissivity. This variable was negatively correlated with the LST and was found to be statistically significant to a level over 95%, since p < 0.05.
The values of R 2 obtained by means of the panel data are given in Table 9. They were adequate, all reaching over 0.97. This finding suggests good concordance between the dependent variable and the independent ones of the methods analyzed, given the level of statistical significance over 97% in every case. Table 9. Linear adjustment coefficients R 2 using the panel data.

Group
Author R 2 Average R 2 by groups The highest R 2 value was obtained with SC, using the method of Jiménez et al. (2014) [16]. In this case, 99.97% of the data variability was explained by the algorithm, the reliability of all variables being over 99%. In turn, the lowest R 2 value corresponded to the method of Weng et al. (2004) [17], where the algorithm accounted for 97.07% of the variability of the data, and the reliability of all variables was over 99%. That is, the independent variables correlated better with the dependent variable in the first algorithm than in the second one.
Therefore, the greater the number of negative correlations between the independent variables and the dependent one of an algorithm, the lesser the variability of the data, hence a smaller value for R 2 . Among the groups of algorithms, SC presented the highest value (R 2 = 0.9997), MW gave an intermediate value (R 2 = 0.9843), and SW presented the lowest value (R 2 = 0.9838), though very close to MW.
The results obtained with respect to the contamination for each LST recovered are given in Table 10. In view of the statistical analysis, the pollution factor presented a significant relationship at 99% (p < 0.01) with all the LST recovery algorithms used. The values obtained for R 2 and the F statistic indicated good agreement between the LST algorithms and pollution for each method used, with an adjustment level up to 99% of significance, since Prob > chi 2 < 0.01.
There are numerous investigations [12,79,85,86] that have related an increase in LST and environmental temperature with higher concentrations of pollution, especially in urban areas. This is due to the fact that short-wave radiation from the sun passes through the atmosphere and heats the earth's surfaces. The waves give off infrared radiation that escapes into the atmosphere and heats it up. High concentrations of pollution prevent the escape of radiation, conserving heat that warms the lower layers of the atmosphere and the land surface [85,86]. Table 10. Results of the statistical analysis between LST and contamination by algorithm.

Group
Author

Discussion
In this research, mean LST differences obtained by algorithms were lower than the ambient temperatures measured near the ground. However, the LST differences derived using algorithms were uneven, SW methods giving values above the mean, while MW values were well below the mean. The SC method gave a value around the mean. Therefore, the values obtained from the SW and SC algorithms largely agreed with data reported previously [40,41,[43][44][45] in studies of other cities or territories using the comparison of temperatures near the ground as a validation method. These studies gave LST results superior to those obtained by environmental temperatures ranging between 0.7 and 1.0 K in the city of Hong Kong [40], 3 and 7 K in studies of cities in Senegal [45], 0.2 and 7.8 K in studies of various cities in Canada [41], the United Kingdom and the USA, where the authors [43] reported differences between 2.58 and 7.62 K. On the contrary, the MW algorithms presented differences that were below ambient temperature. This circumstance was already reported by the authors of [42] in their LST studies using MW algorithms in India, where they obtained differences between 3 and 4 K lower than ambient temperatures.
Evidence of the relationship between pollution in the city of Granada and the LST recovered by all algorithms was found. This relationship has been reported by numerous investigations [12,79,85,86], validating the results presented here. In light of the LST results cited above, and the relationship between LST and contamination, it cannot be said that contamination affects the LST recording using Landsat 8 algorithms.
We found that the temperature curves derived from algorithms MW and SC were similar in shape and intensity. Contrariwise, the curves resulting from SW algorithms differed, above all in the summer season when the air humidity (w) diminishes. A similar finding was reported by other authors, with the same distribution [12][13][14]28,50,72]. It can be attributed to the number of TIRS bands that each algorithm uses-algorithms MW and SC use the data of a single TIRS channel (band 10), while SW algorithms rely on two TIRS (bands 10 and 11) with wavelength intervals that oscillate between 10.6 and 11.19 µm, and 11.50 and 12.51 µm, respectively. The uncertainty that could be generated in band 11, implying possible error in the algorithm of correction, contributes to this situation (SLCA). Although band 11 was in working order in February 2017, numerous authors insisted on the need to conduct further studies on validation [12,13,25,87] to verify its correct functioning.
The panel data for the study of relations between the dependent variable and the independent variables of every algorithm indicated they all yielded values with a level of significance over 95%. This means that the variability of the data accounts for an important percentage of the algorithm. Still, algorithms MW and SC presented R 2 values superior to those of the SW algorithms. The variability in the data of the MW and SC is therefore better explained by their algorithms. Because the data of a single TIRS channel were used for the algorithms MW and SC, the possible uncertainty of band 11 used for SW appears more likely.
In turn, the more marked differences seen in the LST curves of the SW algorithms with respect to the others over the whole period of study can be interpreted as an overestimation of LST with low values of relative humidity (w), and an underestimation of LST with high humidity. Such a situation was described by García et al. (2018) [13], who justified the bias of LTS for bands 10 and 11. This could also explain the standard deviations of the LST averages from season to season: larger values correspond to periods of low humidity (summer), and lower values to periods of high humidity (winter) [12].
The data obtained in the city of Granada reveal an important distinction in LST between industrial and urban zones and suburban and green zones, along with the influence of the seasonal period. The maximum value detected between the two zones was 12.8 K, very similar to figures obtained by other authors [7,12,35,50,88]. Lower values have been found in conjunction with zones having abundant green groundcover [17]. The heterogeneous system of the city studied here allows for the distinction of industrial and urban areas with high LST and high density, plus vast surfaces of waterproof ground, as opposed to suburban areas with lower LST, low density, and limited waterproof surfaces. The lowest LSTs corresponded to surfaces with vegetation. The average LST differences between the urban and suburban areas of Granada increased by 50% in the summer period (as opposed to winter). Such findings are in line with the findings of other authors [10,15,49]. The direct relationship between ground permeability and NDVI implies lower LST values in areas with high NDVI, and the other way around. Numerous studies corroborate this and validate the data obtained [5,7,10,12,[88][89][90]. The results obtained in the determination of the LST of the Spanish cities of Barcelona [49,50] and Madrid [51] showed values that are very similar to those reported here. Thus, the hot zones identified in Barcelona and Madrid were the industrial ones and the urban fabric, in contrast to the colder zones located in the suburban zones and green zones [49][50][51]. The mean differences detected between the two areas in both cities were between 1 and 2 K [49][50][51], very similar to the values obtained in our investigation (1.34 K). The LST standard deviation values obtained in our research (winter = 2.39 K and summer = 2.72 K) are also very similar to those obtained by the authors [49] for the city of Barcelona (winter = 2.3 K and summer = 3.1 K). The LST studies carried out on the city of Bangkok [7] established that the highest LST was obtained in urban areas and not in industrial areas. This circumstance is motivated by the fact that Bangkok has an urban area of approximately 90%, reducing the industrial area to only 10%.
The R 2 determination coefficients obtained by means of the panel data technique for the compared algorithms were above 0.97, which denotes a good agreement of the values.  [11,14,15,17,37]. One study [13] obtained values higher than ours, but this could be because the area under study was completely vegetal. With respect to the MBE data These values denote a good agreement of data for the SC and SW algorithms, and low agreement in the data of the MW algorithms. Additional investigations [11,14,15,17,37] have presented similar values, thus supporting the results obtained.

Conclusions
This is one of the first multi-temporal experiments aimed to determine LST with the Landsat 8 images of a Spanish city. Although such trials are common in large cities, in recent years they are becoming somewhat more frequent in mid-size cities. This circumstance, together with the morphologic and geographic characteristics of Granada, its levels of pollution and high thermal contrast, heighten interest in the results of our investigation.
According to our results, the LST of Granada can be determined using six algorithms under the categories MW, SC and SW, validated by intervening in situ, at five points, with datalogger measurements of temperature and humidity. Quantitatively close distributions and similar behaviors were observed for the city overall. However, important distinctions were seen between urban and suburban zones, and seasonal periods generated heterogeneous distributions. Even so, the results are largely in consonance with figures obtained at other geographic points.
In our study, SW algorithms gave higher LST values than MW and SC. We believe that these results can be applicable to other cities; indeed, they are similar to those obtained in investigations carried out elsewhere. The former uses the data from two thermal bands, the latter two using the data of a single band. The possibility of uncertainty in band 11 data can be corrected by the SLCA algorithm. The panel data method served to confirm that the variability in relations between the dependent variable and the independent variables of each algorithm is greater when a single channel (rather than two bands) is used.
The difference in means between the calculated LST and the LST in situ found here was lower for SC models (0.28 K), intermediate for SW (3.29 K) and highest for MW (−5.90 K). Coefficients of determination R 2 of the LST obtained from algorithms and the LST in situ demonstrated high correlation, over 95%. Nevertheless, the SW algorithms presented larger coefficients than MW. Finally, RMSE and MBE values were inverse to the R 2 , being lower for SW algorithms and higher for MW. No evidence was found regarding the effect of high contamination on the process of obtaining the LST by means of recovery algorithms.
In our research, the SC algorithm gave the values closest to the ambient temperature obtained using dataloggers, and the values of RMSE, MBE and R 2 were considered adequate. Therefore, it is deemed to be most suitable for evaluating the LST of the city of Granada, a medium-sized city with high levels of pollution and notable thermal differences, using Landsat 8 images. Even so, we consider it necessary to carry out more research to validate the precision of these algorithms and determine their possible impact on the contamination phenomenon.
The present contribution is significant in the context of LST determinations in midsized cities by means of Landsat 8 images. In view of the advantages and disadvantages encountered in this research regarding the selected algorithms, the inclusion of six algorithms could prove most useful for future LST studies of other urban areas. This novel, low-cost, fast and accurate approach allows one to monitor variations in the LST of cities in a precise manner, taking into account physical, environmental and climatic variables. Ongoing efforts are aimed to complement this study. Future studies involving the LST of other cities worldwide will underline the importance of validating the precision and maturity of algorithms used in multi-temporal assessments. Ultimately, this should contribute to well-founded decision-making by urban planners and public administrations in view of the distribution of the LST in cities. Data Availability Statement: Data from this research will be available upon request to the authors.