Online Global Land Surface Temperature Estimation from Landsat

Abstract

This study explores the estimation of land surface temperature (LST) for the globe from Landsat 5, 7 and 8 thermal infrared sensors, using different surface emissivity sources. A single channel algorithm is used for consistency among the estimated LST products, whereas the option of using emissivity from different sources provides flexibility for the algorithm’s implementation to any area of interest. The Google Earth Engine (GEE), an advanced earth science data and analysis platform, allows the estimation of LST products for the globe, covering the time period from 1984 to present. To evaluate the method, the estimated LST products were compared against two reference datasets: (a) LST products derived from ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer), as higher-level products based on the temperature-emissivity separation approach; (b) Landsat LST data that have been independently produced, using different approaches. An overall RMSE (root mean square error) of 1.52 °C was observed and it was confirmed that the accuracy of the LST product is dependent on the emissivity; different emissivity sources provided different LST accuracies, depending on the surface cover. The LST products, for the full Landsat 5, 7 and 8 archives, are estimated “on-the-fly” and are available on-line via a web application.

1. Introduction

Land surface temperature (LST) is an important parameter for environmental studies and enables the monitoring of landscape processes and responses [1], such as the surface energy and water balance [2,3,4,5,6,7]. In order to better observe the changes, in climate for example, there is a need for frequent data acquisition to obtain consistent LST time series [8]. Especially at the regional [9] and global scales, the only way to study the surface temperature is with satellite monitoring [10].

The Landsat 5, 7 and 8 satellites carry thermal infrared radiometers and therefore their data are suitable for LST estimation. During the last decades, there is an increasing demand for satellite-derived surface temperature products, to be used in various studies and applications. Landsat-derived LST is widely used for urban temperature studies and particularly for the study of the urban heat island (UHI) phenomenon [11,12,13]. Recently, LST data from Landsat were used for urban energy flux estimation [3]. Landsat-derived LST is also used for monitoring the forested areas, such as the correlation of LST with tree loss or the detection of changes in tropical forest cover [14,15]. The land cover/use changes affect LST as well, as has been demonstrated in different studies, using different types of satellite data [16].

A large number of studies present different algorithms and approaches for LST retrieval from Landsat [17,18,19,20]. Some of these algorithms have also been implemented in software tools. For example, a plugin for the open-source software QGIS [21] calculates LST from Landsat 5, 7 and 8 imagery [22]. Similar tools have been developed for ERDAS [23] software, for processing Landsat 8 thermal acquisitions [24]. However, such tools require software installation and, most importantly, raw satellite data downloads, which can be extremely time-consuming. It should be noted that the U.S. (United States) Geological Survey (USGS) has announced the release of LST products from the full Landsat archive, but these products are not yet available [25].

The objective of this work is to describe a method implemented in a web application, specifically developed to provide fast and easy access to LST for the globe from the Landsat archive. The approach uses a single channel (SC) algorithm [19,20] to analyse Landsat 5, 7 and 8 thermal infrared observations. It is noted that SC approaches are not as efficient as split-window approaches [20,26], but are necessary for satellites acquiring with only one thermal band, like Landsat 5 and 7. Three different sources of emissivity were considered to allow more accurate LST estimates for different regions and conditions: a global emissivity map derived from ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) data, that refers to the period 2000–2008 and is at 100 m × 100 m spatial resolution [27,28]; the MODIS (Moderate Resolution Imaging Spectroradiometer) daily LST/emissivity product at 1 km × 1 km spatial resolution [29]; and emissivity based on the vegetation fraction estimated from NDVI (Normalized Difference Vegetation Index) [30].

The number of studies that use LST is growing, and satellites may provide accurate and up-to-date LST products. Providing easy access to accurate, global and up-to-date LST data is expected to benefit a large number of studies. In the era of information technology, this ever-growing demand for LST data opens the path for new unexplored ways of distributing big data information through web applications. Using online platforms, like the Google Earth Engine (GEE) [31], with access to a huge catalogue of satellite imagery and planetary-scale analysis capabilities, it is now possible to build web applications that do not require any pre-processing or installation of software. An application of this type, capable of “on-the-fly” LST estimation from Landsat thermal infrared observations, for any given geographical area and date range, has been developed and presented in this study. The LST estimation process is performed in the background on the Google cloud computing servers, with direct access to the GEE satellite data catalogue. Therefore, the application needs only a few seconds to produce Landsat LST with no need for any computational resources from the user.

The following section outlines the theoretical basis of the SC methodology, describes the data used for the LST estimation, including the different emissivity sources, describes the algorithm implementation in GEE and the products accuracy assessment approach. In the following, the effect of using different parameterizations of the SC algorithms are discussed and a quantitative analysis is presented, along with the accuracy assessment results. Finally, the web application interface is presented and followed by the conclusion of this study.

2. Materials and Methods

2.1. LST Algorithm

SC is a commonly-used approach to estimate LST from the Landsat thermal infrared observations [19,20,32,33]. Among Landsat 5, 7 and 8, only Landsat 8 carries two thermal bands [34], therefore the SC approach is used in this study for consistency.

Table 1. Landsat 5, 7 and 8 thermal bands: wavelengths, spatial resolution and time period covered.

	Thermal Band(s) #	Wavelength (μm)	Spatial Resolution (m)	Time Period
Landsat 5	Band 6	10.40–12.50	120 (30) 1	March 1984–May 2012
Landsat 7	Band 6	10.40–12.50	60 (30) 1	April 1999–Present
Landsat 8	Band 10 Band 11	10.60–11.19 11.50–12.51	100 (30) 1	April 2013–Present

¹ Thermal band data are acquired at lower resolution and resampled with cubic convolution at a higher spatial resolution before distribution as products by USGS [35,36].

presents the wavelength range and the spatial resolution of each Landsat thermal band, along with the time period of operation. The SC algorithm needs surface emissivity information from external sources.

Using an SC approach, the LST (T_S) can be calculated from the radiance-at-the-sensor in a single band using the radiative transfer equation in the following form:

$$\mathrm{B}\left(\mathrm{LST}\right)=\frac{{\mathrm{L}}_{\mathrm{sen}}{-\mathrm{L}}_{\mathrm{up}}-\mathsf{\tau }·\left(1-\mathsf{\epsilon }\right)·\frac{{\mathrm{L}}_{\mathrm{down}}}{\mathsf{\pi }}}{\mathsf{\tau }·\mathsf{\epsilon }}$$

(1)

where B is the Planck function, L_sen is the radiance-at-the-sensor, L_up is the thermal path radiance, L_down is the downwelling irradiance, ε is the surface emissivity and τ is the atmospheric transmissivity.

In Equation (1) the L_up, L_down and τ need to be known for the specific area for which we calculate the LST. These variables make the process of establishing a single global LST algorithm extremely difficult, since they are prone to changes, and bound only to the area under study [19,20]. They are also heavily dependent on the total water vapour in an atmospheric column (precipitable water (PW)), which is also easier to determine with the use of satellite data as [37]. Jimenez-Munoz et al. [19,20] developed a methodology that uses coefficient tables $\left(C\right)$ through simulations with many different atmospheric profile inputs, capable of parameterizing the local transmissivity, path radiance and downwelling irradiance, to lead to one equation that is more convenient to use on a global scale LST calculation. In this study, LST is estimated using the approach proposed by Jimenez-Munoz et al. [19,20]:

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

(2)

where:

$$\mathsf{\gamma }={\left\{\frac{{\mathrm{c}}_2·{\mathrm{L}}_{\mathrm{sen}}}{{\mathrm{Tb}}_{}^2}·\left[{(\mathsf{\lambda }}^4·{\mathrm{L}}_{\mathrm{sen}}{)/\mathrm{c}}_1+1/\mathsf{\lambda }\right]\right\}}^{-1}$$

(3)

$$\mathsf{\delta }=-\mathsf{\gamma }·{\mathrm{L}}_{\mathrm{sen}}+\mathrm{Tb}$$

(4)

$$\mathsf{\psi }=\mathbf{C}·\left[\begin{array}{c}{\mathrm{PW}}^2\\ \mathrm{PW}\\ 1\end{array}\right]\Rightarrow \left[\begin{array}{c}{\mathsf{\psi }}_1\\ {\mathsf{\psi }}_2\\ {\mathsf{\psi }}_3\end{array}\right]=\left[\begin{array}{c}{\mathrm{c}}_{11}{\mathrm{c}}_{12}{\mathrm{c}}_{13}\\ {\mathrm{c}}_{12}{\mathrm{c}}_{22}{\mathrm{c}}_{23}\\ {\mathrm{c}}_{31}{\mathrm{c}}_{32}{\mathrm{c}}_{33}\end{array}\right]·\left[\begin{array}{c}{\mathrm{PW}}^2\\ \mathrm{PW}\\ 1\end{array}\right]$$

(5)

where Planck’s constant value${\mathrm{c}}_1$ is 1.19104 × 10⁸ W μm⁴ m⁻² sr⁻¹ and ${\mathrm{c}}_2$ is 14,387.7 μm K; λ is the central wavelength of the thermal band of the Landsat sensor in question; L_sen in W sr⁻¹ m⁻² μm⁻¹; Tb is the brightness temperature in Kelvin; C is the coefficients table, with c_ij derived by simulations using different atmospheric profiles [19,20]; and ψ_x is the coefficients weighted with PW. The various datasets that are used in this study to estimate LST from Equation (2) are described in detail below.

2.2. Data Used for LST Estimation

The main source of data in this study concerns Landsat 5, 7 and 8 acquisitions, provided by the USGS and included in the GEE data catalogue. GEE provides easy and instant access to satellite products and the computing resources to process them directly on the platform, with no need for download. Data from the different sources and satellites are organized in image collections that make their combination and simultaneous processing feasible.

Table 2. List of products in the GEE catalogue used to estimate the Landsat LST.

Data	GEE Product Identifier
Landsat 5, Radiance at sensor from band 6	LANDSAT/LT5_L1T
Landsat 5, Brightness temperature from band 6	LANDSAT/LT5_L1T_TOA_FMASK
Landsat 7, Radiance at sensor from band 6	LANDSAT/LE7_L1T
Landsat 7, Brightness temperature from band 6	LANDSAT/LE7_L1T_TOA_FMASK
Landsat 8, Radiance at sensor from band 10	LANDSAT/LC8_L1T
Landsat 8, Brightness temperature from band 10	LANDSAT/LC8_L1T_TOA_FMASK
MODIS Daily average emissivity from bands 31 and 32	MODIS/MOD11A1
NCEP/NCAR 6-hour temporal resolution of the total column water vapour from a single band	NCEP_RE/surface_wv
ASTER 1 Global image with emissivity from 2000–2008 clear-sky pixels from band 14	NASA/ASTER_GED/AG100_003
Landsat 5, Surface Reflectance product	LANDSAT/LT5_SR
Landsat 7, Surface Reflectance product	LEDAPS/LE7_L1T_SR
Landsat 8, Surface Reflectance Product	LANDSAT/LC8_SR
Fmask, from extra band in GEE’s Brightness temperature products	LANDSAT/LT5_L1T_TOA_FMASK LANDSAT/LE7_L1T_TOA_FMASK LANDSAT/LC8_L1T_TOA_FMASK

presents the data used for LST estimation along with their GEE data catalogue identifier. Short explanations of each product are given below.

2.2.1. Landsat Thermal Radiance-at-Sensor

The level 1T (precision Ortho-corrected) products for each Landsat, provided by the USGS, are orthorectified images of the thermal infrared radiance-at-the-sensor. These products are available in the form of an image collection for each Landsat in the GEE catalogue. The images in the different collections contain digital number values, which are converted to radiance-at-sensor via a GEE function using scaling factors. Although the Landsat thermal bands are of different spatial resolutions (Table 1), there is consistency among different Landsat sensors, since all the products are resampled by the USGS to 30 m × 30 m, using a cubic convolution resampling method [35,36].

2.2.2. Brightness Temperature and Cloud Mask

Although the brightness temperature can be easily estimated by directly inverting the Planck function, the GEE catalogue includes an image collection of Landsat top-of-the-atmosphere (TOA) brightness temperature. The Landsat brightness temperature product in the GEE catalogue contains cloud cover information derived by the Fmask approach [38,39]. Fmask is a method for detecting among others, clouds, cloud shadows and water surfaces in Landsat imagery. The brightness temperature and information on the clouds, cloud shadows and water surfaces was used in this study from these image collections.

2.2.3. Landsat Surface Reflectance

The Landsat Surface Reflectance is a higher-level product of surface reflectance information for six bands (0.25–2.35 μm) at 30 m × 30 m spatial resolution, generated by LEDAPS (Landsat Ecosystem Disturbance Adaptive Processing System) and distributed by the USGS [40]. The Landsat surface reflectance is available through the GEE catalogue, in the form of image collections. The red and near-infrared bands are used for the calculation of the NDVI, which is needed in order to estimate the NDVI-based emissivity in the present study.

2.2.4. Emissivity

Emissivity is a critical variable for the LST estimation, since a small uncertainty in the emissivity (1%) can lead to large errors in the LST (up to 1 K) depending on the setting of the sensor, the climatological conditions and geographical setting of the area [41,42]. In this study, three different sources of emissivity are used: ASTER and MODIS emissivity that are available in the GEE catalogue and NDVI-based emissivity estimated from Landsat red and near-infrared data. The reason for using a variety of emissivity sources is to test their strengths and weaknesses, by comparing their impact on the LST retrieval for different land cover types, corresponding to different landscapes and ecosystems.

Table 3. Emissivity sources presented with their corresponding band number, band wavelength and temporal resolution.

Emissivity Source	Band(s) #	Wavelength (μm)	Temporal Resolution	Spatial Resolution (m)
ASTER	10 11 12 13 14	8.125–8.825 8.475–8.825 8.925–9.275 10.25–10.95 10.95–11.65	1 static image with average pixel values from 2000–2008	90
MODIS	31 32	10.78–11.28 11.77–12.27	Daily, 2000–present	1000
NDVI based Landsat 5	3 1 4 1	0.63–0.69 0.76–0.90	About every 16 days, 1984–2012	30
NDVI based Landsat 7	3 1 4 1	0.63–0.69 0.77–0.90	About every 16 days, 1999–present	30
NDVI based Landsat 8	4 1 5 1	0.636–0.673 0.851–0.879	About every 16 days, 2013–2017	30

¹ These are the red and near-infrared bands from which the NDVI is calculated and then the NDVI-based emissivity. Not to be confused with the MODIS (Moderate Resolution Imaging Spectroradiometer) and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) bands in this table, which directly provide the emissivity values.

summarizes the different emissivity sources used in this study.

The MODIS/Terra Land Surface Temperature and Emissivity products (MOD11A1) provide per-pixel temperature and emissivity values in a sequence of swath-based to grid-based global products. This is a tile-based and gridded level three product produced daily at 1 km × 1 km spatial resolution. The surface emissivities in bands 31 and 32 are estimated from land cover types [43,44]. The average emissivity value of bands 31 and 32 is used in this study. Nearest neighbour resampling is used to match the Landsat spatial resolution.

The ASTER Global Emissivity Database (ASTER-GED) was developed by JPL (Jet Propulsion Laboratory). This product is derived from clear-sky pixels for available ASTER data between 2000 and 2008, covering the globe. It has a (nominal) spatial resolution of 100 m × 100 m and includes, among others, the mean emissivity and standard deviation for all five ASTER thermal infrared bands separately from all the clear pixels during the 2000–2008 period. The emissivity value of band 14 (10.95–11.65 μm) is used in this study. Nearest neighbour resampling is used to match the Landsat spatial resolution.

NDVI-based emissivity is estimated from the Landsat visible and near-infrared bands and typical emissivity values, following the approach described in [30]. The fraction of vegetation cover (FVC) [30] is estimated using Equation (6), by assuming the NDVI threshold for non-vegetated (NDVI_nonveg) and vegetated (NDVI_veg) surfaces to be 0.18 and 0.85, respectively [19]. Finally, the emissivity is estimated using Equation (7), assuming a reference emissivity for non-vegetated (ε_nonveg) and vegetated surfaces (ε_veg) to be 0.97 and 0.99, respectively [19].

$$\mathrm{FVC}={\left[\frac{{\mathrm{NDVI}-\mathrm{NDVI}}_{\mathrm{nonveg}}}{{\mathrm{NDVI}}_{\mathrm{veg}}{-\mathrm{NDVI}}_{\mathrm{nonveg}}}\right]}^2$$

(6)

$${\mathsf{\epsilon }=\epsilon }_{\mathrm{nonveg}}·\left(1-\mathrm{FVC}\right){+\epsilon }_{\mathrm{veg}}·\mathrm{FVC}$$

(7)

The NDVI-based emissivity product is of 30 m × 30 m spatial resolution, matching exactly the Landsat thermal data.

2.2.5. Atmospheric PW Content

Information on PW is necessary for the atmospheric correction of thermal observations. The atmospheric correction is indirectly performed in Equation (2). The NCEP/NCAR (National Centers for Environmental Prediction/National Centre for Atmospheric Research) reanalysis project data have among others, PW values of six-hour temporal resolution (00:00, 06:00, 12:00 and 18:00 UTC) and 2.5-degree spatial resolution [45]. The GEE product comes as an image collection of PW. In order to align the NCEP/NCAR reanalysis PW product with the Landsat imagery in this study, a time adjustment was necessary in order to choose for each day the PW product that is representative for the Landsat acquisition time. This is achieved through the timestamps in the metadata of every image and by estimating the average PW between the two PW values, which are close in time. It is worth noting that the atmospheric PW content can also be estimated with the use of satellite data [37,46]. The MODIS PW product for example could be used instead of the NCEP/NCAR reanalysis, but this product is not yet available in the GEE catalogue.

2.2.6. Coefficient Tables for Atmospheric Parameterization

Different coefficient tables, $\mathrm{C}$ in Equation (5), are obtained through simulations with many different atmospheric sounding inputs [19]. These are created with atmospheric profiles for different locations all over the world. For Landsat 5 and 7, five coefficient tables $\mathrm{C}$ are available, created using five different databases of atmospheric soundings and MODTRAN code [19]: a custom database named STD (six standard) constructed with data from six standard atmospheric profiles (tropical, midlatitude summer, midlatitude winter, subartic summer, subartic winter and U.S. Standard); three databases, TIGR1, TIGR2, TIGR3 (thermodynamic initial guess retrieval), with TIGR1 being a reduced version of the original TIGR1, TIGR2 being an extension of the original TIGR1 and TIGR3 an extension of TIGR2; and SAFREE, a combination of TIGR1, TIGR2 and radiosondes from Meteo France and the Norwegian Meteorological services. For Landsat 8 only one coefficient table was available [20], corresponding to the GAPRI (Global Atmospheric Profiles from Reanalysis Information) database by the ERA-interim (European Reanalysis) reanalysis data [47].

2.3. LST Algorithm Implementation in the GEE

The algorithm described in Section 2.1 is implemented in the GEE using the data described in Section 2.2. GEE has two APIs (application programming interfaces), a JavaScript and a Python interface. Both provide libraries of functions for accessing Google’s computing and storage infrastructure. The JavaScript API is used through ‘playground’, which is the web-based IDE (integrated development environment) of GEE. The ‘playground’ provides a ready-to-use editor to test the code and easily see the results. An issue with the GEE ‘playground’ is that it provides limited access through invitation, and thus it is not suitable for open-access application development. The Python API provides the same functionalities and follows the same principles as the JavaScript API. The main difference, besides the programming language, is that the Python API can be used to develop custom web applications accessible by everyone. The SC method described in Section 2.1 was first implemented in the GEE’s ‘playground’ for easy testing. Once operational, the JavaScript code was translated into Python and a web application was developed using the Google App Engine. The application is accessible to everyone and a user-friendly environment allows easy access to the Landsat LST products.

Figure 1 graphically illustrates the implementation algorithm and the flow of data in the application. Initially, the user draws a polygon on the map to mark the area of interest and provides information on the preferred date range, emissivity source (ASTER, MODIS or NDVI) and Landsat satellite (5, 7 or 8). These choices are made in the front end of the application in an interactive environment. Then an AJAX (Asynchronous JavaScript and XML (extensible markup language)) request containing the user’s preferences is sent to the back end of the application and more specifically to the Python API which provides access to the GEE product catalogue and functions. The selected Landsat radiance and brightness temperature image collections are then accessed and filtered for the required geographical area and date range. For MODIS emissivity, a filter for the selected area and date range is applied, since the MODIS emissivity product is daily. For ASTER emissivity, only the area filter is applied, since there is only one emissivity map corresponding to 2000–2008. For the NDVI-based emissivity, calculations need to be performed. A filter is applied to the Landsat surface reflectance product for the date and area and then calculates the emissivity as described in Section 2.2.4. The different image collections (the radiance, brightness temperature with cloud mask and emissivity) are then joined in a single image collection. The PW data is handled differently than other image collections, since it has a 6-h temporal resolution. First the PW collection is filtered for the selected area and date range. Then the average of the two PW values is selected, the first value is the PW that was acquired closest in time before the Landsat image and second is that acquired closest after the Landsat image. The PW is added using this methodology to every image in the joined collection which now contains all the necessary data for the LST estimation of the specific geographical area and date range. It is worth noting the flexibility of GEE to easily combine all the satellite data products from different sources and of different spatial resolutions, using the nearest neighbour resampling method (or other method if requested). Equation (2) is then applied to every image in the joined collection to estimate the LST products for the selected area and date range. The water surfaces, clouds and shadows are removed from the LST products in the final steps by using the Fmask (Section 2.2.2). In the end, the Python API requests from the GEE a download link and a map layer from Google Maps to visually represent every image. The AJAX request returns to the front end with that information so both the visualization and download functions are available to the users for the calculated Landsat LST products.

2.4. Accuracy Assesment

Since validation with in-situ data was beyond the scope of this study and since at the time of the study other operational Landsat LST products were not available for comparison, the accuracy assessment of the LST products was based on indirect verification. Thus, ASTER higher-level LST products [48] were considered as a reference dataset to assess the accuracy of the Landsat-derived LST. ASTER LST products were selected for their high accuracy, which is reported to be of 1.5 K [49]. A thorough analysis of ASTER radiance and temperature accuracies is reported in [50]. Using instrumented buoys in large lakes, for 246 independent validations of surface temperature, they reported that all channels were well within the preflight specification of ±1 K. Comparing a dynamically changing variable such as the surface temperature derived from different satellites is challenging. In this study, a time difference of maximum half an hour is allowed, therefore, ASTER LST products were compared to Landsat 7 LST products. The LST derived from Landsat 5 and 8 in GEE was evaluated by comparison with other Landsat 5- and 8-derived LST from the same observations, but using alternative approaches, such as the approach implemented by the software ATCOR [51]. GEE was used for selecting ASTER images that meet the reference data selection criteria. Using the GEE catalogue was very efficient in matching the acquisition dates of ASTER and Landsat 7 with time differences of less than half an hour. From the resulting matching Landsat and ASTER image pairs, a number of pairs was manually selected to cover different geographic areas around the world, with distinguishable surface cover (desert, dense forest, cities, etc.). Images corresponding to different seasons and a large range of PW values were selected, to ensure testing for a variety of different atmospheric conditions. For the selected dates, the ASTER LST products were provided by JPL (Jet Propulsion Laboratory). For Landsat 5 and 8, no ASTER data with a time difference of less than half an hour in acquisition time were available. Thus, a set of Landsat 5 and 8 LST products estimated with ATCOR for the cities of Heraklion and Basel were used, that were produced by the DLR (German Aerospace Agency). All the reference images that were used are listed in

Table 4. List of reference images used for accuracy assessment, the area they cover, the date of acquisition, the spatial resolution, satellite source of acquisition and the PW content.

Area of the Scene	Satellite	Acquisition Date	Spatial Resolution (m)	Location (Latitude, Longitude)	Area Size (km2)	PW (g/cm2)
The West Virginia mountains	ASTER	11 October 2004	90	38.6621, −80.4140	2.449	0.72
Amazon forest	ASTER	7 September 2010	90	−0.2032, −57.4076	2.321	4.32
Bangkok	ASTER	11 December 2014	90	13.8167, 100.4589	2.660	2.42
Sahara Desert	ASTER	27 October 2016	90	26.2330, 26.2820	2.346	2.20
Dubai	ASTER	29 June 2009	90	25.2831, 55.3724	1.560	0.83
Tokyo	ASTER	10 May 2015	90	35.6729, 139.7488	2.522	1.10
Basel	Landsat 5	31 July 2010	30	47.5608, 7.5846	415	1.54
Basel	Landsat 8	24 April 2015	30	47.5608, 7.5846	415	0.97
Basel	Landsat 8	23 August 2016	30	47.5608, 7.5846	415	2.92
Heraklion	Landsat 5	19 February 2010	30	35.3398, 25.1330	191	1.06
Heraklion	Landsat 5	16 July 2011	30	35.3398, 25.1330	191	2.02
Heraklion	Landsat 8	29 July 2016	30	35.3398, 25.1330	191	2.60
Heraklion	Landsat 8	7 March 2016	30	35.3398, 25.1330	191	1.36

with their acquisition date, spatial resolution, a reference point of the location, area size and PW values.

The comparison process includes the estimation of the Landsat LST as described in Section 2.1, for 13 images, using all combinations of the three emissivity sources, the five coefficient tables for Landsat 5 and 7 and one coefficient table for Landsat 8. In total, 147 Landsat LST products were available for comparison. Post-processing was necessary for the LST product derived by Landsat 7 in order to be comparable to the ASTER LST products, mainly due to their different spatial resolution. The Landsat LST products at 30 m × 30 m spatial resolution were resampled to 90 m × 90 m, using spatial averaging. The RMSE (root mean square error) was used as an accuracy metric in all cases.

3. Results and Discussion

3.1. Accuracy Assesment

The accuracy of the Landsat LST products was assessed as explained in detail in Section 2.4, based on indirect validation with ASTER and Landsat LST products, derived using different approaches. The set of 147 Landsat LST products was compared to the reference data and the accuracy assessment results for different coefficient tables and emissivity sources are presented and discussed below.

3.1.1. Atmospheric Effects: Coefficient Table Comparisons

The first accuracy test concerns a comparison among LST products estimated based on different coefficient tables. The objective was to assess the magnitude of errors for the different areas tested and to select one model to be used in the GEE implementation. Figure 2 shows the RMSE using different coefficient tables for Landsat 5 and Landsat 7. No results are presented for Landsat 8 since only one coefficient table was available (see Section 2.2.6). For Landsat 5, the use of the SAFREE coefficient table resulted in lower RMSE, with an average value of 1.49 °C (bias −0.54), while TIGR3 was very close with an average RMSE of 1.50 °C (bias −0.64). For Landsat 7, TIGR2 performed better than the other coefficient table with an average RMSE of 1.77 °C (bias 0.9). The reference data had a variety of PW ranging from 0.72 to 4.32 g/cm². Our results are in line with [19], where for a PW range between 0.5 and 3 g/cm² the best RMSE was obtained by using the STD, TIGR2 and SAFREE databases.

Consequently, the coefficient table corresponding to SAFREE was adopted in the final methodology for Landsat 5, the TIGR2 for Landsat 7 and the GAPRI for Landsat 8, as it was the only one available at the time.

3.1.2. Emissivity Effects: Comparisons of the Different Emissivity Retrieval Methods

As already mentioned, three emissivity data sources were used in the LST estimation: the ASTER-GED, the daily MODIS emissivity product and NDVI-based emissivity. This allowed the assessment of the emissivity impact on the LST. The reference data were split into two main categories based on the surface cover of the area: natural/isolated landscapes and urban/peri-urban areas. The first category includes areas that are not much affected by human activity and they do not exhibit vast changes over time. In practise, this category includes images from the Amazon forest, the West Virginia Mountains and the Sahara Desert (Table 4). The second category includes images of cities (Table 4) and the surface cover of the areas depicted in the images are, thus, affected by human activity. Figure 3 shows the RSME of the estimated LST compared to the reference data for each of the two main categories and for each emissivity source. The emissivity of natural/isolated landscapes is not expected to present significant changes over time, since the land cover does not change; in contrast, human activity, particularly urbanization and agricultural activities, are responsible for intense land cover changes, which impact emissivity. As thus expected, the ASTER emissivity provides more accurate LST estimates for the natural/isolated landscapes (average RMSE of 1.54 °C, bias 0.66) than MODIS and NDVI-based emissivity. For the urban/peri-urban areas on the other hand, all emissivity sources provide very similar results in terms of LST error. The reason is that these are highly heterogeneous and mixed areas in space.

Overall, from Figure 3, it seems that the ASTER emissivity provides slightly better results, despite the fact that it is a single static image compared to the other frequently-updated emissivity options. This can be explained by the high quality of ASTER products compared to the other emissivity options, in the areas where the average emissivity remains similar, the values from the static ASTER image are still relevant. On the other hand, in urban/peri-urban areas, where land cover/use changes continuously and alters emissivity’s spatial distribution, the more frequent emissivity retrievals from MODIS and NDVI seem to be on par with ASTER accuracy. For future estimations of LST that are far in time from the 2000–2008 ASTER emissivity map, this gap will probably increase and the ASTER emissivity data will become outdated.

The NDVI-based emissivity maps have a better spatial resolution (30 m × 30 m) than the MODIS-maps (1 km × 1 km), while in terms of temporal resolution both products refer to the Landsat acquisition time. For the NDVI-based emissivity, the NDVI_nonveg and NDVI_veg in Equation (6) have fixed values of 0.18 and 0.85 instead of being estimated particularly for each region. In a global study, the terrain is very much diverse in different areas and these NDVI fixed values are not fully representative for all situations. This possibly explains the slightly poor performance of NDVI-based emissivity compared to the MODIS emissivity in the Landsat LST estimation.

Figure 4a–c shows the Landsat 8 LST maps for the area of Basel, Switzerland on the 24 April 2015, using ASTER- (Figure 4a), MODIS- (Figure 4b) and NDVI-based (Figure 4c) emissivity products. Figure 4d is the Landsat 8 RGB (red green blue) image from the same date. Figure 4e–h presents a zoom in the outlined area of the Basel airport (EuroAirport Basel Mulhouse Freiburg) and Figure 4i–k presents the corresponding emissivity.

The LST maps using different emissivity sources in Figure 4a–c show similar results, as expected from the earlier discussion (see Figure 3b). A closer look at the zoomed area reveals some of the strengths and weaknesses of each different emissivity-derived LST. At the bottom half of the images in the first row, the LST with NDVI-based emissivity (Figure 4c) outlined better the vegetated areas than the LST with ASTER and MODIS emissivity. The NDVI-based emissivity captures well the emissivity of vegetation, which is in the peak in this image since it was acquired during April. This is even more evident in the second row, zoomed-in area in Figure 4. Another noticeable difference in the LST refers to the area of the airport. The ASTER emissivity seems to provide better estimates for the airport area, which is expected to be hotter than the surrounding area. In the third row of Figure 4, the NDVI-based emissivity seems to capture better than ASTER and MODIS the spatial variability of emissivity. The ASTER emissivity product is assumed to be more accurate in the case of non-vegetated surfaces, since it is derived from a combination of long-term emissivity data, while the NDVI-based emissivity assumes a constant emissivity value for non-vegetated surfaces.

3.1.3. Overall Accuracy Assessment of the Landsat LST Products

Following the process of comparison that was described in Section 2.4, the average RMSE for all 147 comparisons was initially computed at 1.89 °C. This can be improved by taking advantage of the findings described in Section 3.1.1, using only the SAFREE database when estimating LST from Landsat 5 and the TIGR2 for Landsat 7. Another improvement can be made by choosing the best emissivity source based on the area as described in Section 3.1.2. Therefore, considering all these adjustments, the final RMSE was found to be 1.40 °C (bias −0.31) for Landsat 5, 1.71 °C (bias 0.77) for Landsat 7, 1.31 °C (bias 0.1) for Landsat 8. Furthermore, an overall RMSE of 1.52 °C was estimated for all Landsat sensors used in this study.

3.2. Landsat LST Web Application

The LST products for Landsat 5, 7 and 8 at 30 m spatial resolution are publicly available through a web application that was developed for this purpose and it is accessible at: http://rslab.gr/downloads_LandsatLST.html. This web application automatically updates with the latest available satellite data, since it is linked to the GEE catalogue, which is updated frequently. The availability of data depends on the operational period of the respective satellites. The Landsat 5 archive extends from March 1984 to May 2012, the Landsat 7 archive extends from May 1999 to present and the Landsat 8 archive from April 2013 to present. LST products can be created for each Landsat sensor for the time period of its operation, shown in Table 1, using NDVI-based emissivity. Although products are available for the same periods using ASTER emissivity, their use is suggested based on the guidelines we provide below. Landsat LST products with the use of MODIS emissivity are available for the period 2002 up to present.

Figure 5 shows a screenshot of the web application. The user can choose the desired date range, Landsat satellite and emissivity source; and draw a polygon to identify an area of interest. Then, by clicking the “Calculate LST” button, the calculations are performed based on the user preferences. The LST products are presented in the form of a map layer overlaid to Google Maps. The LST product(s) corresponding to the date range and area of interest are listed under the map and they are available for download either individually, or as a batch file. Due to the computational limitations, a maximum date range of one month is allowed in the current release of the web application. This will be improved in future releases to allow larger datasets to be handled at once.

Using the ASTER emissivity is recommended for cases of natural/isolated landscapes, with not-significant land cover changes, as was indicated by our analysis; or to use it for LST retrievals from Landsat observations acquired between 2000 and 2008, which is the time period that corresponds to the data acquisition of this ASTER product. The NDVI-based emissivity should be used in cases of intense phenological changes in vegetated areas. The MODIS emissivity should be used in cases of areas with extended homogeneous areas that are changing over time, i.e., extended deciduous forests.

The application was built with the use of HTML (hypertext markup Llanguage), JavaScript, CSS (cascading style sheets), JQUERY for the front end and GEE’s Python API for the back end, as explained in Section 2.3. It provides instant access to global, up-to-date, LST data and it is light-weight in terms of computations. This is achieved with the use of GEE, which enables instant access to the large satellite data archive and transfers all the computational and storage needs to Google’s infrastructure. Moreover, the user does not need to install any software or download any data for processing, as opposed to currently available applications [22,24]. The Landsat LST web application can even run from mobile devices with limited power for a quick snapshot of the LST of an area.

Other sources of data, capable of improving the accuracy of our LST estimations, such as temperature and emissivity retrievals from geostationary platforms like MSG (Meteosat Second Generation) [52], will also be taken into account in our analysis, as soon as they are available in GEE. Feedback by the scientific community on the quality of our LST products is highly welcome and we encourage users that own in-situ observations to publish validation studies.

4. Conclusions

LST information from Landsat is not easily calculated, and so far, it is not directly provided as a standard product. Many studies propose methods and algorithms for retrieving LST from Landsat observations, but only experts can easily implement these approaches. Landsat image downloading and ancillary atmospheric data is also needed and the installation of extra tools and software is required in most cases. This study presents an online application for LST estimation from Landsat 5, 7 and 8 that requires the least possible time and effort to provide maps on-the-fly and without the need for extra resources from the user. For global coverage and the whole Landsat archive, this is only possible using advanced cloud computing technologies, i.e., the GEE. The GEE data catalogue includes all the Landsat imagery and it transfers all the heavy computational processes to cloud computing servers. The user only needs to define the time range and area of interest and receives the Landsat LST maps and products in seconds.

The LST estimation method implements a SC algorithm because it allows LST retrieval with the use of one thermal channel. This allowed for algorithm consistency between Landsat 5, 7 and 8. Since for SC algorithms the emissivity is provided from external sources, the developed web application provides the flexibility to choose between three different emissivity sources: ASTER, MODIS and NDVI-based. An accuracy analysis using ASTER and LST estimations from Landsat using alternative approaches as reference datasets, indicated that the LST can be retrieved from Landsat 5, 7 and 8 within an overall accuracy of 1.52 °C (RMSE). The three different emissivity sources did not significantly differentiate LST estimations; however, recommendations on emissivity source selection were made in this study: the analysis revealed that ASTER-derived emissivity performs better for areas with no significant land cover changes; the NDVI-based emissivity captures phenological changes and their spatial pattern well; and the MODIS-derived emissivity is suitable for extended homogeneous areas.

Further improvement of the web application includes adding extra options for users such as a minimum, maximum and average extractor for every image, increasing the maximum date range, as well as further investigation of the results in order to improve the accuracy of the LST products for specific areas at local scales. There is room for improvement in the calculation of the NDVI-based emissivity, which can provide more accurate estimations with varying definitions of NDVI_nonveg and NDVI_veg in Equation (6), such as the NDVI-THM (NDVI threshold method) approach [53]. Improved emissivity retrieval methods for urban environments [41] are considered for future implementation, to ensure accurate estimation of energy budget components at local scale. The possibility of a consolidated emissivity product is also under consideration. The usability of applications like the one developed in this study are not yet mainstream. Therefore, we appreciate also any feedback by users of the web application, based on their specific needs, in order to improve it and to provide tailored solutions to users in future releases.