Harnessing Multi-Source Data and Deep Learning for High-Resolution Land Surface Temperature Gap-Filling Supporting Climate Change Adaptation Activities

1. Introduction

Responding to global warming and adapting to climate change effects such as heat waves and drought is a key priority of stakeholders involved in the definition of Climate Change Adaptation strategies [1]. Most of the largest cities experience profound changes due to urbanization and hence, city administrations are facing challenges in safeguarding high-quality urban growth despite increasingly tight spatial resources. Studies have shown that the degree of surface sealing has a direct impact on the radiation balance – by modifying the surface radiative properties and radiatively active air pollutants [2,3], thus influencing the ratio between sensible and latent heat fluxes. These modifications in the built-up environment make cities (but also smaller municipalities) warmer than their surroundings and more prone to excess heat, leading to an urban heat island effect [4]. Other studies have demonstrated that the temperature within urbanized areas shows a high degree of variability depending on the local urban morphology [5,6,7]. Land use and land cover, along with their changes in general, influence the local climate characteristics [8]. Regional and city administrations make use of this knowledge, aiming to reduce health risks related to climate change and to increase human well-being by implementing heat mitigation measures such as green and blue infrastructure [9]. Altering thermal comfort in urban environments enhances time of exposure to an uncomfortable amount of heat. This can be particularly dangerous for vulnerable individuals, as well as for those performing strenuous physical work in high heat, potentially leading to fatal outcomes [10,11]. Moreover, recent studies have highlighted that climate change in general, and increasing temperatures in particular, pose a significant risk for mountainous areas, affecting alpine communities and their economy (e.g., tourism) [12]. Global and regional warming can further amplify the effect of excess heat [13]. Taking adaptive measures is not only the focus of larger cities but also of any urban area and mountainous regions. Understanding how land use and climate trends impact local climates is essential for decision-makers to develop cost-effective, evidence-based, and consistent solutions for sustainable cities and for rural and mountainous communities.

Land Surface Temperature (LST) is a broadly used proxy to investigate the Surface Urban Heat Island (SUHI) effect [4]. The positive correlation between LST and the degree of surface sealing [14,15] represents the intensity of the SUHI which, identical to the urban heat island effect, is higher in urban environments compared to rural areas [14]. The degree of surface sealing is also related to the Bowen ratio ($\mathsf{\beta }$; ratio between sensible and latent heat fluxes) [16]. Surfaces with higher Bowen ratio ($\mathsf{\beta }>1$) indicate lower soil moisture availability [17], leading to an increase in LST, enhanced heat exchange by convection, and an increase in near-surface air temperatures [18], thus intensifying urban heat [14,15]. Conversely, surfaces with lower Bowen ratio ($\mathsf{\beta }<1$) indicate higher soil moisture availability [17], which decreases LST and favors the evapotranspiration-driven cooling effect, leading to a decrease in near-surface air temperatures [19]. LST is a crucial parameter in numerous fields, including surface energy and water balance, ecology, agriculture, environment, climatology, meteorology, and hydrology [20,21,22], contributing to an overall understanding of Earth’s surface dynamics and the impact of climate trends. Improving our understanding of LST and its interplay with surface sealing, land cover, and meteorological conditions is thus paramount, with a wide range of applications involving Surface Heat Island, urban climate studies [23,24,25], drought monitoring [26], surface soil moisture and evapotranspiration estimation [27,28], and numerical weather prediction, to name a few.

A major advantage of using land surface temperature to investigate SUHI effect is its availability from gridded data (e.g., Earth-observation (EO)-based retrieval), which enables the analysis of local effects depending on the resolution of the satellite sensor. In addition, remote-sensing-based LST dataset has been accepted by the International Geosphere and Biosphere Program as one of the high-priority parameters , and the Global Climate Observing System identified it as an Essential Climate Variable [30]. EO-based LST represents the accumulative radiometric surface temperature of all materials of the surface cover within the sensor’s field of view [31]. Thus, LST estimation from thermal infrared images is complex due to the surface composition, with materials of varying emissivity and geometry [32,33,34]. For example, a densely vegetated area represents the surface temperature of vegetation; whereas in a sparsely vegetated area, the surface temperature includes contributions from vegetation and soil simultaneously [33]. EO-based retrieval of land surface temperature has a long tradition dating back to the 1960s with the launch of the TIROS-II satellite [35,36]. Numerous EO sensors subsequently followed, carried on geostationary and low-Earth-orbit satellites and providing data at coarse spatial resolutions (750 m to 4 km, e.g., GOES, SEVIRI, AVHRR, MODIS, AATSR, VIIRS, Sentinel-3 SLSTR). Medium spatial resolution data (70 to 100 m) is offered through the sensors such as the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), the Landsat Thermal Infrared Sensor (TIRS), and the ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS).

Satellite-based land surface temperature datasets have been extensively used, e.g., for advanced assessments of urban heat island effects – overcoming the challenges of conducting high-resolution air temperatures at similar scale [23], in assessment and mitigation studies elaborating on the spatiotemporal behavior of the urban heat islands [37,38,39,40,41,42,43], and in downscaling studies, providing high-resolution LST maps based on the assumption of a scale-invariant relationship between LST and other influencing parameters (e.g., surface sealing degree, reflectance, spectral indices) [42,44,45,46]. These advances establish land surface temperature as an essential parameter for advanced applications in urban and built-up environments. A particularly interesting approach consists of integrating methods operating at different scales to enhance modelling capabilities in heat assessment studies; namely, the integration of LST datasets with the Computational Fluid Dynamics (CFD)-GIS integrated modelling approach [6,7] – enabling modelling of the SUHI effect at very high spatial resolution, and urban climate models used to downscale and evaluate localized effects of heat or other climate-related parameters [47]. A critical requirement for fully utilizing LST within a multi-resolution and multi-sensor data fusion methodology is access to a comprehensive, high-quality LST dataset. Specifically, developing an inventory of LST data at high spatial and temporal resolution is a timely task. Such an inventory would be an important asset in multi-sensor methodologies at large scales (national, regional), enabling consistent analysis of various built-up environments, from small settlements to large cities.

Over the past decade, deep learning models have become integral to the field of Earth observation, providing cutting-edge tools for the exploration and handling of complex remote sensing datasets. Numerous pivotal reviews highlight the breakthroughs in remote sensing, demonstrating the advantages of deep learning over traditional approaches in applications such as image processing, data fusion, time series analysis, object detection, and land cover mapping [48,49,50,51,52,53]. Most of the developed deep learning methodologies focus on scene classification and image segmentation tasks such as land cover and land use classification, change detection, and urban mapping [54,55,56]. Convolutional Neural Networks (CNNs) have been extensively employed for these tasks due to their ability to learn relationships between neighboring pixels in an image, enabling the extraction and learning of complex spatial features [48,49]. Despite the focus on classification and segmentation tasks, regression-based CNN models utilizing satellite imagery are becoming increasingly prevalent, facilitating the estimation of continuous data in environmental, urban, and agricultural contexts [51]. The varied applications include ice concentration estimation [57], quantification of bacterial bloom in water systems [58], soil moisture prediction [59,60], modelling of water quality indicators [61], and surface and air temperature prediction [53,62,63,64].

In this paper, we develop a Convolutional Neural Network framework for the estimation of land surface temperature from land cover and land change information, cloud-free Sentinel-2 data, and hourly meteorological data. The output of the model are land surface temperature maps at the native resolution of ECOSTRESS, accessible at any arbitrary point in time and location given the availability of meteorological data. The lightweight, regression-based CNN methodology designed to effectively fill gaps in the ECOSTRESS LST observations is the main contribution of this work, providing high-quality land surface temperature data over large spatial regions and with high temporal resolution. This establishes an enhanced, high-quality land surface temperature dataset, crucial for integration and use in sophisticated urban climate modelling approaches or spatial planning activities [6,7,47].

Our paper is organized as follows. In Section 2, we describe the areas of interest and the datasets used in this work. In Section 3, we provide details on the CNN methodology. In Section 4, we present the model results and analyze land surface temperature predictions in spatial and temporal contexts. The discussion and conclusions are provided in Section 5 and Section 6, respectively.

2. Areas of Interest and Datasets

2.1. Areas of Interest

We selected two areas of interest for our study (AREA I and II), covering various land cover classes within the urban, rural, and alpine environment.

AREA I, located at approximately 47.27° N and 11.39° E and covering 19.3 × 16.3 km², includes the city of Innsbruck, Austria, and its surrounding area (see Figure 1). The city region includes the typical heterogeneous urban environments (such as buildings and roads), vegetated areas (parks, lawns), and bodies of water (rivers, lakes). The area outside of the city borders includes rural and agricultural features, as well as the typical alpine features such as mountainous terrain, forests, and meadows. Innsbruck has a humid continental climate with highly variable summers, and the June-July daily mean temperature of 18.9 °C [65].

AREA II, located at approximately 48.21° N and 16.37° E and covering 29.8 × 22.4 km², includes the city of Vienna, Austria, and its surrounding area (see Figure 1). The city region is characterized by the typical urban areas, residential, commercial and industrial areas, as well as vegetated and water areas. The greater region also includes extensive agricultural surfaces and forests. Vienna has a borderline oceanic and humid continental climate with warm summers, and the June-July daily mean temperature of 21.2 °C [65].

The primary application of interest within our study being heat monitoring, we focus on the analysis of land surface temperature trends during summer months. Therefore, we define our time window of interest as the summer months of June, July, and August. For the preparation of datasets described in the following subsections, we focus on the data available from the years 2022 and 2023. The data obtained within such constraints is thus a representation of typical summer conditions in AREA I and AREA II.

2.2. Land Surface Temperature

The primary objective of this study is the development of a neural-network-based model capable of capturing and predicting dense land surface temperature time-series over large spatial regions. Achieving this goal strongly relies on the careful selection of input datasets, a crucial step that underpins the modeling process. In the input data selection, we have been guided by the established theoretical models for land surface temperature, which provide the essential framework for relating the land surface temperature to the various physical parameters characterizing the surface and the surrounding environment.

Land surface temperature is a key environmental metric that measures the radiative skin temperature of the Earth's surface. Unlike air temperature, which is typically measured at weather stations and reflects the atmospheric conditions above ground, land surface temperature encompasses the thermal infrared radiation emitted by the land surface. Consequently, it is strongly dependent on the atmospheric dynamics, global radiation, and the reflective and absorptive surface characteristics [6,66,67,68], with the dependence on these various surface and atmospheric parameters captured by a physical model [66,67]:

$$T_{\mathrm{surface}}=T_{\mathrm{air}}+{\displaystyle \frac{Q+B}{\left(6.2+4.26\times v_{\mathrm{wind}}\right)\left(1+\frac{1}{\beta }\right)}}$$

(1)

Here, $T_{\mathrm{surface}} \left(T_{\mathrm{air}}\right)$ denote the surface (air) temperature, $Q$ represents the net all-wave radiation flux, $\mathrm{B}$ represents the substrate heat flux, $v_{\mathrm{wind}}$ is wind speed, and $\beta$ is the Bowen ratio [16]. The model in Eq. (1) captures the fundamental relationships governing land surface temperature in relation to ambient and surface parameters. However, it is of limited practical use, its utility being constrained by the demand for precise calculations of multiple heat fluxes, a task that is often challenging and spatially limited. Instead, we opt for a more robust approach leveraging a neural-network-based methodology, as described in more detail in Section 3, and we base our selection of input parameters on the physical model given by Eq. (1).

In particular, the parameters $T_{\mathrm{air}}, Q, B, v_{\mathrm{wind}}$ represent variables dependent on ambient meteorological conditions, while the Bowen ratio $\beta$ elucidates the relationship between sensible heat and latent heat, depending on factors such as surface type (e.g., vegetation, water, urban), weather conditions, and time of day [69]. With our choice of input parameters, we carefully address and encompass the elements from the physical model into our neural network model, ensuring a model that is both physically meaningful and computationally robust.

2.3. Input Datasets

The modelling workflow presented in this study makes use of multiple datasets, which are summarized in Table 1 and described in more detail in the following subsections. The thermal data consists of the ECOSTRESS 70 m LST product (Section 2.3.1). The meteorological data is provided by the Integrated Nowcasting through Comprehensive Analysis (INCA) dataset of the Austrian national meteorological service (Section 2.3.2), and the optical multispectral data includes Sentinel-2 imagery (Section 2.3.3). Additional land cover and topographic datasets are described in Section 2.3.4.

2.3.1. ECOSTRESS: Label Dataset

ECOSTRESS (ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station) is an ongoing NASA scientific mission mounted on the International Space Station (ISS). The main instrument in ECOSTRESS is a multispectral thermal infrared radiometer, which collects and provides measurements of the surface temperature [70]. We use the ECOSTRESS Land Surface Temperature and Emissivity (ECOSTRESS-LSTE) Level 2 dataset for the supervised learning approach. The ECOSTRESS-LSTE L2 product is provided on 70 m resolution, with irregular revisit times of one to five days according to the flight pattern of the ISS [71]. Due to the inclined and precessing orbit of the ISS, observation times of the ECOSTRESS instrument vary, with some days providing multiple observations per day [72]. Spatially, the ECOSTRESS observations usually cover a large geographic area; specifically, a single ECOSTRESS-LSTE observation is typically sufficient to reproduce a single scene over our areas of interest. ECOSTRESS-LSTE is the cornerstone dataset of our model, providing labels for the supervised learning of the land surface temperature.

The corresponding ECOSTRESS data acquisition follows via an automated download directly from the NASA search portal [73]. The raw swath data is transformed into gridded single-band images in the Universal Transverse Mercator (UTM) projection using the open-source swath2grid conversion algorithm [74]. Particularly, we extract the land surface temperature (LST) and the quality control (QC) layers. The georeferenced LST layer is passed through several quality assessments and artifact mitigation steps. First, we make use of the intrinsic QC layer to mask out any low-quality pixels. The QC unsigned 16-bit data are stored as bit flags in the layer, with flags related to data quality, cloud, Temperature and Emissivity Separation (TES) algorithm diagnostics, and error estimates [70]. For the purposes of this study, we mask the LST layer by keeping only the best quality, cloud-free pixels, corresponding to the value for QC bits 1&0 = 00. This step is generally sufficient to eliminate most artifacts and clouds appearing in the observations. Second, we perform further manual checks to evaluate (i) the accuracy of georeferencing, and (ii) the quality of the QC masking. (i) In the case of georeferencing offset of more than 50 m, the image is labeled as incorrect and discarded from further analysis. (ii) In the case of insufficient cloud masking (identified by extremely low negative temperatures), custom additional masking is performed.

The LST dataset obtained after filtering suffers minimal issues due to incorrect georeferencing and undetected cloud coverage. As the result of the filtering, the data availability is reduced, leading to aperiodic availability of observations, and a typical 2–5-day revisit period in a specific location. The filtered ECOSTRESS LST data comprises the final dataset supplied to our model as the training reference labels. An example ECOSTRESS LST observation is shown in Figure 2 (a).

2.3.2. INCA

The INCA (Integrated Nowcasting through Comprehensive Analysis) dataset is a temporally detailed meteorological dataset provided by the Central Institute for Meteorology and Geodynamics (ZAMG) in Austria. The INCA data is modeled using various available data sources – station observations, remote sensing data, numerical weather prediction models and a high-resolution terrain model – to produce the analysis of the current state of the near-ground atmosphere [75]. INCA data is provided at a 1 km spatial and an hourly temporal resolution, with historical data available dating back to 2013. The dataset covers five weather-relevant parameters which we consider as the input: air temperature, global radiation, relative humidity, and wind speed in two directions. An example INCA air temperature observation is shown in Figure 2 (b). This comprehensive coverage offers a view of weather conditions at any given location, and it provides many of the essential physical parameters required by the physical model in Eq. (1). Spatially, the INCA dataset primarily covers Austria and its surrounding regions, focusing particularly on the areas with complex terrain, such as the Alpine region [75].

Within our data processing workflow, we select the time window of interest, as specified in Section 2, and acquire the full spatial INCA dataset on an hourly resolution via an automated download directly from the GeoSphere Austria Data Hub [76]. From the raw downloaded data, we extract the observations which are temporally closest to the ECOSTRESS LST observations downloaded and prepared as described in Section 2.3.1, such that each extracted INCA observation corresponds to a single ECOSTRESS LST observation. The extracted images are resampled to the native ECOSTRESS resolution (70 m). Although such resampling does not add new information, it helps to smooth the images and facilitates a more systematic comparison with the ECOSTRESS data.

The INCA data is provided on a much coarser resolution than ECOSTRESS (1 km compared to 70 m resolution). However, in the context of land surface temperature modelling this represents a sufficient input, as meteorological conditions are less spatially heterogeneous than the land surface characteristics. On the other hand, a nearly exact temporal matching is possible, due to the high hourly temporal resolution of the INCA dataset. As such, the INCA dataset meets our goal of relating the land surface temperature to the dynamic meteorological conditions at the time of the LST image acquisition.

2.3.3. Sentinel-2

The high-resolution input data for the land surface temperature mapping consists of the Sentinel-2 satellite data of the European Copernicus Earth observation program [77]. The Sentinel-2 mission provides high-resolution multispectral imagery, with spectral bands ranging from the visible to infrared part of the spectrum. Having a decametric resolution and a 5-day revisit time [77], the Sentinel-2 mission accurately captures the reflectance properties of surfaces around the time window of interest. It thus provides spatial information similar to, and more comprehensive than, the Bowen ratio parameter required in the physical model in Eq. (1). This makes it a crucial input for our model, enhancing its capability to analyze and interpret surface characteristics and providing a link between the land surface temperature and the various land cover classes.

For our data processing workflow, we select the Sentinel-2 images covering the area and time window of interest, and retrieve them locally from the Sentinel-2 Cloud Storage bucket [78]. We consider six spectral bands for the input (B2: Blue, B3: Green, B4: Red, B8: Near Infrared, B11 and B12: shortwave infrared), as well as the Scene Classification Layer (SCL) for the identification of clouds and shadow. An example Sentinel-2 B04 observation is shown in Figure 2 (c). Sentinel-2 imagery is divided into standardized 100 km × 100 km tiles. AREA I considered in our study is covered by the tile 32TPT, while AREA II is covered by the tiles 33UWP and 33UXP, with each tile having the projection of the 32 or 33 UTM zone (UTM/WGS84 projection). As the input for our model, we select and keep only the best available scenes, focusing on (i) mitigating cloud contamination and (ii) reducing data gaps due to swath overpass patterns, as the presence of any such low-quality areas in an image can negatively impact model training. To that end, we remove any of the images for which the metadata information indicates more than 10% cloud coverage. Additionally, we discard all the observations that contain more than 1% data gaps. Such filtering leads to, on average, 2-4 scenes available per tile and per summer.

Even though the filtering described above leads to a modest final number of total observations per tile, for the purposes of our analysis this is an adequate input. More specifically, these filtered images are a good representation of the average surface reflectance properties for the time window around which they were taken. As such, they are deemed sufficient to meet our goal of relating the land surface temperature to the more static surface properties, such as the land cover.

Finally, we match the selected Sentinel-2 observations to the ECOSTRESS LST and the matching INCA observations in several steps. (i) We resample (average) the Sentinel-2 images to the ECOSTRESS native resolution (70 m), (ii) we reproject ECOSTRESS and INCA images to the UTM projections, (iii) we clip ECOSTRESS and INCA images to each available Sentinel-2 tile, and (iv) we match any given ECOSTRESS/INCA tiled observation to the temporally closest Sentinel-2 observation.

2.3.4. Additional Input Datasets

In addition to the Sentinel-2 satellite imagery and the INCA meteorological data, in our approach we utilize several supplementary numerical datasets to enhance the robustness and the accuracy of the model. These include the European digital elevation model (EU-DEM) and derived parameters (slope, aspect), which provide detailed topographical characteristics [79], and several land cover datasets provided by the Copernicus Land Monitoring Service [80]:

(i)

tree cover density, a dataset which provides insights into vegetation distribution,

(ii)

water and wetness index, a dataset which indicates moisture levels and the presence of water bodies,

(iii)

imperviousness, a dataset which highlights the areas covered by artificial surfaces such as roads and buildings.

All the additional datasets are shown in Figure 2 (d–g). Combining these various additional layers with the optical and meteorological data accounts for environmental and land surface characteristics, thereby improving the overall modelling capability of our approach, and providing a link between land surface temperature and surface properties.

3. Methodology

3.1. Convolutional Neural Network Architecture for Pixelwise Regression

In this study we utilize a standard Convolutional Neural Network (CNN) model to perform pixelwise land surface temperature predictions based on the input features described in Section 2.3. Our model is adapted from the CNN configurations presented in [56,81], where a CNN architecture is developed for pixelwise segmentation and classification of Sentinel-2 imagery. Building upon this foundation, we have modified the architecture to address a continuous regression problem, enabling land surface temperature estimation.

The input to the pixelwise CNN consists of input features (described in Sections 2.3.2–2.3.4) spatially limited to a narrow image window (patch) around the pixel of interest. Adapting the approach in [56,81], we select a patch of size 5×5 as the input image to the CNN. At the native resolution of ECOSTRESS (70 m), this corresponds to a window of 350 m × 350 m, incorporating the contextual information surrounding each pixel. This inclusion of surrounding prominent features—such as vegetation, water bodies, agricultural fields, and industrial areas—enhances the model's ability to utilize and leverage contextual data during training.

The model architecture is shown schematically in Figure 3. We consider a network configuration consisting of a sequence of four 2D convolutional layers. Each layer performs a convolution with a 2×2 kernel and a stride of 1×1, using the rectifier activation function between layers. The stability and the performance of the model are enhanced by batch normalization layers incorporated after the second and the fourth convolutional layers [82]. Overfitting is mitigated by a subsequent dropout layer [83], which excludes 10% of the neurons during training. The final regression is done via a fully connected dense layer. The output layer has a size 1, utilizing linear activation to produce an output value that corresponds to the land surface temperature prediction for the central pixel in a patch. The network architecture and the relevant hyperparameters are summarized in Table 2.

The training dataset consists of nearly all summer observations, prepared as described in Section 2.3. One observation per year and per area of interest is excluded from the training dataset and used as the test dataset to evaluate the model’s performance after training. The ECOSTRESS data, detailed in Section 2.3.1, serves as the label data. For each area of interest, the training is carried out over the entire corresponding Sentinel-2 tile, namely, tile 32TPT for AREA I and tile 33UWP for AREA II.

3.2. Preparation of the Input Patches

Several steps are necessary to transform the input and label images into the image patches provided as the input to CNN. In the following we summarize the procedure for the training dataset, with the test dataset preparation being analogous. First, we generate 5×5 patches from all training and label images. For this, a sliding window of size 5×5 is applied, such that each pixel appears in a single patch only. Second, we filter the patches, keeping only those where all pixels in a patch have a valid value. This step eliminates patches with missing data in the ECOSTRESS images [as visible, e.g., in Figure 2 (a)], as well as patches with cloudy or shadowy pixels in the Sentinel-2 images, identified using the Scene Classification Layer (SCL). The input for the CNN is thus a multidimensional array of shape $N_{\mathrm{train}}$×5×5×$N_{\mathrm{features}}$, and the labels are a multidimensional array of shape $N_{\mathrm{train}}$×5×5. Here, $N_{\mathrm{train}}$ denotes the total number of filtered patches across all summer observations considered for the training, and $N_{\mathrm{features}}=15$ is the number of input features.

4. Results

4.1. Model Training

During the training phase, the input data is rescaled in the range [0,1] using the min-max scaling. The $N_{\mathrm{train}}$ input data points are randomly shuffled into the training and validation datasets (used to evaluate the model’s performance during training), with the split ratio being 0.2. The configuration of the CNN and the hyperparameters used during training are outlined in Table 2. Training is performed for 50 epochs, with patches processed in batches of size 512. We utilize Adam stochastic optimization with inverse time decay, starting with an initial learning rate of 0.0001, which decays by a factor of 1 every 1000×($N_{\mathrm{train}}$/512) steps. The model is optimized using the mean absolute error (MAE) as the loss function,

$$\mathrm{MAE}={\displaystyle \frac{1}{N_{\mathrm{train}}}}{\displaystyle \sum }_{i=1}^{N_{\mathrm{train}}}\left|{\widehat{T}}_{\mathrm{surface}}^{\left(i\right)}-T_{\mathrm{surface}}^{\left(i\right)}\right| ,$$

(1)

with ${\widehat{T}}_{\mathrm{surface}}^{\left(i\right)}$ and $T_{\mathrm{surface}}^{\left(i\right)}$ representing the predicted and the true LST value for each patch, respectively.

Once the model training is completed, the LST prediction is carried out over the unseen image in the test dataset. The test data is prepared in the same manner as described in Section 3.2, with the key difference being the use of a 1×1 sliding window to generate 5×5 patches for all non-border pixels in the image. This approach facilitates pixelwise predictions, allowing the creation of output images in a raster format that matches the resolution of the input images.

4.2. Performance Evaluation

We assess the model and its generalization capabilities by generating the LST predictions on the test dataset consisting of the images unseen by the model. More precisely, no pixel patches belonging to the test dataset were provided to the CNN for the training. These predictions are compared with the corresponding ECOSTRESS observations to evaluate the performance.

It is important to highlight a key distinction between the generation of the test and validation datasets. The validation dataset is a subset of filtered pixel patches derived from the train/validation split. In contrast, the test dataset is created by excluding an entire image before the train/validation split. The validation dataset is thus less distinct from the training dataset than the test dataset. As a result, predictive performance of the model on the test data provides a more accurate indication of the model’s generalization ability compared to the validation data. This is particularly relevant for assessing the model’s extrapolation capabilities, where "extrapolation" refers to making predictions on data points that fall outside the range of the training data, such as predicting LST in high-temperature scenarios with the model trained on lower-temperature observations.

For more comprehensive analysis, we assess the model over different training datasets by varying the hour range of the observations to be chosen for the training. Therefore, besides using training datasets with observations spanning the entire day, we also consider training datasets limited to smaller hour ranges at different times of the day, like for examples ranges falling in the morning or in the afternoon. Note that, due to the temporal scarcity of ECOSTRESS data, to prevent underfitting most of the available data (within the chosen hour range) must be used for training.

4.3. Qualitative Analysis of the Predictions

We begin by conducting a visual evaluation of the model predictions. In this qualitative assessment, the results are generally promising, as the model predictions closely resemble the ECOSTRESS measurements (see Figure 4 for examples of predictions over AREA I and AREA II). While some visual discrepancies do occur, they are typically limited to small, isolated patches of pixels, which constitute a minor portion of the overall image and are usually located outside of urban areas. These artifacts may originate from occasional model errors or input data issues (e.g., cloud interference). Despite this, the model predictions effectively capture the spatial features present in the ECOSTRESS observations, often further enhancing spatial detail by increasing contrast — particularly noticeable in features like rivers (cf. ECOSTRESS observations and the corresponding model predictions in Figure 4). The improvement in spatial detail is largely due to the integration of higher-resolution input data, which provides additional fine-scale information.

Additionally, the model exhibits strong ability to fill various gaps in the ECOSTRESS observations, including those caused by cloud cover, image border limitations, and grid-pattern sensor artifacts. These gaps, visible in Figure 4 (a,c), are effectively addressed by the model, as shown in Figure 4 (b,d). A close-up of typical ECOSTRESS grid artifacts is shown in Figure 5 (a), and the corresponding correction by the model, as well as the enhancement of spatial detail, is demonstrated in Figure 5 (b).

Overall, the visual analysis shows that the model not only replicates the spatial features of ECOSTRESS observations, but it also enhances detail and effectively fills various gaps appearing in the observations. These strengths highlight the model's potential for addressing common issues in remote sensing datasets and improving the accuracy of land surface temperature predictions.

4.4. Quantitative Analysis of the Predictions

For the statistical analysis, we consider all the pixelwise LST pairs, consisting of the prediction value and the corresponding ECOSTRESS measurement, within a single Sentinel-2 tile. To ensure valid comparisons, we include only the pairs where both pixels contain valid data. As a result, predictions for pixels lacking valid ECOSTRESS observation data are excluded from the analysis. Each prediction is evaluated based on the prediction error,

$${\mathsf{\Delta }}_i={\widehat{T}}_{\mathrm{surface}}^{\left(i\right)}-T_{\mathrm{surface}}^{\left(i\right)} .$$

(3)

The following metrics are subsequently derived to obtain the statistics and the error distribution within a tile: (i) mean error (ME) and mean absolute error (MAE), (ii) coefficient of determination $r^2$, (iii) minimum (min) and maximum (max), (iv) standard deviation (std), and (v) percentiles.

In our focus on urban applications, we prioritize MAE – a metric used as the loss function during training [see Eq. (2)] – as the key statistic for validating the predictions. While the maximum acceptable MAE varies depending on the specific application, a general guideline is that an MAE should not exceed 3°C, with no more than 25% of the pixelwise absolute error above this threshold. This standard aligns with recent studies investigating various machine learning approaches for estimating land surface temperature [63], as well as a milestone survey referenced therein [68].

Additionally, we compute the $r^2$ score, a scale-independent metric that indicates model fit, ranging from 1 (perfect fit) to 0 (mean fit), with negative values indicating poor fits. While $r^2$ is a valuable metric for evaluating performance, it can be challenging to interpret in isolation. In particular, an acceptable $r^2$ value depends on the explainable data variability and the signal-to-noise characteristics, as well as the intention to model the underlying pattern and avoid overfitting. While there is no universal standard for acceptable $r^2$ values in land surface temperature modeling, remote sensing studies typically consider values above 0.5–0.7 as acceptable, with a preference for higher values up to 0.9–0.95, above which the chance of overfitting increases. Note that there is no direct relationship between MAE and $r^2$; low MAE can coexist with low $r^2$and vice versa. Finally, to ensure a comprehensive analysis, we include minimum, maximum, and percentile values to provide a more complete picture of model performance.

A total of 22 image predictions on test data are computed and evaluated. Table 3 summarizes the statistics for the predictions with the 4 best and 4 worst resulting MAE. Note that the predictions are not all derived from the exact same trained model. Instead, they are derived from different realizations of the model yielded by training on observations spanning different hour ranges (as discussed in Section 4.2). The bottom row in each tile cell displays ensemble statistics, computed for all pixels across all 22 predictions. The ensemble statistics provide a comprehensive assessment of the model's overall predictive performance, as opposed to evaluating individual predictions based on specific datasets or single training instances. It is important to highlight that the ensemble MAE represents the weighted average of the single prediction MAEs based on pixel count. However, this statement does not apply to all statistics, $r^2$and percentiles.

As shown in Table 3, the ensemble MAE for tile 32TPT is 1.93°C, with the 10th and 90th percentiles at -2.59°C and +3.16°C, respectively. The ensemble MAE for tile 33UWP is 1.60°C, with the 10th and 90th percentiles at -2.68°C and +2.27°C, respectively. This means that at least 80% of the pixel errors fall within the ±3°C range, with only a slight deviation for the 90th percentile in 32TPT. The slightly worse performance over the 32TPT tile could reflect the more complex microclimate and surface conditions of the alpine environment in AREA I, which increases variability in the temperature data and leads to a broader error distribution. In contrast, the flatter urbanized conditions in AREA II lead to more stable LST patterns, resulting in a narrower error range. Despite this, the error margins in both locations fall well within acceptable limits. Additionally, the ensemble $r^2$ values of 0.87 for both tiles indicate a robust overall fit and a strong correlation between the model predictions and the ECOSTRESS data, despite some individual predictions having lower $r^2$. Finally, it is important to note that maximum overestimation and underestimation errors at individual pixels can occasionally reach extreme values, with some predicted temperatures falling outside reasonable ranges (refer to the max and min columns in Table 3). Such errors are primarily outliers; they correspond to the artifacts reported in the qualitative analysis in Section 4.3 and have a limited impact on overall predictive performance. They could be mitigated through targeted postprocessing techniques, ensemble learning (i.e., combining predictions from multiple models), or by flagging these data points as invalid.

Figure 6 provides further analysis of the LST distributions for the two prediction ensembles, displayed through histograms and bivariate histograms. Panels (a) and (c) compare the predicted with observed LST distributions for tiles 32TPT and 33UWP, respectively. While the two datasets do not align perfectly, the model generally captures the mayor temperature peaks in terms of both width and height, indicating reasonable agreement. The bivariate histograms in panels (b) and (d) further reveal the relationship between the predicted and the true LST values in 33TPT and 33UWP, respectively. Most data points fall along the diagonal, confirming strong overall agreement, though some outliers are present. These deviations may originate from systemic data issues such as unmasked clouds (thin linear artifacts in the histograms), or from a model tendency to overestimate low LST values in alpine areas [a larger concentration of outliers present in panel (b)].

We considered models trained on different time windows to derive the 22 predictions evaluated in Table 3. In particular, we considered full-day and afternoon models trained on observations spanning, respectively, the entire day and afternoon hours. To assess the impact of the training data time window on the model performance, in Table 4 we compare the ensemble statistics for the predictions which are derived from the full-day model with those derived from afternoon models. While the MAE is slightly lower for the afternoon models, $r^2$ values are similar. This indicates that the full-day model is less accurate than the afternoon models, while they explain data variability equally well. However, the error distribution for the full-day model is significantly more spread, with higher maximum overestimation and underestimation errors. In fact, the full-day model mostly contributes to the lower performance metrics in Table 3. Despite having more data for training, the full-day model thus performs worse overall than the narrower afternoon model. This is likely a sign that the relationship between the considered input features (e.g., air temperature, etc...) and LST changes with the time of the day. The input features include in fact only atmospheric conditions at a single point in time, whereas time-series are likely to determine the LST.

For completeness, in Figure 7 we display the LST distributions for the afternoon model prediction ensembles for tiles 32TPT and 33UWP, which can be directly compared with Figure 6. The histograms in panels (a) and (c) show a similarly reasonable agreement as seen in Figure 6, with the morning-related lower temperature peaks now absent. The bivariate histograms in panels (b) and (d) demonstrate improved correlation and fewer outliers, indicating a stronger alignment between predictions and observations compared to the full-day model prediction ensembles.

Finally, all predictions discussed so far are generated within the same time range as the training data. We also evaluate the model's performance when extrapolating beyond this range, such as generating morning predictions from afternoon models, with the results summarized in Table 5. The accuracy for these extrapolated predictions is significantly lower compared to interpolated ones. This is somewhat expected, given the inherent challenges of extrapolation. Nonetheless, the statistics remain close to the acceptable bounds, with approximately 80% of the pixel errors falling within the ±5°C range. This suggests a robust and consistent relationship between input features and LST, with minimal variations between morning and afternoon predictions. However, for more fine-scale modeling and improved accuracy, such extrapolation is not sufficient.

5. Discussion

5.1. Application of the Data Fusion Approach

In our proposed approach, we employ a CNN model for pixelwise land surface temperature predictions at moderate to high spatial resolution (i.e., 70 m) combining multi- source and multi-resolution input features (i.e., coarse resolution reanalysis data from a meteorological nowcasting prediction system, high-resolution land cover data, high-resolution multi-spectral optical satellite data) to gap-fill and densify remote-sensing-based LST observations. The CNN model is trained based on the ECOSTRESS-LSTE L2 dataset provided via the NASA Earthdata portal. The model makes use of the 5×5 input pixel patches, corresponding to a 350 m × 350 m area at the native resolution of ECOSTRESS, incorporating surrounding contextual information such as vegetation cover and water bodies to enhance prediction accuracy. The advantage of the proposed methodology compared to conventional and purely remote sensing-based techniques is that the model-based approach combines the superior dense temporal sampling of the meteorological model data with the substantially higher spatial detail of the remote-sensing-based LST observations. Having such a model in place allows us to get dense LST estimates at a high spatial resolution to support a variety of applications like supporting local spatial planning authorities in identifying hotspot regions within urban environments, analyzing the context between land use and land cover and associated heat impact, or supporting the development of local, near-real-time heat stress warning systems by exchanging the reanalysis data with numerical weather prediction parameters, to name a few. Having such applications in mind, the prediction accuracy (i.e., observed $r^2$ and MAE) of our model is good and allows to highlight and analyze spatial variability and temperature contrasts between various land cover types.

The advantage of our method compared to other approaches which purely rely on remote-sensing-based LST estimates and data cross-calibration (e.g., [44]) is two-fold. First, due to the strong diurnal variations of LST mainly depending on the meteorological conditions, fusing data from various satellite sensors to gap-fill LST data requires sophisticated inter-calibration techniques to correct for LST differences caused by different overpass times. Thus, using meteorological reanalysis data from a highly resolved numerical weather prediction (NWP) model with similar spatial resolution like coarse resolution satellite-based LST estimates (e.g., Sentinel-3, MODIS, VIIRS) has the advantage of offering the meteorological conditions for any observation time and already incorporating a physical-based treatment of meteorological processes. Second, due to the superior temporal sampling of the NWP model outputs, our trained model can be applied to any time step and is less dependent on satellite flight schedules and revisit times, finally enabling us to create a more representative LST dataset which provides comparability between regions over large scales. Furthermore, the NWP output can also be replaced with the output from Regional Climate Models (RCMs) to analyze and evaluate future climate change impacts. Alternative approaches (e.g., [84]) make use of geostationary satellite imagery (e.g., MSG SEVIRI, GOES) to gap-fill LST data, which offers hourly or even sub-hourly temporal observations. Yet, they typically have a larger Ground Sample Distance (~ 5 km in mid-latitudes) then the above-mentioned polar-orbiting instruments (~ 1 km) and cannot be considered for the RCM approach in the context of climate change adaptation studies.

Other studies have successfully used multi-source data fusion to gap-fill LST by incorporating land surface models [85] and applying deep learning techniques to fuse remote sensing data with in-situ observations [64]. These approaches have demonstrated high spatial and temporal coverage with strong accuracy over specific areas and time windows of interest. Our approach complements these efforts by introducing a lightweight CNN model that is robust across time domains, allowing us to generate LST predictions consistently over a 3-month summer period. This flexibility enables the creation of long-term, high-density LST datasets for any region with available meteorological data, facilitating the development of near-real-time LST monitoring systems that can be seamlessly implemented for real-time applications.

5.2. Limitations of the Deep Learning Model

While the model demonstrates generally good performance, there are several limitations that affect its accuracy and reliability. One of the main issues observed is the occurrence of over- and underestimation errors at individual pixels, as reflected in the minimum and maximum values in the ensemble statistics (see Table 3). Although these errors are infrequent, they contribute to outliers that can affect the model's overall predictive performance. Tracing the exact source of these errors is challenging due to the intricate relationships between input features – such as land cover, meteorological data, and satellite imagery – and land surface temperature, as well as the intrinsic complexities of deep-learning-based models. These errors could originate from the inherent variability in surface conditions (leading to a complex land surface temperature modelling response), from the inconsistencies in the input data (e.g., unmasked clouds, noise, and other artifacts), or from model architecture limitations. While the calculated ensemble MAEs for both tiles remain within acceptable limits, addressing these outliers could further enhance accuracy. Potential approaches that could help improve error management and provide a more robust and interpretable framework include tailored postprocessing techniques, ensemble learning, or hybrid models which incorporate additional physical constraints.

The performance of the model is also impacted when trained on the data spanning an entire day, as demonstrated by the comparison of the full-day and the afternoon ensemble statistics (see Table 4). Specifically, the full-day model tends to exhibit poorer performance compared to models trained on specific time ranges, such as afternoon-only datasets, despite having generally more data available for training. This difference is likely due to the dependency of the LST on the input features changing with the time of the day, as more complex time-series dynamics in such relationship are not considered in our model. Consequently, the full-day model exhibits a more spread error distribution, with higher maximum over- and underestimation errors.

Another limitation stems from the scarcity and uneven quality of ECOSTRESS data available for training. The limited availability of good-quality observations affects the model's ability to generalize, particularly in cases where predictions are made outside the training data time range. For instance, morning predictions from models trained on afternoon data show worse statistics than interpolated predictions, but they remain within acceptable limits. Extrapolation, in general, proves to be more error-prone than interpolation due to the model's inability to fully capture diurnal variations in LST.

5.3. ECOSTRESS Data Quality

We encountered several data quality issues when collecting the ECOSTRESS-LSTE data for both areas of interest. These issues are primarily observed in the form of missing data due to cloud contamination, and artifacts and georeferencing inconsistencies, which can impact the reliability of the data. Applying the associated ECOSTRESS quality mask (provided by the QC layer in the ECOTRESS-LSTE dataset) helps mitigate some of these issues.

The quality issues affecting the ECOSTRESS dataset can be divided into three categories: (i) instrument artefacts, (ii) incorrect georeferencing, and (iii) QC layer insufficiencies. (i) Instrument artifacts are a prominent issue, appearing as stripes (linear or grid-like patterns disrupting data uniformity) and grooves (irregular indentations affecting data consistency), as illustrated in Figure 8 (a). These artefacts are caused by damage to the TIR sensor during pre-launch testing [86]. They are an intrinsic property of the ECOSTRESS dataset, arising from the issues in the sensor itself, and as such, they cannot be eliminated. However, as shown in the visual analysis in Section 4.3, the CNN model is able to correctly identify and mitigate these artifacts effectively (see, e.g., Figure 5). (ii) Another significant issue is inconsistent or incorrect georeferencing, as provided by the ECOSTRESS-GEO geolocation dataset within the swath2grid algorithm (see Section 2.3.1 for details). As illustrated in Figure 8 (a), incorrect georeferencing appears as misalignment between LST data and geographic coordinates, while inconsistent georeferencing refers to variability in spatial accuracy across different observations. This issue is sporadic and unpredictable, and developing an automated correction procedure is beyond the scope of this work. In our processing workflow, we addressed this by reviewing all observations in the training dataset and excluding those with inaccurate georeferencing from further analysis, as the accuracy of the deep learning model depends on the quality of the label data (i.e., ECOSTRESS LST). (iii) The original ECOSTRESS QC mask is intended to filter out unreliable data. However, it is not always sufficient and can introduce further complications such incorrect cloud masking, leading to erroneous temperature readings. Specifically, the QC mask occasionally failed to fully remove cloud cover, leading to unrealistic temperature gradients and misleading thermal readings due to partial cloud obstruction. An example of this issue is illustrated in Figure 8 (b-c), where a cloud is not correctly masked. To address this, we screened the QC data and added additional manual masking in observations with significant cloud masking artifacts.

Despite our efforts to mitigate ECOSTRESS data quality issues, some of the observed prediction artifacts (such as the significant under- or overestimation errors reported in Table 3) may still originate from undetected low-quality data used for the model training. Similar data filtering and quality issues with ECOSTRESS datasets have been reported in previous studies [86,87].

5.4. Future Directions

Planned future directions include efforts to enhance the spatial resolution of LST predictions, specifically by exploring techniques to downscale LST to 10-meter resolution. This approach will explore the scaling effects characteristic to downscaling, reported in previous studies [44,45,46], as well as deep learning based methodologies to overcome these issues. Downscaled LST product provided at high temporal resolution and spatial coverage would provide more detailed insights, particularly in urban environments where fine-scale temperature variations are critical for local climate analysis and planning.

Furthermore, the inclusion of in-situ temperature measurements to validate and complement remote sensing LST data should be considered. Planned ground truth validation activities focus on using data from the black ice monitoring system in Upper Austria, where sensors embedded in roads provide high-precision point measurements of surface temperature. Despite these data being limited to one surface type, they serve as valuable reference points for model validation. By integrating in-situ data, future studies can bridge the gap between modeled predictions and real-world conditions, thus providing more robust datasets for urban heat monitoring. Additionally, model comparisons are planned using satellite-derived datasets from commercial satellite operators like OroraTech [88] and Constellr [89], although alignment between timestamps and resolutions remains a challenge.

Finally, a thorough evaluation of the INCA (Integrated Nowcasting through Comprehensive Analysis) dataset should be conducted, as it is a model-based system. Since it is a crucial dataset in attaining the temporal gap-filling and density of the LST predictions, its performance as a predictor for LST estimation is crucial and could affect overall model accuracy. Ensuring the reliability of input models like INCA is crucial for improving the predictive quality of downstream applications.

6. Conclusions

This study presents a deep learning-based methodology for estimating land surface temperature (LST) using a combination of multi-source and multi-resolution meteorological, land cover, and satellite data, including the ECOSTRESS LST data. The proposed Convolutional Neural Network model demonstrates promising results in generating a gap-filled LST time series over large areas at medium to high spatial resolution (<100 m). We showed that the generated LST predictions lead to at least 80% of the pixel errors falling within an acceptable ±3°C range. Unlike traditional satellite-based techniques, our model leverages high-resolution meteorological data to capture diurnal variations, allowing for more robust LST predictions across different regions and time periods. The proposed methodology is timely and can be applied to a variety of fields including urban planning, climate resilience, and real-time heat stress monitoring, with the goal of supporting spatial planning and climate change adaptation activities.

While the model offers robust performance over extended time periods and large areas, several challenges remain, particularly in terms of error traceability, ECOSTRESS data quality, and limitations related to model-based inputs such as INCA. Future work should focus on addressing these limitations by integrating ground truth measurements for improved validation. Through continued refinement and validation, the approach holds significant potential to support climate adaptation strategies and to improve our understanding of land use and climate interactions.