Online Global Land Surface Temperature Estimation from Landsat

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.


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 satellitederived surface temperature products, to be used in various studies and applications. Landsatderived 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].

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 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.  1 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].
Using an SC approach, the LST (TS) can be calculated from the radiance-at-the-sensor in a single band using the radiative transfer equation in the following form: where B is the Planck function, Lsen is the radiance-at-the-sensor, Lup is the thermal path radiance, Ldown is the downwelling irradiance, ε is the surface emissivity and τ is the atmospheric transmissivity.
In Equation (1) the Lup, Ldown 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 ( ) 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]: where:  [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.

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 presents the data used for LST estimation along with their GEE data catalogue identifier. Short explanations of each product are given below. 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].

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.

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 nearinfrared bands are used for the calculation of the NDVI, which is needed in order to estimate the NDVI-based emissivity in the present study.

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 summarizes the different emissivity sources used in this study. The MODIS/Terra Land Surface Temperature and Emissivity products (MOD11A1) provide perpixel 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 (NDVInonveg) and vegetated (NDVIveg) 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].
The NDVI-based emissivity product is of 30 m × 30 m spatial resolution, matching exactly the Landsat thermal data.

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.

Coefficient Tables for Atmospheric Parameterization
Different coefficient tables, 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 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].

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. In the back end a function joins the GEE image collections of Landsat and emissivity data, whereas a second function adds the PW (precipitable water) information. The joined collection is then used for the LST estimation. The process is the same in GEE 'playground' as well as in the web application.

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 with their acquisition date, spatial resolution, a reference point of the location, area size and PW values. 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. 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.

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. 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 [19], where for a PW range between 0.5 and 3 g/cm 2 the best RMSE was obtained by using the STD, TIGR2 and SAFREE databases.

Atmospheric Effects: Coefficient
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.

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 MODISmaps (1 km × 1 km), while in terms of temporal resolution both products refer to the Landsat acquisition time. For the NDVI-based emissivity, the NDVInonveg and NDVIveg 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.  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 nonvegetated 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.

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.

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. In this screenshot, the LST product corresponding to the first image in said table is shown overlaid on Google Maps. The missing values correspond to clouds, cloud shadows and water surfaces.
Using the ASTER emissivity is recommended for cases of natural/isolated landscapes, with notsignificant 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.

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 NDVInonveg and NDVIveg 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.