Evaluating the Reconstructed All-Weather Land Surface Temperature for Urban Heat Island Analysis

Abstract

With the continuous improvement of urbanization levels in the Lhasa area, the urban heat island effect (UHI) has seriously affected the ecological environment of the region. However, the satellite-based thermal infrared land surface temperature (LST), commonly used for UHI research, is affected by cloudy weather, resulting in a lack of continuous spatial and temporal information. In this study, focusing on the Lhasa region, we combine simulated LST data obtained by the Weather Research and Forecasting (WRF) model with remote sensing-based LST data to reconstruct the all-weather LST for March, June, September, and December of 2020 at a resolution of 0.01° while using the Moderate-Resolution Imaging Spectroradiometer (MODIS) LST as a reference (in terms of accuracy). Subsequently, based on the reconstructed LST, an analysis of the UHI was conducted to obtain the spatiotemporal distribution of UHI in the Lhasa region under all-weather LST conditions. The results demonstrate that the reconstructed LST effectively captures the expected spatial distribution characteristics with high accuracy, with an average root mean square error of 2.20 K, an average mean absolute error of 1.51 K, and a correlation coefficient consistently higher than 0.9. Additionally, the heat island effect in the Lhasa region is primarily observed during the spring and winter seasons, with the heat island intensity remaining relatively stable in winter. The results of this study provide a new reference method for the reconstruction of all-weather LST, thereby improving the research accuracy of urban thermal environment from the perspective of foundational data. Additionally, it offers a theoretical basis for the governance of UHI in the Lhasa region.

1. Introduction

Land surface temperature (LST) is an important parameter for studying surface energy balance and land processes [1], and it serves as a fundamental dataset for exploring the surface urban heat island effect (UHI, the following UHI refers to the surface UHI) [2,3]. With the increasing level of urbanization in recent years, UHI has become a significant environmental governance issue, and obtaining high-precision LST is a prerequisite for accurately exploring UHI [4].

The conventional method to obtain the LST involves meteorological station and remote sensing (RS) observation data. However, the meteorological station distribution is relatively sparse, and the LST can exhibit significant variations within short distances [5]. Obtaining the LST through RS observation techniques offers high spatial coverage and easy accessibility, demonstrating a tremendous potential [6]. With the maturation of satellite thermal infrared technology, accessing LST data has become increasingly easier, and the accuracy of LST retrieval has gradually improved [7]. In particular, Moderate-Resolution Imaging Spectroradiometer (MODIS) product data, which have a short revisit period and global coverage, are commonly used in various research fields [8,9]. MODIS LST products have been widely applied in various research fields, including urban thermal environment analyses, agricultural drought assessments, and hydrological model simulations [10,11]. Nevertheless, the presence of clouds and rainfall can affect satellite TIR measurements, resulting in a lack of spatiotemporal continuity in LST imagery. Therefore, this incompleteness diminishes its value as foundational research data, or renders it entirely unusable.

Currently, there are four commonly used methods for LST reconstruction: spatiotemporal information-based methods, passive microwave methods, physical model-based methods and machine learning-based methods. Spatiotemporal information-based methods, including time interpolation algorithms, spatial interpolation algorithms, and interpolation algorithms that utilize both time and spatial information, involve the use of the geographical principle of similarity and proximity to interpolate clear-sky LST pixels that are temporally and spatially adjacent, thus obtaining pixels with clouds [12,13]. Common interpolation methods include inverse distance weighting [14], spline function [15], and kriging interpolation methods [16]. The accuracy of interpolation methods decreases with increasing range of missing pixels, so they are only applicable for reconstructing images with few cloud pixels. The passive microwave method can be used to directly retrieve all-weather LST data by penetrating cloud layers and obtaining surface information [17]. However, the spatial resolution of the LST data obtained by this method is relatively low, and the microwave signal significantly varies with surface characteristics; hence, the data are generally applied as supplementary data for LST generation [7,18]. The physical model method aims to construct mathematical equations based on physical principles and has become an important method for obtaining remote sensing data products [19]. The widely used physical model methods include the surface energy balance method and its derivatives [20,21], which can provide data that are as close to the true values as possible. However, there is often high uncertainty due to the difficulty in determining certain parameters. In recent years, artificial intelligence technology has continuously advanced and has been extensively applied in various fields, including research on LST reconstruction [22,23,24]. This method aims to establish relationships between internal features based on the initial data or observations to discover predictive patterns. LST reconstruction based on this method is also currently a commonly used approach, as it offers an efficient modeling approach that does not require the consideration of specific relationships between variables. Moreover, it has achieved remarkable achievements in terms of accuracy [24,25,26]. The method employed in this study for LST reconstruction involves the utilization of the Random Forest (RF) model, which is a flexible and stable machine learning model and is commonly used for generating various remote sensing products [27,28].

The new generation of the mesoscale WRF model can be used to simulate the surface skin temperature (TSK) that is very close to the LST through real final reanalysis (FNL) data. In fact, some researchers have even directly used the TSK as a substitute for the LST [29,30] Although the TSK derived from WRF simulations lacks the spatial details and accuracy of the satellite remote sensing-based LST (actual observations), it exhibits the characteristic of all-weather availability. This is also the reason why the WRF-simulated TSK (WRF LST) is used in this study. Furthermore, according to related studies, certain multisource remote sensing observation data can affect the variation in the LST, such as the digital elevation model (DEM) [31], normalized difference vegetation index (NDVI) [32], and latitude (LAT) [33]. Based on the above ideas, this paper combines the WRF model with the Random Forest mode, further combining multisource remote sensing observation data to ultimately obtain all-weather LST.

Urbanization causes significant changes in land within urban areas and frequent human activities, leading to the global issue of UHI, which results in the gradual deterioration of the thermal environment and severely affects urban ecological environment and public health [34]. The UHI has become one of the most frequent and influential climate disasters, and studying the UHI accurately by providing all-weather LST has become a hot topic [4]. Examples include the changes in UHI in Hangzhou City [35], the spatiotemporal patterns of UHI in Beijing [36], and the prediction of nocturnal UHI variations in the United Kingdom [37]. This study takes the Lhasa region in China as the study area, which includes Chengguan District and seven counties: Dazi, Duilongdeqing, Qushui, Linzhou, Mozhu-gongka, Nimu, and Dangxiong. Many scholars have also conducted research on the thermal environment in this area and have reached several conclusions. For instance, it is noted that Lhasa is one of the cities experiencing the most intense warming on the Qinghai–Tibet Plateau [38]. Additionally, urban expansion has a significant impact on the increase in average temperature in this region [39]. In recent years, changes in the urban ecology and climate environment (Figure 1d) in Lhasa have attracted much attention. According to related studies [40], the urban thermal environment issues in this area are becoming increasingly serious, requiring urgent research on UHI effects. Based on these remote sensing data and the WRF LST, we reconstruct the all-weather LST using the RF model, with the MODIS LST serving as a reference in terms of accuracy. Additionally, the reconstructed LST is used to analyze the UHI effect in the study area for the months of March, June, September, and December 2020, which is important for exploring the ecological and sustainable changes in the region.

2. Study Area and Datasets

2.1. Study Area

Lhasa is located in the northern part of the middle reaches of the Yarlung Tsangpo River and its tributary, the Lhasa River basin (Figure 1a). Its geographical coordinates range from 29°14′26″–31°03’45″ N to 89°45′9″–92°37′16″ E (Figure 1c). The total area is 29,528.91 km², including Chengguan District and seven counties: Dazi, Duilongdeqing, Qushui, Linzhou, Mozhugongka, Nimu, and Dangxiong. The permanent resident population is 860,000 people. According to the vegetation zoning of Tibet, Lhasa belongs to the shrub zone of the southern Tibetan mountains, with relatively high vegetation coverage (Figure 1b). With economic development and the background of western development, Lhasa has undergone significant changes as the economic, social, and cultural center of Tibet. In 2022, the gross domestic product of the region was 74.757 billion yuan, with the value added of the primary industry at 2.674 billion yuan, the value added of the secondary industry at 29.125 billion yuan, and the value added of the tertiary industry at 42.958 billion yuan. The consumer price index for residents has been increasing year by year.

2.2. Datasets

This paper utilizes five types of data, which can be categorized into MODIS data, FNL data, and other data. The data sources for each type are listed in

Table 1. Information on data sources.

Name	Source	Time period
MODIS LST	MODIS data website (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 8 June 2023)	March, June, September, and December of 2020
FNL	National Centers for Environmental Prediction (https://rda.ucar.edu/, accessed on 8 June 2023)
DEM	Resource and Environment Science Data Center (https://www.resdc.cn/, accessed on 8 June 2023)
LAT	Based on DEM data
NDVI	MODIS data website (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 8 June 2023)

.

2.2.1. Modis Data

This paper utilizes two types of MODIS data. The first type of data are the LST data from MOD11A1 (daytime product), with a spatial resolution of 1 km and a temporal resolution of 1 day. It is used to provide clear-sky pixels. Temperature data are obtained from the scanning product of MOD11_L2. Relevant research indicates that the data exhibit a high level of accuracy, with an error of 2 K relative to the observed data [41]. Due to its high temporal resolution, it is commonly used in various studies related to thermal environments. The second type of data are MODIS NDVI, which is a ground-standard data product derived from MOD13Q1 data product released by the NASA Data Center in the United States. The spatial resolution of the data is 250 m, and the temporal resolution is 16 days.

2.2.2. FNL Data

FNL data are a global meteorological analysis dataset. Due to its integration from multiple data sources and quality control, it exhibits high accuracy and is widely applied in meteorological research and climate analysis [42]. In this article, we utilized the FNL ds083.2 dataset. Although the spatial resolution of the data is coarse at 1°, the data can be dynamically downscaled using the WRF model and have been commonly used as driving data for the WRF model. Additionally, the dataset exhibits a relatively high temporal resolution of 6 h, allowing for the simulation of variable values at specific time points using the WRF model based on these data.

2.2.3. Other Data

The other data in this study include DEM and LAT data. The DEM data are derived from resampled Shuttle Radar Topography Mission data based on the U.S. Space Shuttle [43]. The spatial resolution of these data is 250 m.

3. Methodology

The study in this paper focuses on reconstructing the all-weather LST for March, June, September, and December 2020 in the Lhasa region and exploring UHI effects in the area using the reconstructed data. Because the daily MODIS LST data in these four months were consistently contaminated by clouds, they are suitable for validating the LST reconstruction method proposed in this paper. Additionally, these four months can individually represent the four seasons, making them applicable for conducting seasonal analysis of UHI in the study area. Therefore, these four months have been selected as the research period. The content consists of three parts (Figure 2). The first part involves using the WRF model to simulate the WRF LST source data for LST reconstruction. To ensure spatial consistency between the simulated WRF LST and MODIS LST, the coordinate system is transformed to WGS84 with a resolution of 0.01°. The second part focuses on generating the all-weather LST. Firstly, based on the physical model WRF, WRF LST is simulated. Then, by combining multisource remote sensing observation data such as NDVI, DEM, and LAT, all-weather LST is reconstructed through a Random Forest model. The final part focuses on examining the UHI effect in the study area based on the reconstructed all-weather LST. The urban heat island intensity (UHII) is analyzed to investigate the spatiotemporal distribution characteristics of the UHI. To ensure alignment with the MODIS LST data, the aforementioned five types of data underwent resampling and reprojection operations, resulting in a final spatial resolution of 0.01°, the coordinate system of WGS84 and the temporal resolution is 1 day. These operations were performed in PIE-Basic 6.0 software, which is widely used for preprocessing various sensor images and is a popular professional software option for remote sensing image processing. The LAT data are obtained by extracting latitude information of each pixel from the elevation data. The spatial and temporal resolutions of these data are the same as those of the DEM data.

3.1. All-Weather LST Reconstruction

3.1.1. WRF Configuration

The simulation results of the WRF model exhibit a downscaled spatial resolution through the use of nested configurations. In this study, a three-level nested domain configuration (Figure 1b) was selected. The domains of each level are 127 × 105, 220 × 202, and 334 × 334. They are named D01, D02, and D03, with the spatial resolution of the final level reaching 1 km. The chosen projection type is the Lambert projection, with the center of projection located at the center of Lhasa (30.126° N, 91.016° E). The latitudes of the intersection points where the projection cone and sphere are tangent are 60° and 30°, respectively. The selection of the WRF scheme can affect the accuracy of the simulation results. Based on relevant research conducted in the region [44], the following WRF physics schemes (

Table 2. Parameter configuration for the WRF physics schemes.

Option Name	Parameterization	ID
Cumulus	Kain–Fritsch	3
Planetary Boundary layer	Mellor–Yamada–Janjic	2
Longwave Radiation	RRTMG	4
Shortwave Radiation	RRTMG	4
Surface Layer	Mellor–Yamada–Janjic	2
Microphysics	New Thompson	8
Land Surface	Noah-MP	4

) were selected in this study.

3.1.2. LST Reconstruction

The reconstruction method for the all-weather LST can be divided into two scenarios: the MODIS data are completely cloud-covered data and partially cloud-covered data. The reconstructed LST has a spatial resolution of 1 km and a temporal resolution of 1 day. Under the scenarios of partial cloud coverage, the MODIS LST data can be used as the training set for the RF model, with a training-to-testing pixel ratio of 7:3. Since both the testing set and the cloud-covered pixels are not involved in the training process, the accuracy of the test set pixels is assumed to be the same as that of the cloud pixels. Therefore, the accuracy of the former matches that of the latter in this study. This is the validation method for cloud-covered pixels accuracy in this paper. Under this scenario, the reconstructed LST can be obtained as follows:

$${\mathrm{LST}}_{\mathrm{REC}}=f({\mathrm{LST}}_{\mathrm{WRF}},\mathrm{DEM},\mathrm{NDVI},\mathrm{LAT})$$

(1)

where LST_REC denotes the all-weather LST, LST_WRF, DEM, and NDVI are the driving data for the RF model, and f denotes the RF model.

Under the scenarios with completely cloud-covered MODIS data, resulting in no clear-sky MODIS LST pixels available as a training set for the RF model or too few clear-sky pixels leading to a low accuracy of LST reconstruction (the reconstructed LST ρ < 0.8; ρ greater than 0.8 indicates a high correlation [2]), the reconstructed LST from temporally adjacent images can be used to obtain the reconstructed LST of the current day. Temporal adjacent images do not have to be chosen only one day apart. The selected date is the one closest to the target date, and the LST can be reconstructed using Equation (1). Accuracy verification in this situation is based on a small number of remaining pixels in the image. If all pixels are covered by clouds, the LST for the target date is subjected to accuracy verification. The weighted method based on temporal adjacent images can be expressed as:

$${\mathrm{LST}}_{\mathrm{REC}}={\displaystyle \sum _{i=1}^M{\omega }_i{f}_i({\mathrm{LST}}_{\mathrm{WRF}},\mathrm{DEM},\mathrm{NDVI},\mathrm{LAT})+\rho {\mathrm{LST}}_{\mathrm{cWRF}}}$$

(2)

$$e_i=(\frac{1}{(1-{\mathrm{COR}}_i)\times T_i})/({\displaystyle \sum _{j=1}^M\frac{1}{(1-{\mathrm{COR}}_j)\times T_j}})$$

(3)

$${\omega }_i=e_i(1-\rho )$$

(4)

where LST_REC denotes the all-weather LST using the weighted method of temporal adjacent images, f is the RF model using temporal adjacent images, COR is the correlation between the MODIS LST and LST_REC for the corresponding date, and ρ is the correlation between the MODIS LST and LST_cWRF for the corresponding date. When the number of clear-sky pixels in the MODIS data is 0, ρ is considered as 1. Moreover, e denotes the weight of temporal adjacent images excluding the current day, and ω denotes the weight of temporal adjacent images including the current day. T and M are the number of days between a certain adjacent day and the current day and the total number of days (4 in this study) [45], respectively.

3.1.3. Accuracy Verification Metrics

The accuracy of the reconstructed LST includes the computation of accuracy for both cloud-covered pixels and all pixels. The accuracy of cloud-covered pixels is the error between test set pixels (not involved in model training) and corresponding MODIS LST pixels, while the accuracy of all pixels is the error between reconstructed LST and all clear-sky MODIS LST pixels. Three metrics, namely, the mean absolute error (MAE), the root mean square error (RMSE), and the Pearson correlation coefficient (ρ), were utilized to validate the accuracy of the reconstructed LST data. The MAE is a measure of the average difference between the predicted values and actual observations, while the RMSE reflects the overall deviation between the predicted values and actual observations. In contrast, ρ is a measure of the linear relationship between the predicted values and actual observations.

$$\mathrm{MAE}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N|x_i-y_i|}$$

(5)

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^N{(x_i-y_i)}^2}}{N}}$$

(6)

$$\rho =\frac{{\displaystyle {\sum }_{i=1}^N(x_i-\overline{x})}(y_i-\overline{y})}{\sqrt{{{\displaystyle {\sum }_{i=1}^N(x_i-\overline{x})}}^2}\sqrt{{{\displaystyle {\sum }_{i=1}^N(y_i-\overline{y})}}^2}}$$

(7)

where MAE, RMSE, and ρ are the accuracy verification metrics, N denotes the number of pixels, and x_i and y_i denote the values of the i-th pixel in two images.

3.2. Analysis of the UHI Effect

To quantitatively analyze the UHI effect in the study area, the UHII is adopted to represent the intensity of the UHI effect. To apply this method, it is necessary to first determine the suburban and urban areas. According to relevant research on the boundary between urban and suburban areas in Lhasa [46], this paper uses the LST of Chengguan District and Duilongdeqing as representative of urban LST, and the LST of Dazi, Linzhou, Nimu, Motuo, and Qushui as representative of suburban LST. The difference between these two values reflects the UHII in Lhasa. A value greater than 0 indicates the presence of an UHI effect, while a value less than 0 indicates a cool island effect. The UHII can be calculated as follows:

$$\mathrm{UHII}={\mathrm{LST}}_U-{\mathrm{LST}}_S$$

(8)

where UHII is the urban heat island intensity, LST_U denotes the average LST in the urban area, and LST_S denotes the average LST in the suburban area.

To analyze the spatial distribution characteristics of the UHI effect, in this study, we further classified the UHI effect into different levels. Currently, there are three main methods for defining the levels of the UHI effect: the equal-interval (EI) method, mean-standard deviation (MSD) method, and Moran’s I index–G coefficient (MG) method. Among them, the second method outperforms the other methods in terms of the temperature variation and spatial distribution differences [34]. By using this method to analyze the spatial distribution characteristics of the UHI in the study area, the differences in the spatial distribution of LST can be effectively presented. Therefore, in this study, the MSD method is adopted to classify the levels of the UHI effect (

Table 3. Classification method for UHI levels.

Heat Island Intensity	Temperature Range	Temperature Range
Cold island	Low	TS < μ − σ
Green island	Medium low	μ − σ ≤ TS < μ − 0.5σ
Normal	Medium	μ − 0.5σ ≤ TS ≤ μ + 0.5σ
Subheat island	Medium high	μ + 0.5σ < TS ≤ μ + σ
Strong heat island	High	TS > μ + σ
Cold island	Low	TS < μ − σ
Green island	Medium low	μ − σ ≤ TS < μ − 0.5σ

).

In Table 3, T_S is the LST, μ is the mean of the LST, and σ is the standard deviation of the LST.

4. Results

4.1. All-Weather LST Reconstruction

According to the abovementioned all-weather LST reconstruction method, the LST for March, June, September, and December 2020 was reconstructed. Figure 3 shows a comparison between the MODIS LST and the reconstructed LST (including pixels from the surrounding area of the study region) for these four months, with 19 March, 5 June, and 22 September as dates completely covered by clouds. It can be observed that the MODIS LST is less affected by cloud coverage in December relative to the other months, where extensive cloud coverage is observed in the MODIS LST images. The results indicate that the reconstructed LST effectively captures the spatial distribution characteristics and exhibits high spatial continuity. Moreover, many spatial details, such as the UHI effect, LST differences, and elevation temperature differences, are suitably captured in the reconstructed LST.

4.2. Accuracy of the All-Weather LST

To further analyze the reliability of the reconstructed LST. First, the accuracy of cloud pixels and all pixels for these four months was calculated. Figure 4a reveals that the MAE of cloud pixels is approximately 2 K, the RMSE is approximately 3 K, and ρ consistently exceeds 0.9, indicating a high accuracy in reconstructing cloud pixels. Similarly, as shown in Figure 4b, the MAE of all pixels is approximately 1 K, the RMSE is approximately 2 K, and ρ consistently exceeds 0.9, indicating a high accuracy in reconstructing all pixels.

The average accuracy for each month is calculated above. Next, daily accuracy analysis of these four months was conducted. Figure 5 shows that the accuracy of the reconstructed LST fluctuates with the day in all 4 months. For example, the accuracy is relatively low on 21, 25, and 26 March, 15, 20, and 21 June, 9, 18, and 20 September, and 21, 25, and 30 December. This is mainly attributed to 2 factors: 1. the high cloud coverage in the MODIS LST data, which results in insufficient pixels; and 2. the use of the weighting method of temporal adjacent images to reconstruct the LST. However, the RMSE values are consistently below 3 K, the MAE ranges from 1.0 to 2.0 K, and ρ remains above 0.9 for these dates. On the other dates, the accuracy is higher, with the RMSE ranging from 1.5 to 3.0 K, and ρ still exceeding 0.9. Overall, the reconstructed LST in this study demonstrates a high level of accuracy, with an average RMSE of 2.20 K, an average MAE of 1.51 K, and a ρ value consistently above 0.9. The resulting data possesses both all-weather characteristics and high precision.

4.3. Analysis of the Urban Heat Island Intensity

Figure 6 shows the daily UHII and LST in the urban and suburban areas for the four months. Both the UHI effect and cold island effect occurred in Lhasa in these four months. The urban LST, suburban LST, and UHII values exhibited the highest variability in March, June, and September (Figure 6a–c). The changes in the urban LST, suburban LST, and UHII values remained relatively stable in December (Figure 6d). Furthermore, March and December exhibited the largest number of days with the UHI effect, accounting for 30 days (96.77%) and 25 days (80.65%), respectively, followed by June and September, with 23 days (76.67%) and 18 days (60.00%), respectively. If we consider March, June, September, and December as the four representative months for the spring, summer, autumn, and winter seasons, respectively, the UHI effect in Lhasa mainly occurs during the spring and winter seasons, with the UHII remaining relatively stable in winter.

Figure 7 shows the spatial distribution of the UHI effect in the four months. In March (Figure 7a), the areas and proportions of cold islands, green islands, normal areas, mild heat islands, and strong heat islands are 6829 km² (23.42%), 3980 km² (13.65%), 9230 km² (31.66%), 4251 km² (14.58%), and 4867 km² (16.69%), respectively. In June (Figure 7b), the areas and proportions are 8265 km² (28.35%), 5291 km² (18.15%), 9778 km² (33.54%), 3036 km² (10.41%), and 2787 km² (9.56%), respectively. In September (Figure 7c), the areas and proportions are 9375 km² (32.15%), 5034 km² (17.27%), 8920 km² (30.59%), 3059 km² (10.49%), and 2769 km² (9.50%), respectively. In December (Figure 7d), the areas and proportions are 3513 km² (12.05%), 3273 km² (11.23%), 9422 km² (32.31%), 6170 km² (21.16%), and 6779 km² (23.25%), respectively. The largest heat island area is mostly observed during the spring and winter seasons, followed by the summer and autumn seasons. Additionally, the southwestern part of Lhasa and the Linzhou area in the central region show significant UHI effects. This occurs because the major urban areas of Lhasa are located along the Lhasa River, which is situated in the southwestern region, resulting in a more significant UHI effect in this area relative to the other regions.

5. Discussion

In this study, we combine the WRF model with remote sensing data to reconstruct the all-weather LST at a resolution of 0.01° while using the MODIS LST as a reference in terms of accuracy. This approach overcomes the limitation of incomplete LST retrieval from satellite TIR imagery due to the influences of clouds and precipitation. The reason for using the RF model in this paper is because it has high flexibility and stability, and uses an “averaging” approach to enhance prediction accuracy and control overfitting. It is commonly used in various LST inversion model studies [27,34,47]. The reconstructed LST in this paper exhibited a high degree of accuracy, with an average RMSE of 2.20 K, average MAE of 1.51 K, and ρ greater than 0.9. Reconstructing the all-weather LST has consistently remained a topic of interest [5,48], and several studies have been conducted. For instance, Zhou, et al. [49] proposed the data interpolating empirical orthogonal functions (DINEOF) method for LST reconstruction, with an RMSE of 4.57 K. Duan, Li and Leng [7] fused TIR and passive microwave data to reconstruct the LST, with an RMSE ranging from 3.50 to 4.40 K. Li, et al. [50] introduced a three-step mixed gap-filling method, with an RMSE of 4.78 K. Long, Yan, Bai, Zhang and Shi [9] utilized a data fusion approach combined with bias correction to reconstruct the LST, with an RMSE of 2.72 to 4.00 K. Compared with previous results [41], the reconstructed LST in this study exhibits a higher precision and hence can provide a valuable reconstruction approach for future research and serve as a fundamental data source for studies related to the thermal environment.

Based on the reconstructed all-weather LST, in this study, we also analyzed the UHI effect in the Lhasa region. The occurrence of the UHI effect could be better captured in March and December, indicating that the UHI effect in Lhasa can mainly be observed during the spring and winter seasons. Moreover, the UHII in winter is relatively stable, which is similar to the conclusions drawn by Ciren, Ga, Bu and Ciren [46] based on meteorological station data. Additionally, from the spatial distribution of the UHI effect in Lhasa, high-UHI areas are more prevalent during the spring and winter seasons than during the summer and autumn seasons. Furthermore, the UHI effect is primarily concentrated in the southwestern part of Lhasa and the Linzhou area in the central part. This occurs because the important urban areas of Lhasa are developed along the Lhasa River [38], which is located in the southwestern region of Lhasa, leading to a more pronounced UHI effect relative to the other regions. In recent years, the Lhasa region has experienced rapid economic development, leading to increased density of buildings. This has resulted in strong diffuse reflection and low wind speed, which hinders heat dissipation [40]. Additionally, with the growing population, artificial heat release in the main residential areas has become significant [51]. As a result, the region has exhibited noticeable HUI. The UHI in the Lhasa region reduces air density and lowers atmospheric pressure in the high-temperature areas. This leads to the continuous influx of various exhaust gases and harmful chemicals into the high-temperature zone, making residents highly susceptible to digestive or neurological disorders. To verify the reliability of the UHI results, we selected dates with minimal cloud and precipitation influences in the MODIS data, which provided a higher coverage of clear-sky pixels in both the urban and suburban areas (needed for UHII calculation). A comparison was made between the urban MODIS LST and reconstructed LST, the suburban MODIS LST and reconstructed LST, and the urban UHII and suburban UHII. The results for both datasets were highly consistent (Figure 8), indicating the high reliability of UHII analysis in this study.

As we utilized the WRF model in this study to reconstruct the all-weather LST, its application requires high computational power and a significant amount of time. Approximately three days are needed to simulate one month of data, resulting in relatively high computational costs [52]. With future advancements in computing power, it will be possible to apply the methodology proposed in this paper to reconstruct long-term LST data for various studies on thermal environments. Although the spatial distribution details of WRF LST are not sufficient and its accuracy is relatively lower compared to MODIS LST, its all-weather characteristics make it suitable for compensating for missing data from MODIS LST [13]. Since this paper is based on the reconstruction of LST from clear-sky pixels, the accuracy of cloud pixels is also validated based on clear-sky pixels. In the future, the accuracy of cloud pixels can be further improved by collecting a large amount of station data in cloud-covered areas to enhance the reliability of validating cloud pixels’ accuracy. This study only utilized some commonly used multisource remote sensing data as driving data for RF. In the future, it can be considered to incorporate more driving data that may affect LST, including water-related data, in order to further improve the accuracy of LST reconstruction. Additionally, we only conducted a simple analysis of the UHII in the research area. In future research, social and natural factors can be incorporated to analyze the driving forces behind the UHII. Thus, the impact of human activities on the UHI effect can be comprehensively explored.

6. Conclusions

In this study, we first combine WRF simulation with remote sensing data to reconstruct an all-weather LST dataset for the study area of Lhasa. This approach addresses the limitation of satellite TIR imagery, which is affected by cloud cover and precipitation, resulting in incomplete LST data. Subsequently, the reconstructed LST is analyzed to assess the UHI effect in the region. The method of reconstructing all-weather LST presented in this paper can provide new ideas for obtaining all-weather LST in the future. Although the urban heat island effect was only briefly analyzed in this study, we plan to reconstruct long-term LST in the future to further explore the thermal environment issues in the study area. The research findings are as follows:

(1)

The all-weather LST reconstructed in this study can effectively compensate for the missing pixels in MODIS LST and exhibits high spatiotemporal continuity. Many spatial distribution details are accurately restored. This provides a new method for the reconstruction of all-weather LST.

(2)

The reconstructed LST in this study demonstrates high accuracy, with an average RMSE of 2.20 K, average MAE of 1.51 K, and ρ greater than 0.9. This provides fundamental research data for studies related to the thermal environment.

(3)

In terms of the seasonal analysis, the UHI effect is most frequently observed during the spring and winter seasons in the Lhasa region, with the winter season showing a relatively stable UHII value. In terms of the spatial distribution, the areas with the highest proportion of strong heat islands are primarily observed during the spring and winter seasons, followed by the summer and autumn seasons.