Urban Surface Temperature Time Series Estimation at the Local Scale by Spatial-Spectral Unmixing of Satellite Observations

The study of urban climate requires frequent and accurate monitoring of land surface temperature (LST), at the local scale. Since currently, no space-borne sensor provides frequent thermal infrared imagery at high spatial resolution, the scientific community has focused on synergistic methods for retrieving LST that can be suitable for urban studies. Synergistic methods that combine the spatial structure of visible and near-infrared observations with the more frequent, but low-resolution surface temperature patterns derived by thermal infrared imagery provide excellent means for obtaining frequent LST estimates at the local scale in cities. In this study, a new approach based on spatial-spectral unmixing techniques was developed for improving the spatial resolution of thermal infrared observations and the subsequent LST estimation. The method was applied to an urban area in Crete, Greece, for the time period of one year. The results were evaluated against independent high-resolution LST datasets and found to be very promising, with RMSE less than 2 K in all cases. The developed approach has therefore a high potential to be operationally used in the near future, exploiting the Copernicus Sentinel (2 and 3) observations, to provide high spatio-temporal resolution LST estimates in cities. OPEN ACCESS Remote Sens. 2015, 7 4140


Introduction
The microclimate of cities is strongly influenced by man-made structures and human activities and thus differs from that of the surrounding natural areas, a phenomenon widely known as the urban heat island [1].Besides the urban heat island estimations based on urban canopy temperature differences between the urban core and the peri-urban areas, the respective land surface temperature (LST) differences are commonly used to both characterize the surface urban heat island and support urban energy budget studies [2], since LST is closely related to human comfort and the city energy use [3].LST affects radiation balance, turbulent heat fluxes and evapotranspiration, as well as other key climatic factors in urban environments.Thermal remote sensing is valuable for assessing these effects [2].Satellite thermal infrared (TIR) observations are widely used to assess the extent and the intensity of the surface urban heat island, because of the well-known advantage of satellite spatial coverage (e.g., [4][5][6][7][8]).Lately, there have been attempts to estimate more parameters related to the microclimate of cities, like the urban energy budget fluxes, with the use of thermal remote sensing [9][10][11].Grimmond et al. [12] in a review on the understanding of urban climate, highlighted, among others, the need to explore the use of new measurement techniques, including satellite systems, for the study of urban microclimate.
Despite the long history of available LST products from space, recently, the scientific community expressed its interest in enhancing and improving the scientific value of satellite-derived LST products through the EarthTemp network, an international collaboration focusing on measuring and understanding the surface temperatures of the Earth [13].In an extended report on the EarthTemp open discussion, Merchant et al. [14] list the objectives for progress towards meeting the societal needs for surface temperature understanding, and they identify twenty-eight specific ambitious steps relevant to these objectives.One of these steps explicitly refers to the use of remotely-sensed LST in understanding urban and suburban temperature distributions, because of the geographically complete coverage.At the same time, they identify the limitations of low temporal coverage and viewing angles and characterize the downscaling of satellite TIR imagery for use in urban climatology as a current scientific challenge.Moreover, their recommendations for future missions [14] highlight the need for multi-band TIR imagery, with high spatial resolution, to improve the fundamental knowledge of spatial heterogeneity in surface temperature and surface emissivity at scales better than those resolved by the meteorological-style sensors, which is particularly required for understanding urban areas.
Although LST is routinely derived by satellite TIR observations, currently, there is no spaceborne sensor capable of providing frequent thermal imagery at the spatial resolution needed in urban studies (local scale, or at least 100 m × 100 m).Landsat series of satellites retrieve TIR data at a resolution between 60 and 120 m, with the recently launched Landsat 8 that measures radiation in two thermal bands of 100-m spatial resolution [15].The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is a multispectral thermal instrument that measures in five TIR bands of 90-m spatial resolution [16].Although these sensors provide TIR measurements of high spatial resolution, they have a low revisit rate (twice per month for Landsat 8 and much less for ASTER).On the other hand, coarse spatial resolution instruments, like the Moderate Resolution Imaging Spectroradiometer (MODIS) [17], the Advanced Along-Track Scanning Radiometer (AATSR) [18] or the Advanced Very High Resolution Radiometer (AVHRR) [19], provide daily TIR measurements of 1-km resolution.Sensors aboard geostationary satellites, on the other hand, like Meteosat Second Generation (MSG) [20] or the Geostationary Environmental Satellite (GOES) [21], have sufficient temporal resolution, but poor spatial resolution to support urban applications.ESA's upcoming high spatial resolution mission, Sentinel 2, will have a better revisit rate of twice per week with two satellites, but it will not include a TIR sensor [22]; and TIR observations will be held from Sentinel 3 at 1-km spatial resolution [23].Therefore, current and forthcoming TIR remote sensing is confronting the trade-off between spatial and temporal resolution.
Thus, synergistic algorithms exploiting information from high-resolution visible and near-infrared (VNIR) measurements to disaggregate the coarse-resolution TIR data have drawn the attention of the scientific community.Zhang et al. [24] recently reviewed the existing disaggregation methods for remotely-sensed LST and noted that the decomposition of pixel-based temperatures is a rapidly growing branch of thermal remote sensing.They distinguish two main categories of disaggregation algorithms, the thermal sharpening, referring to the enhancement of the spatial resolution of thermal images, and the thermal unmixing, referring to the decomposition of component temperatures.They conclude that there is a need for a theoretical framework that includes both thermal sharpening and thermal unmixing for improved LST disaggregation.
A number of studies and developed methodologies for enhancing the spatial resolution of LST with a focus on urban areas appears in the literature.Nichol [25] proposed an emissivity modulation technique to enhance the spatial resolution of one nighttime ASTER thermal image using land-cover information classified from higher resolution data.Combining the Normalized Difference Vegetation Index (NDVI) and albedo, Dominguez et al. [26] developed a high-resolution urban thermal sharpener using statistical regressions and achieved 10-m spatial resolution from ASTER thermal data.Essa et al. [27] applied DisTrad (Disaggregation procedure for radiometric surface temperature) [28], a thermal sharpening method originally developed for agricultural areas, to downscale Landsat 7 thermal data to 30 m for an urban site.Those methods, although successful for single acquisitions, do not provide frequent enough observations to monitor the urban climate.Focusing on the diurnal temperature cycle, on the other hand, a number of studies uses data from geostationary satellites and by different means enhances their spatial resolution.Zakšek and Oštir [29] integrated high-resolution data using a diurnal temperature cycle model and produced LST maps of 1-km spatial resolution and 15-min temporal resolution for central Europe, using the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) data.Bechtel et al. [30] downscaled SEVIRI data by a factor of two for monitoring UHI in Hamburg, Germany, while Keramitsoglou et al. [31] tested different regression algorithms for the same purpose over the area of Athens, Greece.Lastly, focusing on the annual temperature cycle, Stathopoulou et al. [32] sharpened AVHRR thermal images and derived daily LST of 250-m spatial resolution for four cities in Greece, while recently, Weng et al. [33] adapted the STARFM (Spatial and Temporal Adaptive Reflectance Fusion Model) [34] for predicting thermal radiance and LST, by fusing MODIS and Landsat data, and produced daily LST maps of 120-m spatial resolution.
In this study, a synergistic method that unmixes the low-resolution TIR measurements using high spatial information from VNIR sensors for estimating high spatial resolution LST is presented.The method includes the estimation of emissivity in high resolution from VNIR data and ancillary information from spectral libraries and applies a spatial-spectral unmixing method for downscaling the thermal radiance measurements.The spatial-spectral unmixing method respects the energy balance requirement, i.e., preserves the radiometric information between scales.The theoretical framework includes both thermal sharpening and thermal unmixing concepts, since LST downscaling is based on component temperatures in combination with thermal sharpening.Moreover, the proposed method is not only suitable for single acquisitions, but it can be applied to estimate LST time series.Given the current and the upcoming satellite thermal sensors, the proposed method can be used for daily LST retrieval at the local scale (hundreds of meters).Frequent, detailed LST estimates over urban areas are essential for urban studies in order to capture the intra-urban temperature variability.
In this paper, the developed method was applied to MODIS TIR and Landsat 8 VNIR imagery for a test site in Crete, Greece.Daily high-resolution LST products were derived for a period of one year.Independent ASTER LST products were used for evaluating the accuracy of the proposed method.The methodology that was developed to enhance the spatial resolution of thermal measurements and to derive time series of high-resolution LST combined with high-resolution emissivity is explained in detail in Section 2. The test site and the data that were used for a test case are presented in Section 3. Results from the application and evaluation of the method to the test site are presented and discussed in Section 4, and conclusions follow in Section 5.

Methodology
A detailed description of the methodology developed for estimating high-resolution LST for urban areas is presented in the following section.The proposed methodology is a multi-step procedure for enhancing the spatial resolution of satellite TIR observations and estimating high-resolution emissivity, with an ultimate goal of deriving LST in high spatial and temporal resolution (Figure 1).The surface cover information is derived from high spatial resolution VNIR observations (tens of meters), and it is then used to improve the low spatial resolution of thermal measurements (hundreds of meters).The terms "high-resolution" and "low-resolution" are used in a relative sense for the given combination of two images, and the annotations (H) and (L), respectively, will be used from now on to distinguish between them.Emissivity is estimated using the high spatial resolution VNIR data and information from spectral libraries.Estimated emissivity and downscaled TIR observations are combined with atmospheric information in a split-window algorithm, and they provide frequent LST estimations at the local scale.

Surface Cover Characterization and Abundances Estimation
High-resolution VNIR is used to characterize the surface cover.The underlying urban landscapes are assumed to be composed of a few fundamental components.Abundances of these components are estimated by resolving the spectral mixture problem.Mapping the urban environment in terms of its physical components preserves the heterogeneity of urban land cover better than traditional classification, characterizes the urban surface cover independent of analyst-imposed definitions and accurately captures changes through time [35].The surface cover characterization in this study is based on the description of Herold and Roberts [36].Vegetation, built-up, bare surface and water bodies are considered the fundamental components of the urban land surface, and at a second level, they are further subdivided based on their material properties and generic surface characteristics.The surface cover is modeled assuming six surface cover types, i.e., green vegetation, non-photosynthetic vegetation, bright man-made, dark man-made, bare soil and bare rock.This scheme is used to structure the diversity of the urban surface cover component, focusing more on the thermal properties.More detailed levels of categories can be differentiated, like for example, different types of vegetation or different types of soils, depending on the site and the available computational power.Representative spectra for the surface cover types are collected from the high-resolution (H) images and used as endmembers in the spectral unmixing process [37].The surface cover abundances ( ) are estimated from the high-resolution VNIR measurements by multiple endmember spectral mixture analysis, following the approach described by Powell et al. [38].

Thermal Unmixing
Spatial-spectral unmixing is used to enhance the spatial resolution of thermal images by taking advantage of the higher spatial resolution surface cover abundances.The approach adopted here is based on the idea described in the multi-resolution image fusion technique [39].Multi-resolution image fusion techniques merge the spatial information from a high-resolution image with the radiometric information from a low-resolution image.Maintaining the original radiometric information in the resulting image is highly important for this kind of application, in order to fulfil the energy balance requirement, meaning that if the resulting image is degraded again to the original low resolution, it should coincide with the original low-resolution image.
The multi-resolution image fusion technique is based on classifying the image of the high-resolution sensor and then retrieving signals of the low-resolution sensor for the classes recognized in the high-resolution data.The algorithm includes the definition of each surface type contributions to the signal of the low-resolution pixels, a window-based unmixing of the low-resolution pixels and the construction of the resulting sharpened image.In this case, instead of using hard image classification, the abundances ( ) of the different urban surface cover types ( ) are used.The contribution ( ) of each surface cover component to a low-resolution thermal pixel is estimated by: where is the number of high-resolution pixels ( ∈ ) corresponding to each low-resolution one ( ).Each low-resolution thermal pixel (L) is then unmixed, using the contextual information of the neighboring pixels (using a moving window of size ): where ( ) is a vector of the thermal radiances of the pixels in the window, ( ) is a matrix of the contributions of surface types to those pixels, represents the thermal radiances corresponding to each surface cover type in the window and is the model residual.To avoid large deviations, a regularization term is introduced to the optimization problem [39,40] to ensure that estimated spectra present small variations: where are the spectra to be defined, is the window size, is the number of surface cover types assumed, ′ ( ) refers to a priori spectra derived solving Equation ( 2) without the regularization term and without using a moving window and is the regularization parameter.For each low-resolution pixel, is estimated by solving the minimization problem described by Equation (3).Finally, the high spatial resolution thermal image is constructed by applying Equation (2) in the high resolution for every low-resolution pixel and the corresponding estimated spectra :

Urban Emissivity and LST
Urban emissivity is estimated following an approach similar to the one described in [37].A representative spectral emissivity value is assigned to each surface cover type, using information from spectral libraries (i.e., the ASTER Spectral Library (SL) [41]).Samples from the library, representative for the cover types, are selected, and their spectral reflectance is convolved with the sensor's spectral response function, to extract reflectance of each thermal band .Assuming Lambertian conditions in the thermal infrared area, emissivity is computed as 1 − .The mean emissivity value of the selected samples for each surface cover is assumed as representative for the surface cover type.The emissivity ( ) for each high-resolution pixel is estimated by: where is the number of surface cover types, is the representative emissivity value for the surface cover type and ( ) are the estimated abundances of surface cover types.
High-resolution brightness temperature for two thermal bands ( , ) can be estimated from high-resolution radiances ( ) arising from Equation (4), by applying the inverted Plank law.Given ( , ) and the respective emissivity products ( , ) from Equation ( 5), the LST is derived in high spatial resolution using a split-window algorithm: where = − , is the atmospheric water vapor content derived from external sources, = ( + )/2 , = − andare the split window coefficients determined from the algorithm calibration [42].All of the above computations in Equation ( 6) are performed in high spatial resolution.

Test Area and Datasets
The proposed method was tested for a case study in the area of Heraklion, which is the largest city in the island of Crete, Greece.It is a rapidly growing urban area and exhibits a mixed land use pattern that includes residential, commercial and industrial surfaces, transportation networks and rural surfaces.Figure 2 shows the Urban Atlas [43] land cover/use map of the study area overlaid on a Landsat 8 true-color composition (02 May 2013).The urban-related surfaces are shown in blue, outlining the city core in the northern part of the study area.Surfaces related to agricultural use are shown in green.The industrial area of Heraklion (marked with A in Figure 2) is located outside the main core in the central eastern part of the city.Another urban settlement, named Archanes (marked with B in Figure 2), is located in the southern part of the study area.Apart from the Heraklion city core (marked with C in Figure 2) and the Archanes settlement, the rest of the study area features a mixed urban and agricultural land use pattern.The methodology described below was applied to derive time series LST for the broader area of Heraklion shown in Figure 2, using a series of MODIS 1-km resolution thermal data and a series of Landsat 8 30-m resolution VNIR data.Daily MODIS Level 1B (MOD021) data from both Terra and Aqua satellites for a period between April 1, 2013, and April 30, 2013, were acquired, along with all Landsat 8 images corresponding to the same time period.The Landsat 8 acquisition dates are shown in Figure 3.Only Landsat images with less than 50% cloud cover in the study site were used in the analysis, and that is the reason for the missing dates that appear in Figure 3. Three high-level products from ASTER (AST08) were used to assess the accuracy of the MODIS/Landsat-derived LST products.The ASTER LST product is provided at 90-m spatial resolution, and the available products dates are shown in Figure 3.The daily MODIS water vapor product (MOD05) was also used to provide ancillary atmospheric information on water vapor and cloud cover.Finally, reference emissivity information extracted from the ASTER SL version 2.0 [41] was used for estimating emissivity.The ASTER SL provides a comprehensive collection of over 2300 spectra of a wide variety of natural and man-made materials covering the wavelength range 0.4-15.4µm.

Figure 3. Graphical representation of Landsat and ASTER data availability (shaded dates refer to cloudy scenes or missing data).
A pre-processing of the Landsat 8 imagery was necessary before applying the downscaling method.The pre-processing includes reprojection of the images and cloud and sea water masking.All of the Landsat imagery was reprojected to match the MODIS Level 1b product projection system [44,45].No manual image co-registration was performed, since a window-based approach is used here that minimizes the co-registration errors between the two sensors.Landsat cloudy pixels and cloud shadows for each scene were masked out using the fmask algorithm [46].Since Heraklion is a coastal city, water is an important component of the scene, but dark pixels, corresponding to sea water, are highly degenerated and cannot be accurately modeled with spectral unmixing techniques.Thus, sea water was treated separately using the Urban Atlas data [43].Finally, MODIS pixels regarded as cloudy in the MODIS water vapor product were masked out of any computations.

Accuracy Assessment
The accuracy of the proposed method was assessed by comparing the resulting high-resolution LST maps to ASTER LST products [16].The ASTER sensor is aboard the Terra satellite, along with a MODIS sensor, and thus, concurrent acquisitions were ideal for comparison.ASTER has a multispectral thermal sensor, and the ASTER LST products are estimated by the temperature emissivity separation (TES) algorithm; their accuracy is estimated to 1 K.The accuracy of the resulting LST products was assessed for three available ASTER scenes (Figure 3) using the following accuracy measures: where is the estimated LST versus the reference .Wilcoxon tests were performed, as well, to check if estimate and reference LST come from the same distribution.

Application to MODIS and Landsat Data
The method described in Section 2 was applied to the series of daily MODIS and the respective Landsat 8 data for a test site in the broader area of Heraklion, Greece.A time series of high-resolution LST (90 m) was derived for the case study.Landsat 8 VNIR and SWIR bands (1-7) were used for estimating the surface cover abundances, as described in Section 2.1.Representative endmembers assigned to the surface cover types were selected for the test site following the endmember selection process described in [37], using higher resolution imagery from Google Earth for visualization purposes.The surface cover type abundances were estimated, by applying multiple endmember spectral mixture analysis [38].The extracted surface cover information was then used to enhance the spatial resolution of MODIS thermal bands (31 and 32), as described in Section 2.2.The window size and the regularization parameter in Equations ( 3) and (4) were set to three and 0.1, respectively, after performing a set of trial tests.Different and values ranging between 3-7 and 0.1-1.0,respectively, were tested and gave very similar results, with the errors slightly increased when increasing the values of and ; thus, the smaller ones were used.Representative emissivity values for each surface cover type were extracted from the ASTER SL [41] following the procedure described in [37], and high-resolution emissivity maps were estimated as described in Section 2.3.High-resolution LST was then derived through the split window algorithm described in Section 2.3 and using low-resolution water vapor information derived from the MOD05 product.MODIS pixels regarded as cloudy were not accounted for in the analysis.Figure 4 shows an example of the methodology application for a cloud-free day (24 April 2013) for which the ASTER LST product was also available.Figure 4a shows radiance values (W m −2 μm −1 sr −1 ) of the MODIS thermal Band 31 (10.78-11.28µm) for the study site in its original spatial resolution (1 km), while Figure 4b shows the respective downscaled values (90 m) estimated following the procedure described in Section 2.2. Figure 4c shows the high-resolution emissivity corresponding to MODIS Band 31 (10.78-11.28µm) and Figure 4d the high-resolution LST, as they were estimated following the procedure described in Section 2.3.Furthermore, the respective MODIS LST product (MOD11) is shown in Figure 4e, whereas the ASTER LST product (AST08) is shown in Figure 4f.

Accuracy Assessment Using Independent LST Products
Since the comparison with in situ LST measurements is extremely difficult in the urban environment, given the inhomogeneity and the three-dimensional structure of the urban surface, ASTER LST products were considered as reference datasets to assess the accuracy of the developed downscaling approach.In the case of Terra, ASTER and MODIS sensors aboard the same satellite, their acquisitions are simultaneous, and thus, their LST products are comparable.Three ASTER LST products at 90-m spatial resolution were available for the time period examined in this study (Figure 3).The reference ASTER LST products were compared to the downscaled MODIS products, and the accuracy measures were estimated in all cases as shown in Table 1.Values of the error metrics for the test site indicate a good agreement between the estimated and the reference LST of less than 2 K in all cases, whereas the p-values from the Wilcoxon tests suggest that both distributions are equivalent.Error distribution boxplots are shown in Figure 5.The red line in the center of the box refers to the median value, while inside the box lie 50% of the observations.The tips outside the box include the rest of the observations, and red crosses outside the box tips are considered outliers.Figure 4d

Daily LST Products for Urban Studies
Hereafter, we present some examples of the use of the daily local-scale LST estimates for urban areas.Figure 6 shows the mean LST (K) estimated for July, 2013, as an example of the method application for time series.The spatial distribution of monthly mean LST was estimated from daily local-scale LST maps by averaging LST in time for every pixel.Cloudy pixels were masked out of every computation.The monthly mean LST has been therefore calculated using the local-scale daily LST estimates for July, 2013.LST hot spots are clearly depicted in the upper-right part of Figure 6, where both the main industrial area and the airport of Heraklion are located.Secondary maxima in the lower-central part of the image correspond to light industrial activity in this area.The LST minimum pattern caused by the green infrastructure around the city wall is also evident in the upper-central part of Figure 6 (located within the urban core).It is therefore obvious that the developed method is capable of providing the intra-urban spatio-temporal LST variability.Furthermore, the combination of such information on urban LST, along with information on the density of urban structures, provides insights into the temperature variations in urban areas and can be useful both for urban planning purposes, as well as for urban climate modeling.Table 2 shows the mean LST values for each Urban Atlas class [43], as calculated using the time series of daily LST estimates at the local scale, shown in Figure 6.In this example, land use classes corresponding to dense urban and industrial areas (i.e., continuous urban fabric and industrial, commercial, public, military and private units) are found to have higher mean LST for July, 2013 (almost 2 K), than the one of rural classes (i.e., agricultural, semi-natural areas and wetlands), which is an indication of urban heat island [1,2].Similar analyses highlight the use of local-scale LST time series to assess the temperature differences in urban and peri-urban areas.Table 2. Mean LST (K) corresponding to July, 2013 (Figure 6), for each class of the Urban Atlas land use data for the test site of Heraklion (Figure 2).  Figure 7 shows the time series of daily differences in local-scale LST means between continuous urban fabric (CUF) and agricultural and semi-natural (ASN) land use/cover classes derived from the Urban Atlas product [43].The daily differences were estimated as: where and are the average LST values corresponding to CUF and ASN, respectively, within the day , excluding missing data corresponding to cloud cover.In the majority of days, these differences are positive, as was expected, being higher than 2 K in several cases.The capability of the developed approach to monitor the intra-urban LST variation at the local scale is clear from Figure 7, revealing its potential to support microclimatic zoning and surface energy budget studies in cities.

Adaptation to Sentinels 2 and 3 Synergy
The application of the downscaling method described in this paper was demonstrated through an example using MODIS and Landsat imagery.Nevertheless, the method was developed independent of the data, and it is suitable for downscaling low-resolution thermal imagery, given higher resolution information on the surface cover [39].It can therefore be applied to Sentinel 2 and Sentinel 3 data, once they are available, to provide frequent LST in the local scale over urban areas.With the launch of the Sentinel 2 satellites, the surface cover characterization is expected to be improved significantly.Sentinel 2 satellites will routinely deliver high-resolution optical images globally, carrying a Multispectral Instrument (MSI), with visible, near-infrared and shortwave infrared sensors comprising 13 spectral bands: four bands at 10-m, six bands at 20-m and three bands at 60-m spatial resolution of 12-bit radiometric resolution (<5% radiometric accuracy) [22].Sentinel 2A and Sentinel 2B will provide frequent acquisitions with a revisit of five days at the Equator, and in combination with the large swath of 290 km, the potential of updating the surface cover information is significantly increased.Moreover, among the Sentinel 3 instruments will be a Sea and Land Surface Temperature Radiometer (SLSTR) measuring in eleven spectral channels, two of which are in the TIR at 1 km, and an Ocean and Land Color Instrument (OLCI) measuring in the VNIR spectrum [23], sharing some common spectral bands with Sentinel 2 MSI.The Sentinels are developed for synergies [47]; thus, algorithms that exploit the common bands of OLCI and MSI may be developed to increase the accuracy of surface characterization and emissivity estimation.

Conclusions
Detailed, frequent and accurate satellite-derived LST estimates are necessary for the study of urban climate.To this aim, a methodology was developed that benefits from higher spatial resolution VNIR observations for emissivity and LST estimation.The method is based on the decomposition of surface component temperatures, but rather than discretizing the surface cover, estimates of the surface cover abundances are used.A spatial-spectral unmixing technique is applied for thermal sharpening, and LST downscaling is based on unmixing the component temperatures.The proposed method can be successfully used for time series LST estimates.
The method was tested for a test site in Heraklion, Greece.Daily LST maps of 90-m spatial resolution were produced, and they were compared to three independent ASTER LST products.The results indicated a good agreement, with RMSE being less than 2 K in all cases, and highlight the potential of the method to support urban climate studies.The analysis made it clear that the developed LST downscaling approach is capable of providing the intra-urban spatio-temporal LST variability, and since it is satellite-based, it is easily transferable to any city for local-scale LST monitoring.
The developed method can be adapted to the Copernicus Sentinels, which will have a high revisit rate to provide high-resolution LST products for urban areas by the synergy of Sentinel 2 and Sentinel 3, exploiting their common spectral channels and their large swaths.The long-term operation of the Sentinels series guarantees the future supply of satellite observations, providing the means for the developed approach to become operational.In this case, the developed method will also become capable of supporting both routine urban planning activities and urban planning activities related to climate change mitigation and adaptation in cities.The synergy between more satellite sensors and other data sources could be exploited in the future to achieve a better surface characterization.The ultimate goal is to assist the study of the urban microclimate and to quantify the thermal pattern of local climate zones [48] in cities using Earth observation data.It is therefore expected to advance the current knowledge of the impacts of the intra-urban LST variability on urban energy budget and, hence, on both urban heat island and energy consumption in cities.

Figure 2 .
Figure 2. The test site, the broader area of the city of Heraklion, Greece.Urban Atlas land use data [43] overlaid on Landsat (02 May 2013) true-color composition.The industrial area of the city is located in A, a small urban settlement (Archanes) is located in B, while the Heraklion city core in C.

Figure 4 .
Figure 4.The high-resolution products (b-d) of the methodology application for 24 April 2013 are presented.The low-resolution ones (a,e), as well as the reference ASTER LST (f) are shown for comparison.
and 4f show the spatial distributions of estimated (d) and the reference (f) LST for 24 April 2013.The spatial pattern of the estimated and reference LST is the same, with lower values in the shaded areas, outside the urban structures.A similar behavior appears in the other two cases of the evaluation (11 June 2013 and 30 August 2013), and thus, the respective images are not shown here.

Figure 5 .
Figure 5. Boxplots of the error distributions, estimated minus reference LST (K) for the three days where reference ASTER LST was available, (a) 24 April 2013, (b) 11 June 2013 and (c) 30 August 2013.

Figure 6 .
Figure 6.Mean LST (K) for July, 2013, calculated using the time series of daily LST estimates at the local scale (90 m × 90 m).

Figure 7 .
Figure 7. Time series of mean daily LST differences between continuous urban fabric (CUF) and agricultural and semi-natural (ASN) land use/cover classes.

Table 1 .
Accuracy measures arose by comparing resulting high-resolution LST (K) to ASTER LST product (MODIS and ASTER are both aboard the Terra satellite).