Retrieving Land Surface Temperature from Satellite Imagery with a Novel Combined Strategy

Land surface temperature (LST) is a key parameter for land cover analysis and for many fields of study, for example, in agriculture, due to its relationship with the state of the crop in the evaluation of natural phenomena such as volcanic eruptions and geothermal areas, in desertification studies, or in the estimation of several variables of environmental interest such as evapotranspiration. The computation of LST from satellite imagery is possible due to the advances in thermal infrared technology and its implementation in artificial satellites. For example, Landsat 8 incorporates Operational Land Imager(OLI) and Thermal InfraRed Sensor(TIRS)sensors the images from which, in combination with data from other satellite platforms (such as Terra and Aqua) provide all the information needed for the computation of LST. Different methodologies have been developed for the computation of LST from satellite images, such as single-channel and split-window methodologies. In this paper, two existing single-channel methodologies are evaluated through their application to images from Landsat 8, with the aim at determining the optimal atmospheric conditions for their application, instead of searching for the best methodology for all cases. This evaluation results in the development of a new adaptive strategy for the computation of LST consisting of a conditional process that uses the environmental conditions to determine the most suitable computation method.

assess the potential of geothermal areas [8,9], or estimate the evapotranspiration of vegetated areas for different environmental studies.
The increase in earth surface temperature is a global dilemma as indicated by one of the sustainable development goals (SDGs) promoted by the United Nations [10]. Thus, monitoring this phenomenon in order to assess changes through time is key and should be performed in any part of the earth for large geographical areas. This is possible due to remote sensing satellite imagery. In fact, remote sensing is positioned as the main alternative for quantification and monitoring of the LST [11].
There are different approaches developed to calculate LST from satellite imagery using the thermal infrared (IR) spectrum. The difference between them lays on the amount and type of data used by their algorithms, especially on whether the algorithm includes the computation of the surface emissivity parameter or if it is introduced as a measured parameter [12]. In order to calculate LST, surface emissivity (also known as land surface emissivity (LSE)) can be known or unknown. In the latter, different methods are used as follows: (1) stepwise retrieval methods, such as the classification-based emissivity method (CBEM) [13]; (2) methods for simultaneous retrieval of LST and LSE with known atmospheric parameters, such as grey-body emissivity (GBE) [14] and the temperature emissivity separation (TES) algorithms [15]; and (3) methods for simultaneous retrieval of LST, LSE, and atmospheric profile, as done by artificial neural networks (ANNs) [16]. If the LSE is known, the existing methods differ on the number of IR thermal bands used. Thus, for the case of using a single IR thermal band, single-channel methods based on the resolution of the radiative transfer equation and the Planck equation are applied [17]. These methods attenuate the scarce terrain information by introducing several atmospheric parameters, such as indirect measurements of air temperature and humidity [18]. On the other hand, for the case of using more than one thermal band, split window [17] or multi-angle [19] algorithms are applied. These methods differ in the way input data is required. Split window requires differential absorption data while multi-angle requires imagery acquired at different viewing angles. All these methods require additional atmospheric corrections. These corrections can be based on images, as in the case of DOS correction [17], or on the modelling of the atmosphere, as in the codes of radiative transfer 6S [20] or MODTRAN [21].
Since 2018, the United States Geological Survey has proposed a systematic treatment for all American states named Landsat Analysis Ready Data (ARD), a 30-m spatial resolution LST product from Landsat 8 and the United States area [22]. In this regard, there are no direct LST products available anywhere else in the world with the same spatio-temporal resolutions. Regarding the most similar options, the European Space Agency has launched daily LST products as Level 2 data from Sentinel 3 satellites [23], but their spatial resolution is 1000 m, 33 times lower than data from Landsat 8. In the Sentinel 3 case, the methodology used is split-window based, thus being applicable only to sensors acquiring TIR (Thermal InfraRed) data in several bands, while ARD LST data is calculated with single-channel methodology. ASTER sensors also acquire TIR data and allow for the generation of an LST product with a resolution of 90 m. On the other hand, ASTER sensors from NASA and the Japan Ministry of Economy, Trade, and Industry offer this product for specific data collections [24] without a real global coverage.
Thus, this paper aims to find the optimal atmospheric conditions for the application of two existing single-channel methods for estimating LST from Landsat 8 data, but applicable to other satellite platforms, to reproduce this result anywhere in the world. For this reason, the Landsat ARD LST product is considered as a reference to evaluate the results obtained with the strategy proposed. In addition, the fact that the ARD LST products have been widely evaluated and validated reinforces their use as reference products [25,26].
Based on this, the present work seeks to make an in-depth study of the existing single-channel methodologies for computing LST based on satellite thermal data. It differs from similar works, such as [27], in the fact that the focus is set on the development of an adaptive strategy to apply the best methodology based on the atmospheric conditions of each case study, instead of searching for the best method in a general way. In this case, Landsat 8 and MODIS sensors onboard Terra were chosen, but the process can be extrapolated to any other satellite sensor that offers thermal, visible, and near-infrared data, such as ASTER. The results support the possibility of applying a combined strategy, which selects automatically the appropriate method based on the environmental conditions of each pixel.

LST Retrieval Methods
There are two different types of methodologies to calculate LST from satellite images, as shown in Figure 1; on the one hand, there are methods that use previously "known" or rather approximate emissivity, and on the other hand, there are methods that consider emissivity as unknown [12]. This paper focuses only on those methods in which emissivity was previously known, and more specifically on single-channel methods. The reason for working with these methods was the fact that they can be applied to many different satellites, from those that acquire in only one range of the thermal IR band, to those that acquire images in two or more bands in the thermal IR range. When working with Landsat 8 in this article, the use of thermal band 10 was recommended, since band 11 presented problems in its calibration as indicated by different investigations [18].

Single-Channel Methods
Single-channel methods include those that use a single thermal band. In our case studies in which the Landsat 8 satellite was used, band 10 was chosen.
Within the single-channel methods there were two different methodologies considered in this paper: Direct inversion of the radiative transfer equation: the radiative transfer equation and the Planck equation are known, such that being: L sensori : the spectral radiance obtained by the sensor, measured in W·m 2 ·sr -1 ·µm -1 . B i (T s ): the measured radiance if the surface were a black body at the surface temperature T s . τ i : the total transmissivity of the atmosphere for the channel I of the sensor. ε i : the surface emissivity. L ↓ atm i : the ascending atmospheric radiance. L ↑ atm i : the descending atmospheric radiance. λ: the wavelength of the band i. k: the Boltzmann constant, k = 1.3806 × 10 -23 m 2· kg/s 2 ·K. h: the Planck constant, h = 6.6261 × 10 -34 m 2 ·kg/s. c: the speed of light in a vacuum, c = 2.9979·108 m/s. C 1 : the first radiation constant equivalent to 1.19104·10 8 W·µm 4 /m 2 ·sr. It is calculated from C 1 = 2π·h·c 2 . C 2 : the second radiation constant calculated as C 2 = h·c/k with a value of 1.4388·10 4 K·m.
In this case it would be necessary to know the emissivity of the earth's surface and the different atmospheric parameters: transmissivity and atmospheric radiance both ascending and descending. These parameters can be obtained from atmospheric profiles and radiative transfer codes. As an example of the first, there are radiosoundings, reanalysis, and satellite sounders [28]. In the case of radiative transfer models, some examples are MODTRAN, RTTOV, and ATCOR [21].
Generalized Single-Channel Algorithm: These algorithms were developed to simplify the radiative transfer equation. An attempt was made to avoid dependence on atmospheric profiles due to the fact that they are not always available or are too complex to acquire, and instead algorithms were developed that use other atmospheric parameters such as water vapor or air temperature. The original algorithm that all single-channel methods have in common was developed in [29], and its formula is as follows: C2: the second radiation constant calculated as C2 = h·c/k with a value of 1.4388·10 4 K·m.
In this case it would be necessary to know the emissivity of the earth's surface and the different atmospheric parameters: transmissivity and atmospheric radiance both ascending and descending. These parameters can be obtained from atmospheric profiles and radiative transfer codes. As an example of the first, there are radiosoundings, reanalysis, and satellite sounders [28]. In the case of radiative transfer models, some examples are MODTRAN, RTTOV, and ATCOR [21].
Generalized Single-Channel Algorithm: These algorithms were developed to simplify the radiative transfer equation. An attempt was made to avoid dependence on atmospheric profiles due to the fact that they are not always available or are too complex to acquire, and instead algorithms were developed that use other atmospheric parameters such as water vapor or air temperature. The original algorithm that all single-channel methods have in common was developed in [29], and its formula is as follows: where ϒ and δ, the parameters obtained from the Planck function, are defined as follows: Tb: Brightness temperature in the sensor in Kelvin (K). ψ1, ψ2, and ψ3: the atmospheric functions. The values of these last functions are the difference from one method to the other. For this paper the generalized single-channel algorithm was used. There are many single-channel algorithms, but two of them are exposed in detail, for which they are called single-channel 1 and single-channel 2, the main difference of which lays in the atmospheric functions used.

Single-Channel 1
This method is based on [17], with the pipeline shown in Figure 2. and δ, the parameters obtained from the Planck function, are defined as follows: Tb: Brightness temperature in the sensor in Kelvin (K). ψ 1 , ψ 2 , and ψ 3 : the atmospheric functions.
The values of these last functions are the difference from one method to the other. For this paper the generalized single-channel algorithm was used. There are many single-channel algorithms, but two of them are exposed in detail, for which they are called single-channel 1 and single-channel 2, the main difference of which lays in the atmospheric functions used.

Single-Channel 1
This method is based on [17], with the pipeline shown in Figure 2. To execute this method, four different images are needed: • Image of the red band (0.62-0.70 µm). Band 4 for Landsat 8.

•
Preprocessed products of column water vapor (W) given in g/cm 2 . For this article, products of the MODIS sensor were used, specifically of the Terra satellite, since this follows a similar path to that of Landsat 8, with an approximate lag of 30 minutes. Specifically, the MOD05_L2 product was used [30].
For the single-channel 1 method, the values of the atmospheric functions were obtained with a second degree polynomial in which the water vapor content is related, since, as indicated previously, this is one of the parameters with higher contribution to the atmospheric effect in the region of the thermal infrared.
According to [17], LST errors for this method become unacceptable for high water vapor contents (for contents greater than 2.5 g/cm 2 ). It also indicates that this problem can be solved by including air temperature as an input, as was done in more recent works [29].

Single-Channel 2
The working scheme is the same as for the single-channel 1 method, so the same inputs were needed. The difference is that this method was based on [21] and works with a third-degree polynomial to obtain the atmospheric functions. The original method published in [31] but corrected for Landsat 8 was analyzed in this paper. Therefore, the atmospheric functions in this case are The atmospheric functions incorporate spectral functions, which are defined for each wavelength as shown in Table 1.

AF Spectral Functions
Taking into account that this work used band 10 of Landsat 8 and that its average wavelength was 10.8 µm, the atmospheric functions are defined as According to [17], the most important source of error in single-channel algorithms (in general, as regards the single-channel 1 and single-channel 2) is the introduction of the atmospheric effects, which leads to an error in the LST between 0.2 K and 0.7 K, and the uncertainty of the emissivity of the earth's surface, which produces an additional error in the surface temperature between 0.2 K and 0.4 K.

Methods
Once the most recommended methods for calculating LST by means of satellite images from only one thermal band were determined from a theoretical point of view, their evaluation was carried out in order to determine the most appropriate conditions for the application of each of them. To do this, throughout this section we proceed to the automation of the processes and the visualization of the results making combined use of the following two tools with great potential: Matlab ® as a programming language for process automation, and QGIS ® as a Geographic Information System for geospatial data processing. As a result of this combination, it was possible to perform an evaluation of the different methods mentioned, both from a numerical and a visual point of view. The products offered by NASA called analysis ready data (ARD) were considered as reference for the evaluation of each method under study in order to check their viability regarding the most official temperature product.

Landsat Analysis Ready Data (ARD)
ARD data was published in 2018 and created by NASA. Their bases were Landsat Level-1 products, including TOA (Top-Of-Atmosphere) reflectance, satellite brightness temperature, surface reflectance, and provisional terrestrial surface temperature.
Landsat ARD LST products are currently only available for the US. Images from the continental United States and Hawaii are generated in the Albers Equal Area (AEA) conical cartographic projection, using the WGS84 (World Geodetic System 1984) as datum.
Each ARD LST product has a spatial resolution of thirty meters, just like the original Landsat 8 products, but the size is smaller. Landsat ARD LST has 5000 × 5000 pixels versus approximately 6300 × 6000 pixels in the other products from Landsat 8.
The ARD LST products are different from those calculated in this article, mainly due to the processing procedure. Regarding the method for atmospheric correction of the optical spectrum, Landsat ARD generates surface reflectance using the Landsat Surface Reflectance Code, which uses the coastal spray band, MODIS auxiliary weather data, and a unique model of radiative transfer [32], while this paper implemented DOS 1 atmospheric correction [17], based on the MODIS atmospheric water content product. In addition, ARD LST products are generated using ASTER GED (Global-Emissivity Dataset) for surface emissivity and MODTRAN simulations powered by the atmosphere proposed by the NARR reanalysis [21], while the single-channel approaches used in the combined strategy retrieve surface emissivity by the NDVI threshold method to calculate surface emissivity.

LST Retrieval Methods: Evaluation and Comparison
The objective of the evaluation was to determine if the final result of each method, that is, the LST value obtained, approximated the value considered as reference. As mentioned before, the provisional surface temperature image of ARD was considered as reference.
For this comparison, sample images were randomly selected from locations distributed as homogeneously as possible in the USA and on different types of soil, including cases studies of desert, urban, crop, and forestry areas (Figure 3), taking into account case studies both in warm and cold seasons ( Table 2). The images were cut in order to avoid the presence of urban areas, the particularities of which regarding temperature distribution and emission of thermal infrared radiation made necessary the adaptation of the existing LST retrieval methods [33,34].  This means that three different locations and two different dates are studied for each area. For this, 72 Landsat 8 images (24 sets of bands 4, 5 and 10) and their metadata will be downloaded. In addition, 24 MODIS water vapor products and the corresponding 24 ARD LST products are required for the study. Dates selected where those with the lowest percentage of cloudiness possible (below 2%). Tables 3 and 4 show the numerical results obtained from the comparison between ARD LST products and the LST images resulting from the application of single-channel 1 and single-channel 2 methods. RMSE was calculated per pixel for each image resulting from both methods. Then, the mean RMSE was calculated per image and for the pixels of all images from the same method. Table 3. Results of the mean square error between the analysis ready data (ARD) LST products and the different methods of LST retrieval in warm season.

Crop Area
Forest Area For each sample area, LST was calculated using the different methods considered for this work (single-channel 1 and 2), calculating the difference with the ARD LST product, and obtaining its root mean square error. The surface temperatures and average water vapor content of each area were also indicated.

Determination of Conditions for the Optimal Application of the Methods
The results obtained in Tables 3 and 4 were sorted according to the average water vapor content and the methods with better and worse results, as seen in Table 5. The results could also be ordered taking into account their average terrestrial surface temperature, as shown in Table 6 and Figure 4. Table 6. LST retrieval methods in the thermal spectrum that are closest to the results of ARD LST products for the case studies. Images are sorted by the average surface temperature ARD of each case study.  Analyzing Tables 5 and 6 and Figure 4, the following conclusions could be reached. In the case of water vapor content, the results are as follows:

Mean LST Range (K) Number of Images Best Methodology
For water vapor content below 1.2 g/cm 2 , the most approximate method compared to the results of NASA is the single-channel 2.
For water vapor content greater than 1.8 g/cm 2 , the most approximate method is the single-channel 1.
For intermediate water vapor content values, the results vary between both methods. The same happens if the analysis focuses on the surface temperature, as follows: • For surface temperatures less than 295 K the most approximate method is single-channel 2.

•
For temperatures greater than 302 K the closest method is single-channel 1.

•
For intermediate values there is also variation.
The analysis of Tables 3 and 4 verified that, as a general rule, the higher the temperature, the higher the humidity, so that the characteristics necessary for the choice of the method are shared, the results of which are the most similar to those of the ARD LST products. With this information the combined strategy can be developed. Tables 5 and 6 show a trend in the results; for images acquired with high levels of water vapor content (1.8 g/cm 2 ), the most accurate method was single-channel 1 (assuming as more accurate the results that were closer to the values of the ARD LST products). However, when water vapor content was low (1.2 g/cm 2 ), the best method for LST computation was single-channel 2. It should be noted that the results presented high RMSE in cases where water vapor content was very high. This fact goes in coherence with studies such as [13], where results of LST when water vapor content was over 2.5 g/cm 2 were considered as not reliable. For intermediate values of water vapor content, no clear trend was found, with the best results oscillating between single-channel 1 and single-channel 2 methods.

Results: Combined Strategy
Taking into account the results of the analysis, a combined strategy was developed that analyzed the image and computed LST pixel per pixel, according to the rule shown in Figure 5. Thus, the first parameter under study was water vapor content (W), with the following decision-rule: W > 1.8 g/cm 2 : single-channel 1 method, W < 1.2 g/cm 2 : single-channel 2 method.
For intermediate values of W (between 1.2 and 1.8 (1.8 g/cm 2 ), the second parameter under study was satellite brightness temperature (Tb), which was closely related to LST. The logic applied was the following: Tb > 295 K: single-channel 1 method, Tb ≤ 295 K: single-channel 2 method.
In this way, the combined strategy consisted of an adaptive pixel-based LST retrieval methodology. Thus, the strategy analyzed the TIR image from which LST values needed to be retrieved, pixel by pixel. For each pixel, the water vapor content was analyzed, and the single-channel 1 or single-channel 2 method was used for the computation of LST. If the water vapor content was not decisive (that is, the value was between 1.2 and 1.8 g/cm 2 ), then the surface temperature was studied. Since the study was made pixel-by-pixel, both single-channel 1 and single-channel 2 could be applied in the same image, in different pixels each.
The combined strategy was subjected to both quantitative and visual evaluations, with the aim of determining its performance.

Quantitative Evaluation
The results of LST obtained with the combined strategy were compared to ARD LST products considered as reference product.
This comparison was performed through the automation of the image processing and the application of the combined strategy in Matlab ® . The images used for the validation of the combined strategy were the sample of 24 images used in the previous sections of the paper. Results are shown in Tables 7-10.   Table 9. Results of root mean square error (RMSE) of the combined strategy with respect to ARD LST products in desert areas.  Results showed that the mean RMSE for all areas of the combined strategy was 1.20 K, while the mean RMSE of the optimal single-channel method for each image were 1.51 K and 2.04 K, for single-channel 1 and 2, respectively. Thus, the combined strategy implied a mean improvement of 20% and 41%, respectively.

Visual Evaluation
Provided that both the inputs and the outputs of the computation of LST from satellite imagery present image-nature, visual analysis of the results is required. This way, for example, the coincidence of location of the coldest/hottest areas can be detected, as well as the identification of the magnitude of the error (the same RMSE can be obtained for a small difference for all pixels in the image, and for a large difference for a small number of pixels, and this can be identified visually). Tables 7-10 show results of improvement, worsening, or selection of the wrong method. These are all evaluated visually in this section. Figure 6 shows an image where the combined strategy presents the same RMSE as the optimal single-channel method, which corresponds to image N from a desert area.  Figure 6 shows that images resulting from the single-channel 2 method (c) and the combined strategy (d) are equal. In addition, the images resulting from these methods are more visually similar to the ARD LST product (a) than the image resulting from the single-channel 1 method (b). This similarity can also be evaluated through the generation of images of difference of temperatures, as shown in Figure 7. Next, a case study where the combined strategy does not improve the results of the best single-channel method is shown. This case corresponds to image C, corresponding to a crop area. According to Table 3, the best method is single-channel 1, with an RMSE of 1.07 K. Figures 8 and 9 show that the RMSE of the combined strategy is higher than that of single-channel 1, but results are better than those from single-channel 2 (2.06 K).  Figure 8 shows no radical difference between temperatures obtained with the single-channel 1 method (b) and the combined strategy (d). It is also not clear that the single-channel 2 method is the worst (c). Figure 9 shows this comparison more clearly. Figure 9 shows that the RMSE of the combined strategy is a combination of the RMSE of single-channel 1 and 2 methods. If the image of water vapor content is studied (Figure 9d), the coincidence of the pixels with lower vapor content with the RMSE of the single-channel 2 method is unequivocally identified. Thus, this image is a clear example of the pixel-by-pixel analysis basis of the combined strategy.  Figure 10 shows a case study where the combined strategy presents poor results, with mean RMSE over 2 K. The reason for this is that the single-channel method did not present accurate results (Figure 10b,c). This is due to the high water vapor content in the area of over 2.5 g/cm 2 . This case study serves as a confirmation that values of water vapor content over 2.5 g/cm 2 are the limit for the application of the combined strategy proposed in this paper.

Application of Combined Strategy to Other Sensors
The combined strategy was developed and tested with images from the Landsat 8 satellite. However, its application is not limited to this satellite, being usable by other satellites that incorporate the visible, near infrared, and thermal infrared bands needed. In fact, original single-channel methods were implemented for satellites earlier than Landsat 8 [12,18,21].
Some studies focus on the computation of LST from other satellites such as ASTER [35,36]. These use the same single-channel methods analyzed in this paper, adapted to the wavelength of the bands in ASTER. This implies that the combined strategy can also be applied to ASTER images, with the corresponding adaptation of the algorithms.
In order to evaluate the applicability of the combined strategy to the retrieval of LST from ASTER images, the resulting LST image was compared to AST 08V003 product [37], which is the ARD LST generated by ASTER. The results obtained are shown in Table 11. AST 08V003 was used as reference in this case in order to make a proper comparison of the combined strategy proposed with the official product generated from the same raw data, avoiding the effect of the different day and hour of acquisition that would have been present if the reference had been the ARD LST product from Landsat 8. In this case, the validity of the combined strategy proposed implies a simplification of the methodology for LST retrieval, since the official methodology used is a multi-channel temperature-emissivity separation (TES) algorithm that requires a minimum of three bands in the TIR spectrum [38], and the strategy proposed only needs one TIR band. In this case, the combined strategy presents the same RMSE as the single-channel 1 method, consequently not worsening the results of none of the existing methods, and improving the results by 30% with regard to the worst method ( Figure 11 shows the visual results of the application of each method). The results of the analysis of the application of the combined strategy to ASTER images show that the method presents better results than the methodologies officially used, also for satellites other than Landsat.

Conclusions
The need to perform global measurements of LST was the main reason for the evaluation of the viability of satellite imagery for the task. However, the number of methodologies available for the task is constantly increasing, the single-channel methods being the ones most widely used. The reason for the latter is that single-channel methods allow the computation of LST from one image in the thermal infrared band only, provided that a set of additional parameters (atmospheric conditions, surface emissivity) is previously known or calculated, in such a way that the method can be applied to any satellite with at least one thermal infrared band and is not limited to those sensors including several acquisitions within the band. However, existing single-channel methods differ from each other by the procedures applied to obtain the required additional parameters, which can also be different for each method, and there is no clear determination of the most precise and accurate method for each situation.
Thus, this paper performed an analysis of the most widely used single-channel methods, resulting in the development of a "combined strategy". The combined strategy uses both single-channel 1 (second grade polynomial including water vapor content (W) as parameter) and single-channel 2 (third grade polynomial with water vapor content as parameter) methods. The combined strategy was developed by comparing the LST images resulting from each method (single-channel 1 and 2) with the official ARD LST images for each sensor (Landsat and ASTER), which were considered as reference products. Results showed that when W was over 1.8 g/cm 2 the best method was single-channel 1, while when W was below 1.2 g/cm 2 single-channel 2 generated better results. For values of W between 1.2 and 1.8 g/cm 2 , satellite brightness temperature (Tb) was the parameter that influenced the performance of the LST computation methods. Thus, the analysis performed showed that Tb over 295 K was better dealt with using the single-channel 1 method, while Tb below 295 K had better results with the single-channel 2 method.
Thus, the combined strategy developed in this work consists of the automatic application of one method (single-channel 1 or single-channel 2) for each pixel in the satellite images, according to the different values of W and Tb in them. The evaluation of a sample of 24 images showed that the RMSE of the combined method was 1.20 K, and the RMSE of single-channel 1 and single-channel 2 methods were 1.51 K and 2.04 K, respectively. This means that the combined strategy implies an improvement of 20% and 41% regarding each single-channel method (1 and 2). However, the performance of the combined strategy can be worse in the case of W values over 2.5 g/cm 2 . Future work will deepen the analysis of the most extreme situations regarding humidity in order to determine and develop better methods for the computation of LST.

Patents
The developments of this paper have been subjected to Intellectual Property rights registration, in the name of the Universidad de Salamanca (Spain), with reference SA-113-19. The name of the software is TEMISAT-Temperature from Satellite Image.