Cloud-Free Land Surface Temperature Reconstructions Based on MODIS Measurements and Numerical Simulations for Characterizing Surface Urban Heat Islands

Land surface temperature (LST) data in the thermal infrared (TIR) band measured by the moderate-resolution imaging spectroradiometer (MODIS) instrument are critical for studying surface urban heat islands (SUHIs); however, these acquired TIR LST data are contaminated by clouds, so it is crucial to develop a method to generate cloud-free LST products. In this article, employing Tianjin as the research area, we combined the weather research and forecasting model with a random forest and a spatial optimization algorithm to propose a cloud-free MODIS-like model (WRFFM). The model can reconstruct cloud-free MODIS-like LSTs and SUHIs are studied. The spatial patterns of the WRFFM LSTs and the MODIS LSTs are consistent; the correlation coefficients in July and December range from 0.8 to 0.91 and 0.8 to 0.93, respectively, and the root mean square errors range from 0.5 to 3.8 K and 0.4 to 1.8 K, respectively, indicating that the modeled results are accurate. We use these WRFFM LSTs to study SUHIs and evaluate the deviations between the MODIS SUHIs and WRFFM SUHIs. When the proportion of clear-sky pixels is below 30%, the deviation is above 3 K, and when the proportion of clear-sky pixels is above 80%, the deviation is below 0.6 K. The results indicate that the developed model can be applied to improve the study of SUHIs and that the number of clear-sky pixels for a city is an important factor that affects the bias relative to the actual SUHI .


I. INTRODUCTION
W ITH the rapid growth of populations and economies, impervious surfaces continue to replace natural surfaces, thus leading to large-scale urbanization [1]. From 2001 to 2018, the global urban area expanded by 1.68 times, reaching a total area of 802 233 km 2 [2]. The thermal properties of impervious surfaces and ground radiation can change significantly in response to the enhanced heat storage and thermal conductivity of urban building materials, changes in urban canopy structures, and the release of heat by human activities [3], [4]. The heat balance differences between urban and suburban areas lead to higher surface and atmospheric temperatures in urban areas than in the surrounding suburbs, resulting in a phenomenon known as "urban heat sources" [5] or "urban heat islands" (UHIs). Changes in urban thermal environments can lead to extremely high temperatures, increase climate risks, affect the human quality of life and well-being in urban areas [2], [6], [7], [8], [9], affect the quality of the ecological environment, and alter ecological functions [10], [11]. Therefore, long-term monitoring of UHIs is necessary to understand their characteristics and the urban climatology that underlies them.
Urban land surface temperature (LST) is an essential quantitative indicator of UHIs [12]. In remote sensing studies, LST is often used to define the surface urban heat island (SUHI) effect [13]. In recent years, LSTs retrieved from satellite thermal infrared (TIR) sensors have strongly promoted SUHI research [14], [15], [16], [17], [18]. Although remote sensing data cannot fully represent surface temperatures, these two metrics are related [19], [20], and SUHI data can be used to clearly distinguish between the heat island effects measured by surface models and those measured by air temperature models [13]. However, in the remote sensing data retrieval process, the TIR signals from the surface cannot penetrate the cloud layer; thus, when a pixel is contaminated by clouds, the resulting satellite observation information may be mixed with or even completely obscured by This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ clouds, resulting in missing values in the obtained data [21]. It is difficult to obtain spatiotemporally continuous LST data [22]. Li et al. [23] found that in the contiguous urban areas of the USA, the valid data in moderate-resolution imaging spectroradiometer (MODIS) LST products accounted for 28% of the total data. Thus, the severe absence of valid LST products limits their scientific application potential.
The retrieval accuracies of remote sensing are high, their data sources and coverages are wide, and their periods are relatively stable [24]. Therefore, remote sensing data have great advantages for studying the spatiotemporal heterogeneity of surface physical processes. However, the indirect detection of ground objects by satellites is easily affected by cloud contamination, and the spatial distribution of the obtained data is poorly continuous. Therefore, the LSTs obtained by remote sensing depend heavily on the utilized cloud-screening algorithm [25]. For MODIS LST product data, an important cloud-masking algorithm has been established that uses a series of visible and infrared thresholds to determine whether the observed region of the Earth's surface was obscured by clouds. If a pixel is covered completely or partly by clouds, its LST is not available (http://modis-atmos.gsfc.nasa.gov/). For this reason, many scholars have used a variety of reconstruction techniques to recover missing LST information resulting from cloud coverage. These methods can be divided into five categories [26]. One category is space-based gap-filling methods; LST has good spatial correlation and continuous spatial distribution characteristics. Through some spatial interpolation methods [27], [28], and geostatistics to fill in cloud-covered pixels, space-based gapfilling methods can be used to reconstruct surface temperature; however, as these methods consider only spatial information, the accuracy is unsatisfactory in cases of concentrated or large clouds. The second category is time-based gap-filling methods. These methods fill cloud-contaminated pixels with the data from valid pixels at adjacent times; examples include linear time interpolation methods [29], time-series and analysis methods [30], and harmonic fitting methods [31], [32], [33]. These methods have high reconstruction performance and are suitable for large cloud coverage areas, but they ignore adjacent geographic information and are very sensitive to land cover changes and sudden changes in natural disaster conditions [34]; the filling accuracy is susceptible. The third category is gap-filling methods based on spatiotemporal information. These methods integrate the spatiotemporal information of TIR remote sensing LST, and can offer the advantages of both time and space gap-filling methods. For example, some spatiotemporal fusion techniques and machine learning models have been introduced to estimate continuous LST in cloudy conditions [35], [36]. The accuracy of these methods mainly depends on the quality of the LST data, and they are suitable for filling large-scale data, but they are still not fully utilized with regard to the available spatiotemporal information. The fourth category is multisource fusion-based gap-filling methods, such as passive microwave measurement (PMW). Compared with TIR remote sensing measurements, these methods entail wider orbital scanning gaps and lower spatial resolutions [37], which create limitations in explaining the spatiotemporal variation in microwave emissivity [38]. Wu et al. [39] proposed to generate gapless all-weather LST by combining PMW and TIR measurements, but this method is more costly than traditional methods and has limited applicability. The fifth category is surface energy balance based (SEB) gap-filling methods, which calculate the difference between clear-and cloudy-sky pixels based on the land SEB equation to estimate the actual LST under cloudy conditions [40]. For example, Jin [41] proposed an SEB-based theory to calculate missing LST information by interpolating adjacent pixels; however, this method does not take into account the terrain, which affects the accuracy of surface temperature reconstruction in mountainous areas. Yu et al. [42] introduced satellite radiation products into SEB theory, reducing the error caused by input variables, and thereby improving the generality and reliability of the SEB method. However, this method makes assumptions about environmental variables and requires specific meteorological and hydrological observations to calculate the LST difference between clear-sky and cloudy areas, which introduces errors [43].
The majority of the methods discussed above have a flaw in that the LST under cloudy conditions is represented by the LST theoretical value under clear-sky conditions. Although the filling methods based on multisource fusion can generate the actual LST, their verification accuracy may not be high. PMW can be regarded as one of the better methods of this kind, but it needs to overcome the rough spatial resolution and the differences among different LST data sources. Therefore, the technology is not quite mature at present. The weather research and forecasting (WRF) model can also obtain the LST under cloudy conditions, and WRF technology is more mature than PMW technology. WRF has five main modular physical schemes, namely, microphysics, shortwave radiation, longwave radiation, near-surface layer, and boundary layer schemes, which can simulate meteorological values under all weather conditions and surface skin temperature (TSK) with high temporal and spatial resolution [44], similar to MODIS. The LST retrieved from remote sensing observation data characterizes the integrated temperature of the surface skin with the surface thickness equal to the penetration depth (0.1-10 times the wavelength, detection of thermal radiation at a depth of several micrometers on the surface of the TIR observation spectrum) [45], [46]; the values of TSK and LST are very close, and can even be regarded as the same kind of data (TSK is treated as WRF LST in this article). Although the LST data simulated by WRF have spatial continuity and high spatiotemporal resolution, the simulation accuracy of this method is limited. In this article, considering that MODIS LST has high retrieval accuracy but is contaminated by clouds and has poor spatial continuity, numerical simulation with the WRF model, MODIS retrieval data, and the random forest (RF) model are integrated to give full play to the advantages of the WRF model and remote sensing data and generate a spatiotemporally complete LST. The advantage of the research method and source data in this article is the simulation of LST by WRF based on real historical observation data (final analysis data), which allows the method to simultaneously obtain LST with high accuracy in both clear-sky and cloudy conditions. First, to ensure an adequate number of valid pixels entering the RF model and improve the training accuracy of the model, we combine Terra/MODIS and AQUA/MODIS data, make full use of the spatiotemporal information of both, and create a "merged" dataset (Terra-Aqua/MODIS) by filling missing Aqua/MODIS data on the premise that Terra/MODIS data are available [47]. We then generate as many MODIS pixels as possible and use the RF model as a bridge to build a model of the relationship between MODIS LST and WRF LST. This procedure can improve not only the continuous availability of retrieval data but also the accuracy of the simulated data. We then obtain the actual LST in cloudy conditions to study SUHIs in Tianjin at night and during the day and analyze the deviations between cloudy-and clear-sky MODIS SUHIs. The purpose of this article is to combine the advantages of remote sensing and meteorological techniques, consider the spatial properties of meteorological simulations, and compensate for the disadvantage of remote sensing data, i.e., its limitation by cloud contamination, which has hindered its scientific application.
This article is structured as follows: "Study area and data collection" introduces the study area and data sources; "Methods" introduces the data processing methods and technical routes; and the section "Results" presents the WRF model simulation results, the RF optimization results, and the SUHI results. "Discussion" describes the strengths of this article, problems with the model simulation and optimization methods, and potential future improvements; and the section "Conclusion" presents a summary of the findings of this article.

A. Study Area
The study area in this article is Tianjin, a provincial administrative area of China that covers a total area of 11966.45 km 2 and has a population of approximately 13.86 million people (Fig. 1). Tianjin is located between 116°43 and 118°04 E longitude and between 38°34 and 40°15 N latitude, with flat topography, a temperate continental monsoon climate, high temperatures, rainy summers, and cold and dry winters. The annual precipitation in this area is approximately 550 mm, and the annual average temperature is 13.7°C. During the urbanization process, great temperature differences have arisen between urban and suburban areas of Tianjin, leading to an obvious heat island effect [48]. Studying the UHI effect in the area will help policy-makers formulate corresponding climate change measures and promote sustainable urban development.

B. Data on LSTs
MODIS is a 36-band sensor onboard the Terra [Earth Observing System (EOS) AM] and Aqua (EOS PM) satellites, which cover the Earth's surface four times a day, twice during the day (at 10:30 and 13:30 local time) and twice at night (at approximately 22:30 and 1:30 local time) [24]. Two kinds of daily satellite-derived MODIS LST products were used in this article, namely, the MOD11A1 data of the Terra/MODIS satellite collected in summer and winter from 2013 to 2018 and the MYD11A1 data of the Aqua/MODIS satellite. The MODIS dataset was derived from the National Aeronautics and Space Administration provided EOS data (EOSDIS, https://earthdata. nasa.gov/) at a spatial resolution of 1 km. The MODIS LSTs were generated from two TIR bands, bands 31 (10.78-11.28 μm) and 32 (11.77-12.27 μm), using a generalized split-window algorithm [49], [50]; the accuracy of these data on a uniform surface was above 1 K [25], [50]. The obtained MODIS data were spliced and reprojected to the world geodetic system-84 coordinate system and resampled to a resolution of 0.01°× 0.01°using the Arcpy tool.

C. Meteorological Data
The final operational global analysis data used in this article (http://rda.ucar.edu/datasets/ds083.2/) were reanalysis data produced and provided by the National Center for Environmental Prediction (NCEP); these data can be used to initialize the WRF model. The utilized dataset had the characteristics of multiple phases, strong continuity, and a high resolution [51]; the temporal resolution was 6 h, the spatial resolution was 1°, and the format was gridded binary (GRIB2). The weather station-recorded data used to verify the WRF simulation results were obtained from the National Oceanic and Atmospheric Administration (NOAA) and National Center for Environmental Information website (https://www.ncei.noaa.gov/), and the temporal resolution was 30 min.

D. Data on Land Cover and Regions
Land use/cover changes (LUCCs) are considered to be the main factor affecting the climate system and are among the key parameters used in coupled atmosphere-land surface models. Land-atmosphere interactions are closely related to the regional climatic environment background, topography, vegetation, and so on [51], [52]. The built-in MODIS land-use MCD12Q1 data (2001) in the WRF model with a spatial resolution of 500 m were used in this article. UHI research requires the division between urban areas and suburbs. To improve the accuracy of the division scheme, we used Landsat 8 land use data obtained from the website of the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences at a spatial resolution of 30 m.

A. Technical Process
The technical process of this article is shown in Fig. 2, and is divided into the simulation of LST with the WRF model, the amalgamation of valid MODIS data, and the retrieval of cloud-free MODIS-like LST by the WRFFM model and characterization of Tianjin's SUHIs based on the model results. First, the WRF model was used to simulate the generation of surface temperature data in Tianjin. For the acquired AQUA/MODIS data, if there were not enough pixels in the study area due to cloud contamination, the pixels were filled by Terra/MODIS data from the same day. The next step was to combine the WRF model, RF model, and spatial optimization algorithm to generate cloud-free MODIS-like LST data in the case of partial cloud contamination in MODIS; if the MODIS image was completely cloud-contaminated or Terra/MODIS also had insufficient pixels, the time-adjacent model was used. Finally, the SUHIs in Tianjin were quantitatively analyzed based on the cloud-free MODIS-like LST data generated by the above steps.

B. WRF Model and Its Parameter Configuration
In this article, we used version 4.2.1 of the regional WRF model; this model was developed as a joint effort among several agencies, including the National Center for Atmospheric Research, NCEP, and NOAA. A new generation of small-and medium-scale weather forecast models has been developed [53]. The main purpose of this model development is to study and predict weather conditions [54]. The utilized model has many modular schemes that can be used to describe land surface processes during the regional climate simulations [55], [56] and has been widely used in the meteorological field. The urban canopy model (UCM) can use the morphology of urban buildings to calculate the SEB [57], [58]; this model can accurately simulate LSTs at spatial resolutions of 1-10 km under cloudy and clear-sky conditions and surface energy fluxes [59]. In this article, we combined the WRF and UCM models using threelayer unidirectional nested domains with spatial resolutions of 9, 3, and 1 km (d01, d02, and d03, respectively, in Fig. 1); the pixel dimensions of these three domains were 97 × 97, 160 × 160, and 223 × 223, respectively. The WRF model has a physical scheme containing five modules: cloud, microphysics, radiation (longwave and shortwave), boundary layer, and surface modules [43]. Herein, we adopted the physical scheme described in Table I. According to the data recorded at the two automatic weather stations listed in Fig. 1(b), the WRF simulation results were verified, and the modeling performance was evaluated.

C. Amalgamation of Valid MODIS Data
The overpass times of the Aqua and Terra satellites are temporally close (within approximately 3 h). Assuming that most of the data are cloud contaminated, the information from these two data sources can be combined to create a "merged" dataset that complements the existing Aqua LST product and thereby reduce data loss. First, an adjustment method was applied to account for the different overpass times of the two platforms. Based on multiyear (2013-2020) MODIS LST data, the average Aqua LST and Terra LST values on Tianjin's image grid were determined according to seasons (winter: December, January and February; summer: June, July, and August). The mean Aqua-Terra LST differences were added to the daily Terra LSTs to remove the mean deviations related to the transit times, resulting in a merged Terra-Aqua/MODIS dataset. With independent offsets for each grid cell representing each season, it was possible to control the first-order statistics for factors, such as land cover, elevation, terrain slope and aspect, latitude, season, and snow cover during the diurnal surface temperature cycle [47]. The specific calculation formula can be expressed as follows: where T i is the summer or winter mean LST of the pixel grid of Terra/MOD11A1 data, A j is the summer or winter mean LST of the pixel grid of the Aqua/MYD11A1 data, n is the number of years (2013-2020), M difference is the 8-year summer or winter mean bias, and (Terra − Aqua) i is the pixel LST value of the combined Terra and Aqua datasets.

D. Retrieval and Optimization of LSTs Under Partial Cloud Conditions
The traditional RF model is a modern and highly flexible machine learning method that integrates multiple decision trees and can handle large amounts of missing data to provide high prediction accuracy [60]. In this article, two strategies were used for MODIS-like LST reconstruction under cloudy conditions based on the number of effective pixels in MODIS images (including Aqua/MODIS and the amalgamation Terra-Aqua/MODIS when cloud contamination was severe).
The reconstruction process when the MODIS LST is partially cloud contaminated (0% < cloud coverage < 99%) is shown in Fig. 3. After the simulations with the WRF model were obtained, the RF model was combined with the WRF model, where the input variables of the RF model are MODIS LST and WRF LST, the MODIS LSTs derived in July and December 2020 and the corresponding WRF LST data were taken as test data, and the image pixels input into the model were randomly divided into a training set and verification set at a ratio of 7:3. A model of the relationship between the MODIS and WRF data under clear-sky pixels was constructed, and RF LST images were obtained after the model was trained using the RF model. As the RF model is affected by the "regression mean" [61], salt and pepper noise arose; thus, it was necessary to use a spatial optimization algorithm to reduce this noise and smooth the images. Here, we combined the image prediction and optimization methods into a model named the WRFFM model and then used the root mean square error (RMSE) and Pearson correlation coefficient (Cor) to evaluate the accuracy of the WRFFM model outputs. Images with low Pearson correlations (Cor < 0.8) were regarded as exhibiting complete cloud contamination; in this way, high-quality cloud-free LSTs similar to those derived by MODIS could be obtained under cloudy conditions.

E. Retrieval and Optimization of LSTs Under Completely Cloud Conditions
In the case where the MODIS LST data were completely contaminated by clouds or the correlation of RF model training was low (Cor<0.8), the time-adjacent WRFFM model and the WRF image of the day were used to retrieve the LST. Fig. 4 presents a flowchart of this method. The idea is to incorporate the WRF simulated images of the current day into the WRFFM model of adjacent dates and introduce the Cor value of each adjacent date verified by the WRFFM model to obtain a MODIS-like LST with weight parameters for the corresponding date. The calculation formula can be expressed as follows: where LST c is the MODIS-like LST image, f i is the temporalneighbor training model, and n is the number of adjacent image dates used to estimate LST; in this article, n was set to 10. R is the correlation coefficient between the MODIS LST image and the simulated image, T is the number of days between image acquisition dates, LST cWRF is the WRF image of the current day, Cor is the Pearson correlation coefficient describing the relationship between the MODIS image and the LST cWRF of the current day, and ω is the weight of the image.

1) Calculated Indicators for SUHI:
Based on the spatial analysis capabilities of geographic information systems, in this study, we used high-spatial-resolution LUCC data from Landsat 8, reclassified the data, and extracted the impervious surface area (ISA) of Tianjin to identify urban areas. We then combined these data with the established buffer zone, which stratified the interior of the city according to ISA density; a 25% ISA threshold was used to define polygons based on thematic Landsat data, as studies have shown that the 25% ISA threshold can be used as a suitable boundary to determine the distance between a city and low-density residential land [62]. In this article, Heping District, Hedong District, Hexi District, Nankai District, Hebei District, and Hongqiao District were found to be the main urban areas in Tianjin with ISA densities >75%; 15-and 20-km buffer zones were built around these main urban areas [63]. Three concentric polygons were established, and the ISA density of the outermost polygon, which included Beichen District, Ninghe District, Jinnan District, Xiqing District, and Dongli District of Tianjin, was <25%; this delineation is similar to that utilized by Zhang et al. [64], who studied the continental USA. The indicators used to study the SUHI effect were consistent with the following simple UHI calculation formula: where SUHI is a surface urban heat island, Urban LST is the mean value of the urban LST, and Suburban LST is the mean value of the suburban LST. If an SUHI was negative, it was considered to be a surface urban cold island (SUCI). According to Chen and Wang [65], by using the equal spacing method to define the intensity grade of a UHI, UHIs were divided into seven grades, as shown in Table II .

2) Deviation Between MODIS and WRFFM SUHI:
Due to the influence of cloud contamination, the use of MODIS LSTs causes a lack of clear-sky pixels in the city and suburban areas in the resulting images. If the number of pixels in the area is too small, the estimated SUHI is unreliable due to the large potential deviations from the actual situation [66]. To quantify the bias (BMS) between the SUHI intensities calculated based on MODIS LST and WRFFM LST under different weather conditions and its influencing factors, in this article, we calculated the urban and suburban clear-sky coverage ratio (CSCR) and the biases of the ratio of urban clear-sky pixels to suburban clear-sky pixels under cloud contamination conditions and the ratio of urban pixels to suburban pixels under cloud-free conditions (BCF) CSCR = P clear−sky P all (7) BCF = P fre−urban P fre−suburbans − P con−urban P con−suburbans (8) where P clear−sky is the number of clear-sky pixels in the city and suburban areas, P all is the total number of pixels in the city and suburban areas without cloud contamination, P fre−urban and P fre−suburbans represent cases without cloud contamination, and P con−urban and P con−suburbans are the numbers of clear-sky pixels in cities and suburban areas, respectively, under cloud contamination conditions.  Fig. 1(b)], and the temperatures around urban areas were relatively low; that is, the WRF LSTs effectively captured the characteristics of UHIs, and they exhibited a similar spatial pattern to the MODIS LSTs of the corresponding date. To verify the difference between the WRF LSTs and the measured data, we compared the WRF LSTs obtained on July 6, July 29, December 20, and December 31, 2020, with the hourly temperatures recorded at the automatic weather stations to obtain the daily average temperature variation time series (Fig. 6). The daily maximum temperature occurred at 2:30 pm local time (UTC+8) and the lowest temperature occurred at approximately 6:00 pm local time. The Pearson correlation coefficient (Cor) values of the two data sources on July 6, July 29, December 20, and December 31, 2020 were above 0.8, with values of 0.818, 0.964, 0.809, and 0.816, respectively, and the corresponding RMSEs were 1.388, 0.896, 1.637, and 1.658, respectively, indicating that the WRF model optimally simulated the diurnal variations in air temperature, which is reasonable to some extent.

B. Generating Aqua-Terra/MODIS LSTs
Due to the influence of cloud contamination, the obtained Aqua/MODIS LST series had discontinuous data spaces due to missing pixels. However, sufficient effective pixels are required when fitting the WRF LSTs and MODIS LSTs using the RF model. In this article, we used Terra/MODIS clear-sky pixels to fill in the cloud-polluted pixels of the Aqua/MODIS data. Fig. 7 shows the average summertime 1) and wintertime 2) deviations in the Terra-Aqua data from 2013 to 2020. In Heping District, Hedong District, Hexi District, Nankai District, Hebei District, and Hongqiao District, the summertime surface temperature deviations were more obvious than those in winter. Taking the daytime Aqua/MODIS data obtained on July 5, 2020, as an example (Fig. 8), the pixels that filled the blank areas in the  Aqua data with Terra data were scattered, and the effective pixel coverage of the Aqua data was subsequently greatly increased, especially in the main urban area. Table III shows the percentages of clear-sky LST pixels before and after the merging of the Terra and Aqua data in July and December 2020 in Tianjin. Generally, there were more clear-sky pixels in the MODIS data in December than in July. By supplementing with Terra data obtained at the corresponding times, the effective pixels in the LST data obtained by the Aqua satellite in daytime and nighttime in July increased from 25.18% to 29.45% and from 34.44% to 35.56%, respectively, exhibiting effective clear-sky pixel coverage increases of 4.27% and 3.71%, respectively. In December, the daytime and nighttime values improved from 46.02% to 49.73% and from 46.46% to 49.59%, respectively, representing increases of 3.59% and 3.13%, respectively. These improvements played an important role in improving the fitting accuracy of the RF model.

C. WRFFM Model Results and Validation
In this section, we presented the simulation results of WRFFM model of this article. To further illustrate the results of our article, we presented simulations results of RF model, TRIMS LST proposed by Zhou et al. [67], and MKF-M LST images generated by Xu et al. [68]. TRIMS LST is generated by the satellite TIR remote sensing reanalysis data integration method based on the new surface temperature time decomposition model [69], [70], [71]. MKF-M LST product is generated based on MKF fusion TIR and PMW LST. Both TRIMS LST product and MKF-M LST product were obtained from the National Qinghai-Tibet Plateau Data Center http://data.tpdc.ac.cn/ [72], [73], and we compared them in spatial mode and quantitative accuracy (Fig. 9). First, under the condition that the MODIS data (collectively referred to as the Aqua and Terra-Aqua merged data in this article) were partially cloud-contaminated, the MODIS LST and WRF LST pixels were randomly divided into a training set and a verification set at a proportion of 7:3 for the RF model regression to fill in the blank MODIS LST area. Fig. 9(a), (f),  (k), and (p) shows the MODIS LST images representing July 6, July 22, December 12, and December 18, 2020. Due to cloud contamination, the pixels in these images are spatially discontinuous and exhibit blank areas. After the RF model was regressed, the blank areas were filled in, but a "salt and pepper noise" phenomenon was apparent; therefore, we used the spatial optimization algorithm to optimize the results. Fig. 9 shows that the WRFFM LSTs were relatively smooth following this filtering and exhibited very similar spatial distribution characteristics as the MODIS LSTs. Fig. 9(d), (i), (n), and (t) shows that MKF-M LST images are relatively "hazy." Fig. 9(e), (j), (o), and (u) shows that TRIMS LST images have relatively clear spatial distribution details which were similar to the MODIS LST.
In order to quantitatively compare the accuracy of WRFFM LST, MKF-M LST, and TRIMS LST images, a performance measure RMSE and Pearson correlation coefficient Cor were used to verify them (Fig. 10). The relationship between the MODIS original images and three product LSTs can be seen directly from the scatter density map of the numerical. The correlation coefficients corresponding to July 6, July 22, December 12, and December 18, 2020 were ∼0.875, ∼0.841, and ∼0.848. The correlation coefficients between TRIMS LST and MODIS LST are above 0.75, and the correlation coefficients of MKF-M LST are smaller than the former two. The RMSEs of WRFFM LSTs were ∼1.82, ∼3.06, ∼0.90, and ∼0.83 K, respectively. Compared with MKF-M LSTs and TRIMS LSTS, RMSE is smaller, indicating that the WRFFM LSTs were more consistent with the MODIS LSTs. The observed deviations can be explained by the upward sensible heat flux provided by the WRF model and the offset of solar shortwave radiation [74].
The MODIS images contaminated by clouds in July and December 2020 were filled, and the numbers of days in which valid data were input to the WRFFM model during the daytime in July and December were 22 and 28, respectively. Fig. 11 shows the changes in the Cor and RMSE values between the LSTs and MODIS LSTs from before to after the regression and spatial optimization algorithm were applied. The July RMSEs 1) were generally higher than those of December 2). The RMSEs of July were between 0.5 and 3.8 K, while the December RMSEs were between 0.4 and 1.8 K. The RMSEs of the optimized images were reduced compared to those obtained before optimization. The distribution ranges of the Cor values derived for the two months before and after optimization were similar. The Cor values derived before optimization were almost entirely below 0.8; after optimization, Cor in July and December ranged from 0.8 to 0.91 and 0.8 to 0.93, respectively, indicating that the model results were similar to the MODIS LSTs and could thus be used as MODIS-like cloud-free images to further research the SUHI effect.

D. SUHI Intensities Under Cloudy Conditions
Based on the cloud-free MODIS-like LST images optimized by the WRFFM model, the daytime and nighttime SUHIs derived in July (a) and December (b), 2020, were obtained, as  shown in Fig. 12. According to the grades of the SUHI effect on surface cities listed in Table II, on most days in July, SUHIs were more intense during the day than at night, and half of the days exhibited "weak SUHIs." From July 28 to 31, the daytime heat island effect was weaker than the nighttime effect. Even on July 19 and 30, "weak SUCI" phenomena appeared, and among nighttime observations, the number of days with "no SUHI" effect was highest, followed by the number of days with "weak SUHI" effects. On the whole, in December, nighttime UHIs were more intense than daytime SUHIs, and most identified SUHIs were nighttime "weak SUHIs"; almost no heat islands were observed in the daytime, with only one daytime "weak  SUHI" appearing, on December 8. In July, the average heat island temperature in daytime was 0.43 K higher than that at night, and in December, the average heat island temperature at night was 1.76 K higher than that in the daytime, indicating that the SUHI effect in Tianjin undergoes great seasonal differences.
To the best of our knowledge, current SUHI research that is based on LSTs obtained by MODIS is limited to clear-sky conditions, so the results cannot accurately reflect the actual SUHI intensities on certain days. To quantify the biases caused by performing MODIS-guided SUHI research only under clearsky conditions, we investigated the relationship between the CSCR and SUHI intensity bias (BMS); the results are shown in Fig. 13(a). As the number of filled clear-sky pixels increases, the overall SUHI intensity bias decreases, when the proportion of clear-sky pixels is below 30%, the deviation is above 3 K, and when the proportion of clear-sky pixels is above 80%, the deviation is below 0.6 K. When the pixel percentage is between 45% and 65%, the BMS value does not vary extensively.
According to the SUHI calculation index utilized in this article, if no cloud contamination is present in a city or suburbs, the ratio of the number of pixels representing cities and suburbs must be approximately 1:19; that is, each urban pixel needs to correspond to 19 suburban pixels. Under these conditions, the bias of the derived SUHI will not be excessively large. However, because the locations of cloud-contaminated MODIS LSTs are random, this ratio cannot be fixed, so researching the bias between the urban-to-suburban-pixel ratio under cloud-free and cloudy (BCF) conditions and the corresponding SUHI intensity bias (BMS) is necessary. Fig. 13(b) shows that as the BCF increases, according to formula (8), the proportion of clear-sky pixels in suburbs becomes smaller than that in urban areas; that is, when the ratio of urban pixels to suburban pixels exceeds 1:19, the BMS increases accordingly, with a maximum difference of 7.41 K. In summary, the BMS is greater under cloudy conditions than under cloud-free conditions, indicating that the MODIS SUHI has a large deviation from the actual SUHI under cloudy conditions, and that the number of clear-sky pixels in a city is an important factor influencing the bias relative to the actual SUHI.

A. Comparisons With Other Studies on LST Reconstruction Methods
When MODIS data were used to study the SUHI effect, spatial discontinuities caused by cloud contamination arose in the image data when the sensor obtained surface information, and blank areas appeared when pixels were missing. Such conditions can greatly affect research on the SUHI effect. Several studies [75], [76], [77] have addressed this issue by quantifying only nighttime heat islands, as daytime solar radiation increases the instability of the boundary layer and makes cloud formation more likely [13]. In addition, in some studies, only cloud-free LST images [78], [79], [80] have been selected for short-term case studies, or the mean has been used for long-term series studies; however, these approaches are not suitable for monitoring large areas. In this article, we combined remote sensing and meteorological methods to give full play to the advantages of both. By introducing WRF LSTs and MODIS LSTs into the RF model-fitting and subsequent optimization processes, clear-sky MODIS-like LST images with high spatial continuity and high accuracy were obtained. We then researched the nighttime and daytime SUHIs in Tianjin city. Like the WRFF model established by Fu et al. [43], the LSTs under cloudy conditions were reconstructed based on the WRF model and the RF model, but in the case of complete cloud contamination, LSTs were retrieved by using temporally adjacent images, which may not hold when the search for temporally adjacent images of the target date is extended over a longer period of time. In this article, the WRFFM model integrated the information of the two satellite MODIS images of the day to ensure the validity of the independent variables entering the RF model. Second, in the case of complete cloud contamination, we used the time-adjacent model for weighted calculation, while instead of temporally adjacent images, the results are more plausible. Compared with other filling methods, the model in this article has an advantage in filling pixels contaminated by clouds in that it considers land surface and atmospheric conditions when forming a coupled WRF/UCM simulation framework [81]. Furthermore, Terra and Aqua satellite-derived data were used to ensure that a sufficient number of effective pixels were available to input into the RF model for fitting, thus greatly improving the accuracy of the model retrieval results, and the WRFFM model can obtain LST under cloudless and cloudy conditions at the same time with high accuracy. In model verification (Fig. 10), the optimized RMSEs in July 2020 were found to be in the range of 0.5-4 K, while the December RMSEs ranged from 0.4 to 1.8 K. Compared with LST reconstruction algorithms of previous articles, such as PMW (yielding RMSEs between 2 and 5 K), and SEB model to infer LST under cloudy conditions (yielding RMSEs between 3 and 5 K), and the gap-filling method (yielding RMSEs of approximately 3.35 K) [82], [83], [84], [85], which mostly rely on spatiotemporal information [82], [86], the WRFFM model has strong overall performance, and the RMSEs are within an acceptable range.

B. Applications in SUHI Researches
The WRFFM model simulation results show that the daytime and nighttime SUHIs derived in Tianjin in July and December of 2020 exhibit obvious seasonal characteristics. Wang et al. [87] pointed out that urban-suburban surface thermal differences, which result in daytime UHIs, are driven by the seasonal solar heat radiation driver. In July, which is in summertime, the daytime solar thermal radiation of the Sun is intense, and heat is emitted from buildings and via transportation and the daily activities of residents throughout the city; in contrast, in suburban areas, the vegetation is denser and has a larger surface heat capacity, which can reduce the LSTs [88]. In winter, sunset occurs early, sunrise occurs late, the days are short, and the nights are long; these conditions are conducive to the formation of inversion layers [89]. In addition, in wintertime, there is no solar radiation at night, the suburban vegetation coverage is low, and the heat capacity of the land surface is relatively small [90], [91]; thus, the nighttime heat islands are more intense than the daytime heat islands in wintertime. A comparison of the WRFFM SUHIs with the cloud-contaminated MODIS SUHIs revealed that, the number of urban/suburban clear-sky pixels was not significantly negatively correlated with the SUHI deviation bias; this result was observed because cloud layers have a mitigating effect on SUHI intensity [92] and because the strength of this effect differs between day and night. Due to the randomness of the spatial distribution of clouds, the ratio of the number of urban pixels to the number of suburban pixels is greater than 1:19, so high numbers of urban and suburban clear-sky pixels can be observed. However, the calculated SUHI BMS deviations were large. Romanov [93] studied the influence of cloud cover on UHIs and found that cloud cover was more expansive in summer than in winter, which is consistent with our findings (Fig. 10). The height of the boundary layer is generally higher in urban regions than in suburban areas due to the higher surface temperatures of urban areas; thus, low-level coupling on the urban surface can affect the formation and location of clouds [94]. Cloudy urban environments affect SUHI research [95], and remote sensing can play an optimal role only under clear and cloudless conditions; thus, the results obtained with the model developed in this article can provide a reference for long-term SUHI monitoring.

C. Limitations and Future Perspectives
Because MODIS does not yield the correct LST under clouds, our model predicted the LST under clouds based on clear-sky pixels. However, as neither the validation set nor the pixels under clouds participated in the model training, the validation set can be regarded as the pixels under clouds. The accuracy attained with the two sets is the same, and the accuracy of the validation set is included in the overall accuracy, so the LST filled by this model can be regarded as the LST under clouds. This article has some limitations. The LSTs simulated by the WRF model developed in this article were more sensitive to the utilized physical scheme than the results obtained using traditional models [96], although a complete analysis of the WRF parameter changes was computationally expensive [97]. However, as cloud computing resources become available, the proposed WRFFM model has high application potential. In future article, data assimilation techniques can be introduced to use multisource remote sensing data to drive the WRFFM model and thereby estimate LST under cloudy conditions with greater accuracy. In addition, in this article, we used the results of the WRFFM model to quantify SUHIs in July and December 2020. In our future article, we will continue to investigate in detail the changes in derived long-term SUHI series and the driving factors.

VI. CONCLUSION
In this article, we took Tianjin as the research area, and by combining remote sensing and meteorological methods for application under cloudy conditions, merged the WRF model and RF method to construct a highly practical model (WRFFM) for generating high-quality MODIS-like LST data to research SUHIs. This article shows that the WRF simulation results are consistent with the hourly temperature data monitored by local automatic weather stations, and the extracted four-day verification results show that the correlation coefficients were above 0.8 and the RMSEs were between 0.8 and 1.7 K. To use Terra/MODIS-derived surface temperatures to fill gaps in Aqua/MODIS data caused by cloud contamination, we adjusted for the different overpass times of the two satellite platforms at different overpass times. The effective pixel coverage of the Terra-Aqua/MODIS data was compared with that of the Aqua/MODIS data, and the coverage of the Terra-Aqua/MODIS results during the daytime and nighttime were increased by 4.27% and 3.71%, respectively, in July and by 3.59% and 3.13%, respectively, in December. The WRFFM model results show that the WRFFM LSTs and MODIS LSTs exhibit very similar spatial patterns and can capture the UHI effect well; the correlation coefficients describing the relationship between these two datasets were above 0.8. The RMSE ranges were 0.5-3.8 K and 0.4-1.8 K in July and December, respectively; thus, the model accuracy was high. The SUHIs identified in Tianjin in July and December based on the WRFFM model results suggested that July was characterized as "high in the day and low in the night" SUHIs, and the average daytime heat island was 0.73 K warmer than the average nighttime heat island. December was characterized "low in the day and high in the night" SUHIs, and the average nighttime heat island was 1.81 K warmer than the average daytime heat island, revealing obvious seasonal differences. When the clear-sky pixel percentage was below 30%, the bias was above 3 K; when the clear-sky pixel percentage was above 80%, the bias was below 0.6 K. Thus, filled city pixels are an important factor affecting the measured SUHI intensity.