Reconstruction of All-Weather Daytime and Nighttime MODIS Aqua-Terra Land Surface Temperature Products Using an XGBoost Approach

Generating spatiotemporally continuous land surface temperature (LST) data is in great demand for hydrology, meteorology, ecology, environmental studies, etc. However, the thermal infrared (TIR)-based LST measurements are prone to cloud contamination with missing pixels. To repair the missing pixels, a new XGBoost-based linking approach for reconstructing daytime and nighttime Moderate Resolution Imaging Spectroradiometer (MODIS) LST measurements was introduced. The instantaneous solar radiation and two soil-related predictors from China Data Assimilation System (CLDAS) 0.0625°/1-h data were selected as the linking variables to depict the relationship with instantaneous MODIS LST data. Other land surface properties, including two vegetation indices, the water index, the surface albedo, and topographic parameters, were also used as the predictor variables. The XGBoost method was used to fit an LST linking model by the training datasets from clear-sky pixels and was then applied to the MODIS Aqua-Terra LSTs during summer time (June to August) in 2017 and 2018 across China. The recovered LST data was further rectified with the Savitzky–Golay (SG) filtering method. The results showed the distribution of the reconstructed LSTs present a reasonable pattern for different land-cover types and topography. The evaluation results using in situ longwave radiation measurements showed the RMSE varies from 3.91 K to 5.53 K for the cloud-free pixels and from 4.42 K to 4.97 K for the cloud-covered pixels. In addition, the reconstructed LST products correlated well with CLDAS LST data with similar LST spatial patterns. The variable importance analysis revealed that the two soil-related predictors and the elevation variable are key parameters due to their great contribution to the XGBoost model performance.


Introduction
Land surface temperature (LST) is a key parameter in environmental change, ecological processes and land-atmosphere energy exchange for different spatial scopes [1][2][3][4]. It enables monitoring of agricultural drought, urban heat, surface energy fluxes, and hydrological and meteorological processes [5][6][7][8][9][10][11], etc. In situ measurement is a credible source to acquire accurate LST data, but it is hard to obtain spatiotemporally continuous measurements for large-scale areas [12]. Satellite remote sensing has provided the unique approach to obtaining LST with satisfactory revisit cycle and global coverage [13,14]. Thermal infrared (TIR)-based methods, such as single-channel, split-window, and mono-window algorithms [13,[15][16][17][18], have been widely used to retrieve satellite-based LST data. Passive microwave (PMW)-based approaches have also been used for LST estimation with the advantage of better measurements through clouds [19,20]. Compared with LST retrieval from PMW sensors, the TIR-based LST data has attracted more attention due to its finer surface temperature (SST) and soil moisture (SM)) on LST or their interrelation has not been evaluated.
To address the abovementioned limitations, the purpose of this study is to develop a more practicable approach for reconstructing all-weather LST for cloud-covered pixels of the MODIS Aqua-Terra LST products. Considering that incident solar radiation is the key factor for the equilibrium between incoming and outgoing energy and determines the final LST values, Zhao et al. [26] used a modified random forest (RF) linking model proposed by Zhao et al. [58] to reconstruct daytime MODIS LST product. The RF linking model was conducted to link daytime LST with the solar radiation factor and other explanatory factors. However, the incident solar radiation made a small contribution in establishing the RF linking model of Zhao et al. [26], thus more appropriate variables need to be adopted to make the model more convincing. On the other hand, their developed RF linking models are only suitable for the reconstruction of daytime LST. In this study, in addition to the solar radiation factor, the instantaneous soil-related properties (SST and SM) were also considered as the key linking variables. The eXtreme Gradient Boosting (XGBoost) machine learning (ML) model, with its capacity of large-scale data modeling, was used to construct the relationship between LST with soil properties and other land surface properties (vegetation indices, water index, surface albedo, solar radiation and topographic parameters). Then, the built model was applied to the missing pixels to acquire the LST estimates and used for spatially continuous LST mapping over China. Finally, the reconstructed results were evaluated against LST measurements from China Meteorological Administration (CMA) and ground measurements in the Heihe River Basin, and additionally compared with China Land Data Assimilation System (CLDAS) LST data.
The major objectives of this study are to: (1) introduce an XGBoost-based algorithm to repair the missing information of daytime and nighttime MODIS Aqua-Terra LST data simultaneously; (2) assess whether the reconstructed LST results present the reasonable patterns for different land-cover types; (3) evaluate whether the accuracy of the reconstructed LST products meets the requirement compared with other studies, and (4) explore variable importance of each predictor and prove the rationality of the variable selection.

Study Area
China is selected as the study area ( Figure 1). With an area of approximately 9.6 million km 2 , China lies between latitudes 3 • N and 54 • N, and between longitudes 73 • E and 136 • E. The map of land-cover types in China is shown in Figure 1. The 300 m land-use cover product in 2018 generated from European Space Agency (ESA) Climate Change Initiative (CCI) project was used in this study (http://maps.elie.ucl.ac.be/CCI/viewer/index.php). The land-cover data from ESA was resampled to a 1 km resolution by mode resampling and reclassified into six types including water, bare soil, built-up areas, grassland, woodland, and cropland, for the purpose to simplify the evaluation of the impact of different land cover on LST reconstruction performance. China features high spatial heterogeneity in land-cover types, where the plains and basins account for~33% of the land area, while mountainous areas, hills and plateaus account for the other 67% of area. Due to the complicated terrain, the climate in China varies greatly in space and is mainly dominated by dry seasons and wet monsoon [59].

MODIS Data
In this study, the MODIS Aqua-Terra daily LST products (i.e., MOD11A1/MYD11A1) were selected to reconstruct cloud-covered LSTs. The auxiliary data from MODIS products used as the predictive variables were the normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI), the normalized difference water index (NDWI), and surface albedo.

MODIS Data
In this study, the MODIS Aqua-Terra daily LST products (i.e., MOD11A1/MYD11A were selected to reconstruct cloud-covered LSTs. The auxiliary data from MODIS pro ucts used as the predictive variables were the normalized difference vegetation ind (NDVI), the enhanced vegetation index (EVI), the normalized difference water ind (NDWI), and surface albedo.
The 16-day NDVI and EVI products were provided by the MODIS Aqua-Terra V etation Indices datasets (MOD13A2/MYD13A2). The NDWI were calculated from the s face reflectance product (MOD09A1) at an 8-day/500 m resolution and are given by where nir ρ and ρ swir denote band 2 and band 5 of the MOD09A1 product, resp tively. The surface albedo data were derived from the albedo product (MCD43A3) at an day/500 m resolution. The NDWI and surface albedo datasets were upscaled to 1 km r olution by pixel averaging to match MODIS LST data. All MODIS products were acquir from June to August in 2017 and 2018, and were obtained from the NASA Land Proces Distributed Active Archive Center (LPDAAC) (https://la sweb.modaps.eosdis.nasa.gov/). The 16-day NDVI and EVI products were provided by the MODIS Aqua-Terra Vegetation Indices datasets (MOD13A2/MYD13A2). The NDWI were calculated from the surface reflectance product (MOD09A1) at an 8-day/500 m resolution and are given by where ρ nir and ρ swir denote band 2 and band 5 of the MOD09A1 product, respectively. The surface albedo data were derived from the albedo product (MCD43A3) at an 8-day/500 m resolution. The NDWI and surface albedo datasets were upscaled to 1 km resolution by pixel averaging to match MODIS LST data. All MODIS products were acquired from June to August in 2017 and 2018, and were obtained from the NASA Land Processes Distributed Active Archive Center (LPDAAC) (https://ladsweb.modaps.eosdis. nasa.gov/).

Reanalysis Data
The CMA (China Meteorological Administration) China Land Data Assimilation System (CLDAS) started to release the official products (Version 2.0, 0.0625 • × 0.0625 • ) from 2008, and these products cover East Asia (0-65 • N, 60-160 • E). Meteorological fields used in the CLDAS land surface modeling (LSM) include hourly air temperature, air pressure, humidity, precipitation, solar radiation, surface temperature, shortwave radiation, soil surface temperature, and soil moisture, etc. The CLDAS reanalysis data are derived from LSMs using data fusion and assimilation methods [60,61].
In this study, the hourly CLDAS shortwave radiation, soil surface temperature (0-10 cm), and soil moisture (0-10 cm) were selected as the predictors of LST. As the CLDAS LST product proves to be well consistent with in situ LST [44], the CLDAS LST product was used to evaluate the reconstructed LST. In addition to the shortwave radiation parameter obtained during the daytime, CLDAS soil parameters at 11:00 a.m., 22:00 p.m., 2:00 p.m., and 2:00 a.m. (Beijing time) were selected to match the overpass time of MODIS Aqua-Terra products. These hourly reanalysis data in NetCDF format are publicly accessible from the China Meteorological Data Service Center (CMDC) at http://data.cma.cn.

Topographic Parameters
Topographic parameters affect variations in LSTs because of the varying atmospheric forcing environment and with the effect on the dispersion of incident solar radiation [26,62]. The DEM data was from the Shuttle Radar Topography Mission (SRTM) at a 90 m resolution and obtained from the Geospatial Data Cloud (http://www.gscloud.cn/). Then the slope data was calculated from the SRTM data, and the two topographic datasets were upscaled to 1 km resolution by pixel averaging.

Ground-Measured Data
Two ground-based datasets were used for the evaluation. One in situ measured LST dataset was obtained from "China's surface climate data daily value data set (V3.0)" (http://data.cma.cn) provided by CMA. This dataset, which has been strictly qualitycontrolled, contains data from 699 benchmarks and basic automatic weather stations in China and provides daily maximum, minimum, and average land surface temperature (0 cm) measured by platinum resistance temperature sensor. The maximum permissible error is ±0.2 • C (≤50 • C) or ±0.5 • C (>50 • C). The advantage of this dataset is that the stations cover the entire study area and can be used to evaluate the overall spatiotemporal characteristics of the reconstructed MODIS LST data. In order to make the station distribution more uniform, 200 automatic weather stations were randomly selected for verification by using the "Create Fishnet" module of ArcMap with a 2 • × 2 • cell ( Figure 1).
Another dataset is the ground-based measurements from the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) [63,64], a comprehensive eco-hydrological experiment conducted in the Heihe River Basin. This dataset measured by Net Radiometer CNR4 with a maximum nonlinear error of 1% has been used to assess the accuracy of the all-weather MODIS LST [49,65]. This dataset also has been quality-controlled. Five HiWATER sites (i.e., DSL, AR, HZZ, HH, and SDQ) provided by National Tibetan Plateau Data Center were used for verification (http://data.tpdc.ac.cn) in this study. Some details of these sites were displayed in Table 1. The ground-measured LSTs at the five sites were calculated from the upwelling and downwelling broadband hemispherical radiances by Stefan-Boltzmann's law: where T s is the estimated LST, F ↓ is the downwelling longwave radiation, F ↑ is the surface upwelling longwave radiation, σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), and ε b is the surface broadband emissivity. The ε b was estimated from MODIS narrowband emissivity products using an empirical linear relationship [66]. Additionally, the "3σ-Hampel identifier" was adopted to reduce the influences from outliers [65,67].  Table 2 gives an overview of the datasets used in this study, including the variable types, dataset name/source, and the spatiotemporal resolution. Ground-based surface temperature CMA Daily/point [69] In situ longwave radiation measurements HiWATER 10 min/point [70][71][72][73][74] 3. Methodology

Theoretical Context
LST is comprehensively affected by multiple factors, such as land-cover types, topographic parameters, soil parameters, atmospheric forcing conditions, etc. [26,[55][56][57]. Thus, many studies have considered multiple related factors to express the relationship between LST and these auxiliary factors, and these built relationships have been further used for LST reconstruction [26,50,58], LST downscaling [75], etc. The key of building these models is to find the variables that change synchronously and are closely related to the LST parameter, such as the linking variable (i.e., incident solar radiation) in Zhao et al. [58] and Zhao et al. [26], and the Annual Cycle Parameters (ACPs) in Sismanidis et al. [76]. In this study, we evaluated the utility of several linking variables as proxies for instantaneous LST over China, and reconstructed daytime and nighttime LSTs from MODIS observations.

Extreme Gradient Boosting (XGBoost) Model
Several studies have successfully used an ensemble learning method (i.e., RF method) to reconstruct cloud-covered LSTs [26,51,52]. The XGBoost algorithm belongs to the Gradient Boosting Decision Tree (GBDT) algorithm, which is an iterative decision tree algorithm consisting of multiple decision trees [77], and the final decision is made by iterating multiple trees together. XGBoost represents an efficient GBDT algorithm enabling gradient boosting "on steroids" [78]. It combines advanced optimization techniques to produce excellent results while using less computing resources than other methods [79]. This study involved LST modeling over a large area (i.e., China), and it was feasible to use XGBoost algorithm modeling with vast amounts of LST pixels.

LST Reconstruction Based on XGBoost Linking Model
The RF linking model proposed by Zhao et al. [58] and Zhao et al. [26] is based on the assumption that incident solar radiation dominates the morning warming process and this factor can accurately represent the close relationship between the daily LST and surface solar radiation [26,58,80,81]. In this study, for daytime MODIS Aqua-Terra LST's reconstruction, the incident solar radiation was also taken into account and was estimated by cumulating the CLDAS shortwave radiation data on surface warming process from sunrise time to satellite sunrise time. However, the cumulative incident solar radiation (CSR) is around zero during nighttime, thus the CSR cannot be used as an input to build a nighttime LST linking model, but the instantaneous CLDAS SST and SM data at an hourly resolution can be used as good proxies for daytime and nighttime LSTs due to their close relationships. Combined with other predictors based on Zhao et al. [26], including NDVI, EVI, NDWI, surface albedo (ALB), solar radiation, DEM, and slope (SLP), the proposed LST linking model based on XGBoost is expressed as follows: T act_day = T est_day + e day T est_day = F(NDV I, EV I, NDW I, ALB, CSR, SST, SM, DEM, SLP) T act_night = T est_night + e night T est_night = F(NDV I, EV I, NDW I, ALB, SST, SM, DEM, SLP) where Equation (3) is the daytime LST linking model, and Equation (4) is the nighttime LST linking model with the CSR factor removed. T act_day and T act_night denote the actual daytime and nighttime LST respectively, T est_day and T est_night denote the estimated daytime and nighttime LST respectively, e denotes the estimation error of the linking model, and F is the function of an established XGBoost model. The process of the LST reconstruction is described as follows: (1) Data Preprocessing: Before building the linking model, the quality control (QC) data of the MODIS Aqua-Terra LST products was used to remove pixel errors larger than 2 K to reduce the uncertainty derived from original LST data [26]. The NDVI or EVI data were compiled to generate two new time-series datasets with an 8-day resolution by combing the 16-day MOD13A3/MYD13A3 products. Then, all the MODIS-based datasets with an 8-day resolution (i.e., NDVI, EVI, NDWI, and ALB) were linearly interpolated to generate daily data. As for the hourly CLDAS-based predictors (i.e., SST and SM) with a 0.0625 • resolution, they were downscaled to 1 km resolution using a bilinear interpolation method.
(2) LST Linking Model Building: The XGBoost model was adopted to fit the LST linking model and to construct the relationship between clear-sky LSTs and the predictors. It should be noted that the CSR variable was eliminated for the nighttime LST linking model. We used the XGBoost package in scikit-learn, and most parameters were set to the default except that the parameter "max_depth" was tuned from 10 to 20 and the parameter "n_estimators" was tuned from 100 to 250 following the changes of "max_depth" to avoid overfitting according to the root-mean-square error (RMSE) of the training dataset. The mean squared error (MSE) was selected as the loss function because it is commonly used in regression issue. The selected optimal parameters were determined by the trade-off of model accuracy and time consumption for the XGBoost training model.
(3) LST Reconstruction: Among the input datasets for fitting the linking model, the predictors of Equations (3) and (4), except for CSR, SST, and SM, are mostly invariant within a day. The three instantaneous predictors (i.e., CSR, SST, and SM) are used as critical proxy variables to link the clear-sky LST pixels and cloud-covered LST pixels via the established linking model. To produce the reconstructed LST, the built model was applied to the region with missing pixels.
(4) Seamless LST Postprocessing: MODIS data are influenced by undetected clouds and poor atmospheric conditions, which caused discontinuous data existing in MODIS data [82,83]. The change curve of original LST often demonstrates sudden plunges due to cloud and snow cover, which caused sudden drop points that were not in a good agreement with the overall trend [33]. The SG filtering method has been used to smooth and reconstruct the LST time-series data to reduce the noise data. Two key parameters must be determined. The first parameter is m that denotes the half-width window, and the other parameter is d that denotes the degree of the smoothing polynomial. In this study, the two parameters were set to 2 and 3 respectively, and the SG filtering was only applied to the outlier data with values lower than the overall trend [33].
The flowchart of the whole process for daytime LST reconstruction is shown in Figure 2, and the process for nighttime LST reconstruction is similar with the daytime LST reconstruction. and snow cover, which caused sudden drop points that were not in a good agreement with the overall trend [33]. The SG filtering method has been used to smooth and reconstruct the LST time-series data to reduce the noise data. Two key parameters must be determined. The first parameter is m that denotes the half-width window, and the other parameter is d that denotes the degree of the smoothing polynomial. In this study, the two parameters were set to 2 and 3 respectively, and the SG filtering was only applied to the outlier data with values lower than the overall trend [33].
The flowchart of the whole process for daytime LST reconstruction is shown in Figure 2, and the process for nighttime LST reconstruction is similar with the daytime LST reconstruction.

Validation
In this study, in situ measured LST data were used to evaluate the reconstructed LST results. Besides, the CLADS LST data were also adopted to evaluate the spatial patterns

Validation
In this study, in situ measured LST data were used to evaluate the reconstructed LST results. Besides, the CLADS LST data were also adopted to evaluate the spatial patterns and LST magnitudes of the reconstructed results. Three commonly used criteria were selected for the validation, that is, the coefficient of determination (R 2 ), the root mean square error (RMSE), and the Bias. The three statistical metrics are calculated as follows: where Y i denotes the reconstructed LST value and O i denotes the corresponding in situ measured LST value or CLDAS LST value; Y and O are mean values of Y i and O i , respectively.

Building of the LST Linking Model
We reconstructed the daytime and nighttime MOD11A1/MYD11A1 data across China during summer time (June to August) in 2017 and 2018. Separate datasets including the LST data for each satellite overpass time and corresponding auxiliary data were independently imported into an XGBoost-based LST linking model. To verify the accuracy and stability of the model, a hand-out cross validation method was employed by randomly splitting the training datasets into 90% of the total datasets, and the remaining 10% of the datasets as the training dataset and test dataset, respectively. The parameter of "max_depth" was set to 18 and "n_estimators" was set to 200 as the best parameters according to the minimum MSE value for the test dataset and the computation time. Table 3

Demonstration of Reconstructed Daytime and Nighttime LSTs
To facilitate the ease of reading, the original daytime LST and nighttime LST of MOD11A1 data, their corresponding reconstructed daytime LST and nighttime LST as well as reconstructed daytime LST and nighttime LST after SG filtering are labeled as MODLSTD_Raw, MODLSTN_Raw, MODLSTD_Rec, MODLSTN_Rec, MODLSTD_Rec_SG, and MODLSTN_Rec_SG, respectively hereafter, and so were the MYD11A1 data.
The original cloud-covered LST data of MODLSTD_Raw and MODLSTN_Raw on DOY 167 (June 16) of 2018 were displayed in Figure 3a,b, and their corresponding gap-filled LST data and reconstructed LST data with SG filtering were displayed in Figure

Reconstruction Effects for Different Land-Cover Types
In order to show the reconstructed results in more detail, we selected two representative sub-regions, shown in Figure 5, as examples based on the reclassified land-cover types. Figure 6 shows the original MOD11A1 LST data and the corresponding reconstructed data on DOY 167 of 2018 for sub-region 1 scattered with water bodies. The daytime MODIS LST for pixels covered by water is usually lower than the surrounding pixels, and by contrast, the LSTs of nighttime water pixels are higher than the LSTs of their neighboring pixels. The reconstructed results (Figure 6c,d) can clearly reflect the terrain and land-cover information of the area with no obvious overestimation or underestimation. When the SG filtering was performed (Figure 6e,f), some underestimation and blurred information was rectified, though the improvement is not very evident for this sub-region.
The original cloud-covered LST data of MODLSTD_Raw and MODLSTN_Raw on DOY 167 (June 16) of 2018 were displayed in Figure 3a,b, and their corresponding gapfilled LST data and reconstructed LST data with SG filtering were displayed in Figure 3c,d and Figure 3e,f, respectively. Likewise, the original cloud-covered MYD11A1 LST data and their corresponding reconstructed results were shown in Figure 4. Figures 3 and 4 demonstrated that the XGBoost linking model can effectively fill a lot of missing LST pixels in the original LST images. The overall spatial patterns of the reconstructed LST images were consistent with the original LST data. More examples of reconstructed MODIS LST data on other days were provided in Figures A1-A3 in the Appendix A.

Reconstruction Effects for Different Land-Cover Types
In order to show the reconstructed results in more detail, we selected two representative sub-regions, shown in Figure 5, as examples based on the reclassified land-cover types. Figure 6 shows the original MOD11A1 LST data and the corresponding reconstructed data on DOY 167 of 2018 for sub-region 1 scattered with water bodies. The daytime MODIS LST for pixels covered by water is usually lower than the surrounding pixels, and by contrast, the LSTs of nighttime water pixels are higher than the LSTs of their neighboring pixels. The reconstructed results (Figure 6c,d) can clearly reflect the terrain and land-cover information of the area with no obvious overestimation or underestimation.   Figure 7a,b show that the values of LST pixels covered with bare or sparsely vegetated land are higher than that of other land-cover types for both daytime and nighttime. Because of the large blank area in Figure 7a,b covered by bare soil or sparse vegetation, the surrounding relatively low LST pixels were used as input to the XGBoost model and then to estimate those missing pixels. Therefore, the corresponding reconstructed LSTs of Figure 7c,d were not so accurate with the underestimation effects. With the process of SG filtering via using the temporally adjacent LST information, the underestimated pixels were well rectified (refer to Figure 7e,f) by combining with the information of adjacent time-series LST data. Additionally, the step of SG filtering can restore the LST information in areas covered with other land-cover types (e.g., woodland) (Figure 7e vs. Figure 7c). When the SG filtering was performed (Figure 6e,f), some underestimation and blurred information was rectified, though the improvement is not very evident for this sub-region.

Validation with CMA Ground-Measured LST Data
The CMA ground-measured LST (0 cm) daily maximum, minimum, and average LST data was used to assess the reconstructed MODIS LST products. In this study, the daytime MYD11A1 LST and nighttime MYD11A1 LST were taken as proxies for the daily maximum and minimum LST, respectively. Additionally, we average the reconstructed daytime and nighttime MOD11A1 and MYD11A1 LST for each day to match the ground-measured average LST, and the average LST for both before and after SG filter processing were labeled as MODLSTAvg_Rec and MODLSTAvg_Rec_SG, respectively.
The validation results of the reconstructed MODIS LSTs (from June to August in 2017 and 2018) for both before and after SG filter processing were shown in Figures 8 and 9. The distribution of density scatter points in Figures 8a and 9a indicates that most of the LST values of the MYDLSTD_Rec data were far lower than the values of in situ measured maximum LST data (Bias = −17.00 K to −16.30 K). This suggests that the daytime MYD11A1 LST at local 2:30 p.m. does not give a good approximation of the daily maximum LST with relatively low R 2 and high RMSE (R 2 = 0.40 to 0.42, RMSE = 18.62 K to 19.48 K). Figures 8b and 9b indicate that the nighttime MYD11A1 LST at local 2:30 a.m. correlated relatively well with the in situ measured daily minimum LST (R 2 = 0.62 to 0.65, RMSE = 5.74 K to 6.09 K, Bias = −3.32 K to −3.18 K). The result of Figures 8c and 9c indicate that the reconstructed daily average MODIS LST also has a relatively good correlation with the ground-measured mean LST (R 2 = 0.62 to 0.63, RMSE = 7.07 K to 7.45 K, Bias = −5.37 K to −3.32 K). By comparison with the results of reconstructed LSTs with SG filter processing, we found little differences between  Figure 9b) because the SG filter processing mainly improved the abnormally low LST values in the reconstructed LST products. From the results of the daily average LST, the R 2 value was slightly reduced from 0.63 to 0.61 (Figure 9f vs. Figure 9c) or remained unchanged by comparing Figures 8f and 8c. Thus, the step of SG filtering has little influence on the accuracy of the preliminary reconstructed LST product on the whole. Figure 5. Two representative sub-regions for comparison in more detail.  To further evaluate the accuracy for different land-cover types, Table 4 lists the validation of the reconstructed products for 6 land-cover types. The correlations for the three land-cover types of built-up areas, grassland, and cropland were higher than the other three land-cover types (i.e., water, bare soil, and woodland). In particular, for the reconstructed products of MYDLSTN_Rec_SG and MODLSTAvg_Rec_SG with the three land-cover types of built-up areas, grassland, and cropland, the R 2 ranges from 0.59 to 0.80 during summer time in 2017 and 2018, the RMSE ranges from 2.41 K to 8.18 K, and the Bias ranges from −6.24 K to −0.91 K. The R 2 for the other three land-cover types were relatively lower, while the RMSE ranges from 2.30 K to 9.21 K and the Bias ranges from −8.58 K to 0.56 K. As for the MYDLSTD_Rec_SG product, the poor validation results for all land-cover types were attributed to the large bias between the reconstructed MYDLSTD_Rec_SG and the in situ measured maximum LST data. the corresponding reconstructed LSTs of Figure 7c,d were not so accurate with the und estimation effects. With the process of SG filtering via using the temporally adjacent L information, the underestimated pixels were well rectified (refer to Figure 7e,f) by co bining with the information of adjacent time-series LST data. Additionally, the step of filtering can restore the LST information in areas covered with other land-cover types (e woodland) (Figure 7e vs. Figure 7c).

Validation with CMA Ground-Measured LST Data
The CMA ground-measured LST (0 cm) daily maximum, minimum, and average L data was used to assess the reconstructed MODIS LST products. In this study, the daytim MYD11A1 LST and nighttime MYD11A1 LST were taken as proxies for the daily ma mum and minimum LST, respectively. Additionally, we average the reconstructed da time and nighttime MOD11A1 and MYD11A1 LST for each day to match the groun measured average LST, and the average LST for both before and after SG filter processi were labeled as MODLSTAvg_Rec and MODLSTAvg_Rec_SG, respectively.

Validation Using HiWATER Data
The accuracy of the all-weather LSTs were also evaluated against the HiWATER LST measurements at five sites and the results are shown in Figures 10 and 11. The RMSE of the original cloud-free LST varies from 3.95 K to 5.53 K in 2017 and from 3.91 K to 4.75 K in 2018 at all sites. Duan et al. [49] pointed out that the spatial scale inconsistency between the pixel-wise all-weather LSTs and the point-wise in situ LSTs caused the relatively large RMSE. Another reason for this relatively low accuracy stemmed from the influence of cloud contamination. The accuracy evaluation used nearly all QC data rather than the best QC data (QC = 0) [49]. Thus, the accuracy of the all-weather LST is obviously lower than the nominal accuracy (1 K) of the MODIS LST product. From Figures 10 and 11 vs. Figure 8a-c or Figure 9d-f vs. Figure 9a-c), and the R value was improved from 0.62 to 0.65 (Figure 8e vs. Figure 8b) or from 0.65 to 0.68 (Figure 9e vs. Figure 9b) because the SG filter processing mainly improved the abnormally low LST values in the reconstructed LST products. From the results of the daily average LST, the R 2 value was slightly reduced from 0.63 to 0.61 (Figure 9f vs. Figure 9c) or remained unchanged by comparing Figure 8f and Figure 8c. Thus, the step of SG filtering has little influence on the accuracy of the preliminary reconstructed LST product on the whole.   To further evaluate the accuracy for different land-cover types, Table 4 lists the validation of the reconstructed products for 6 land-cover types. The correlations for the three land-cover types of built-up areas, grassland, and cropland were higher than the other three land-cover types (i.e., water, bare soil, and woodland). In particular, for the reconstructed products of MYDLSTN_Rec_SG and MODLSTAvg_Rec_SG with the three landcover types of built-up areas, grassland, and cropland, the R 2 ranges from 0.59 to 0.80

Validation Using HiWATER Data
The accuracy of the all-weather LSTs were also evaluated against the HiWATER LST measurements at five sites and the results are shown in Figures 10 and 11. The RMSE of the original cloud-free LST varies from 3.95 K to 5.53 K in 2017 and from 3.91 K to 4.75 K in 2018 at all sites. Duan et al. [49] pointed out that the spatial scale inconsistency between the pixel-wise all-weather LSTs and the point-wise in situ LSTs caused the relatively large RMSE. Another reason for this relatively low accuracy stemmed from the influence of cloud contamination. The accuracy evaluation used nearly all QC data rather than the best QC data (QC = 0) [49]. Thus, the accuracy of the all-weather LST is obviously lower than the nominal accuracy (1 K) of the MODIS LST product. From Figures 10 and 11

Comparison with CLDAS LST Data
Taking the reconstructed results on DOY 167 of 2018 as an example, we compared them with the corresponding CLDAS LST products with an hourly/0.0625° resolution. The left column of Figure 12 shows the four CLDAS LST products and the right column shows the reconstructed results on DOY 167 of 2018. It is shown that the CLDAS LST products contain some missing pixels and noise data to the west of China. The reconstructed LST products depicted more detailed spatial information than the CLDAS LST, especially in areas with complex terrain (e.g., Tibetan Plateau). Overall, the reconstructed LST data had similar spatial variation with the CLDAS LST data. Differences in the LST magnitude are seen in some areas. For example, the reconstructed daytime MODIS LST data ( Figure  12b,f) in northwest China were higher than those of the daytime CLDAS LST data ( Figure  12a,e). In contrast, the LST estimates of Figure 12b,f in southeast China were lower than those of Figure 12a,e. For the nighttime LST data, the differences between the reconstructed MODIS LST data (Figure 12d,h) and CLDAS LST data (Figure 12c,g) were much smaller compared to the daytime LST data.

Comparison with CLDAS LST Data
Taking the reconstructed results on DOY 167 of 2018 as an example, we compared them with the corresponding CLDAS LST products with an hourly/0.0625 • resolution. The left column of Figure 12 shows the four CLDAS LST products and the right column shows the reconstructed results on DOY 167 of 2018. It is shown that the CLDAS LST products contain some missing pixels and noise data to the west of China. The reconstructed LST products depicted more detailed spatial information than the CLDAS LST, especially in areas with complex terrain (e.g., Tibetan Plateau). Overall, the reconstructed LST data had similar spatial variation with the CLDAS LST data. Differences in the LST magnitude are seen in some areas. For example, the reconstructed daytime MODIS LST data (Figure 12b,f) in northwest China were higher than those of the daytime CLDAS LST data (Figure 12a,e). In contrast, the LST estimates of Figure 12b,f in southeast China were lower than those of Figure 12a,e. For the nighttime LST data, the differences between the reconstructed MODIS LST data (Figure 12d,h) and CLDAS LST data (Figure 12c,g) were much smaller compared to the daytime LST data.
The reconstructed LST products were aggregated to 0.625 • and matched-up to the CLDAS LST pixels by average aggregating to further understand the correlation between the reconstructed LST products with the CLADS LST products. Figures 13 and 14 show the density scatter plots for the reconstructed LST data and the CLADS LST data from June to August in 2017 and 2018. From Figures 13 and 14, some high outliers are seen in both daytime and nighttime LST estimates from CLDAS. The correlation between these two LST datasets varies (R 2 = 0.33 to 0.80), while the correlation of nighttime LST data (R 2 = 0.77 to 0.80) is significantly higher than that of the daytime data (R 2 = 0.33 to 0.51). The RMSEs are relatively large for daytime estimates (RMSE = 8.05 K to 10.10 K), while the absolute Bias is smaller than 5 K for both daytime and nighttime estimates. Remote Sens. 2021, 13, x FOR PEER REVIEW 18 of 28  CLDAS LST pixels by average aggregating to further understand the correlation between the reconstructed LST products with the CLADS LST products. Figures 13 and 14 show the density scatter plots for the reconstructed LST data and the CLADS LST data from June to August in 2017 and 2018. From Figures 13 and 14, some high outliers are seen in both daytime and nighttime LST estimates from CLDAS. The correlation between these two LST datasets varies (R 2 = 0.33 to 0.80), while the correlation of nighttime LST data (R 2 = 0.77 to 0.80) is significantly higher than that of the daytime data (R 2 = 0.33 to 0.51). The RMSEs are relatively large for daytime estimates (RMSE = 8.05 K to 10.10 K), while the absolute Bias is smaller than 5 K for both daytime and nighttime estimates.

Variable Importance Analysis
The XGBoost algorithm enables assessment of variable importance to rank their contributions to the established model. The average importance of each predictor for the built XGBoost models during summer time in 2017 and 2018 in this study is shown in Figure  15. Figure 15 demonstrates that the variable importance differs with satellite overpass

Variable Importance Analysis
The XGBoost algorithm enables assessment of variable importance to rank their contributions to the established model. The average importance of each predictor for the built XGBoost models during summer time in 2017 and 2018 in this study is shown in Figure 15. Figure 15 demonstrates that the variable importance differs with satellite overpass time. For daytime LST estimation, predictors such as SST, SM, DEM, and NDVI provided significant contributions to the XGBoost models (Figure 15a,c). For nighttime LST estimation, the three variables of DEM, SST, and SM exhibit more importance than the other variables (Figure 15b,d). The variables of SST, SM, and DEM are key inputs in the XGBoost models for both daytime and nighttime LST estimates due to their large contributions. Some studies have explored the relationship between SST and LST to retrieve SST [84,85], and Zhang et al. [51] even used SST data to validate the fused all-weather LST. Sun et al. [56] estimated the effects for SM on the retrieval accuracy of LST, and Jiang et al. [86] estimated hourly SM by using downscaled LST over different urban surfaces. For DEM data, Zhao et al. [26] pointed out there is a notable decreasing trend for surface elevation with LST due to the altitudinal temperature gradients. In the study area, surface elevation is also significantly correlated negatively with LST for both daytime and nighttime.

Comparison with Other Studies
The reconstructed MYD11A1 LST data did not agree well with the CMA LST measurements for daytime data. One reason is that there was a difference between the satellite overpass time of MODIS LST data across China and the measurement time for CMA LST data. The gap time inevitably led to the large difference. On the other hand, Jiang et al. [87] compared the 0 cm LST measured by mercury thermometer with the LST retrieved by upward and downward longwave radiation in Tibetan Plateau. Based on the results,

Comparison with Other Studies
The reconstructed MYD11A1 LST data did not agree well with the CMA LST measurements for daytime data. One reason is that there was a difference between the satellite overpass time of MODIS LST data across China and the measurement time for CMA LST data. The gap time inevitably led to the large difference. On the other hand, Jiang et al. [87] compared the 0 cm LST measured by mercury thermometer with the LST retrieved by upward and downward longwave radiation in Tibetan Plateau. Based on the results, the highest difference in daytime is over 16 K so the results in this study were consistent with those of Jiang et al. [87]. The reason is that the measurement error of mercury thermometer is discrete and non-proportional, thus the representation of the LST measured by mercury thermometer is poor, even making the average or integral process [87]. It is essential to popularize the surface temperature measurement method based on the surface radiation to meet the demand of higher accuracy and reliability [87]. However, although this CMA dataset also has inherent systematic bias, the reconstructed nighttime MYD11A1 LST data correlated well with the CMA LST measurements, which demonstrated the good spatiotemporal continuity and relative accuracy of the reconstructed nighttime MODIS LST products.
Validation against radiation-based LST revealed stable performance at all five HiWA-TER sites, and the accuracy assessment was in line with previous studies. For example, the NP-based approach used by [42] was evaluated with two sites in Africa with the RMSE ranging from 5.11 K to 5.55 K. Duan et al. [49] merged TIR and PMW measurements in China with the RMSE varying from 3.5 K to 4.4 K. The study of Zeng et al. [12] showed the Mean Absolute Error (MAE) ranged from 3-6 K at six sites in America. Pham et al. [88] obtained the gap-filling MYD11C1 LST over Australia with the RMSE ranging from 2 K to 6.4 K at four stations in America. Here, the reconstructed MODIS data using the XG-Boost method reached an RMSE between 3.91 K to 4.97 K, demonstrating an overall good performance with respect to the absolute accuracy.
Zhao et al. [58] proposed an RF linking model to reconstruct MODIS LSTs over a 120 km × 200 km study area (24,000 km 2 ), and Zhao et al. [26] developed this initial RF linking model to recover daytime MODIS LSTs over a region of Southwest Europe (approximately 600,000 km 2 ). Like RF, XGBoost can also handle complex nonlinearity relationships and address the multi-collinearity effect with high accuracy [89]. However, presently, it is hard for RF to process with tens of millions of training data by using existing RF packages, while XGBoost combines software and hardware optimization techniques to make it feasible and easier to study such large-scale study areas (i.e., China). Besides, the solar radiation factor variable, as the unique linking variable in Zhao et al. [58] and Zhao et al. [26], only makes a small contribution to the built models and is limited to the daytime as well. In this study, the variable importance analysis demonstrated the rationality of variable selection for the two soil-related linking variables due to their great contributions to the proposed XGBoost models.

Advantages and Limitations of the Proposed Method
TIR-based LST products are often contaminated by frequent cloud cover with missing pixels in large areas. Thus, we develop an XGBoost-based model to generate spatially continuous MODIS Aqua-Terra LST datasets over a large study area. First, the reconstructed results demonstrated the applicability and effectiveness of the XGBoost linking model to recover both daytime and nighttime LSTs. This reconstructed model is convenient to implement because it does not need much complex and prior parameter knowledge. Second, the XGBoost method had a strong fitting and predictive ability as no obvious fluctuation phenomenon existed in the reconstructed LST values. The proposed model can accurately restore LST with some fine object features such as water bodies (lakes and rivers). Third, the variable importance analysis also revealed the validity for the use of the two instantaneous soil-related variables (i.e., SST and SM) as the linking predictors between cloud-free LSTs and cloud-covered LSTs. The usage of them can facilitate the consideration of soil-related variables when restoring the TIR-based LST parameters.
However, there are still some limitations in this study. First, the preprocessing step for auxiliary variables had some inherent uncertainty. The disaggregation of 8-day MODIS-related auxiliary variables was performed using linear decomposition. Differences exist between the disaggregated results and their true values, though the differences may be ignored. In addition, the soil-related predictors were downscaled by a bilinear interpolation method, which led to relatively smooth results of reconstructed LSTs in areas where a large number of pixels were missing. Second, the SG filter processing cannot fully address the underestimation in the original MODIS LST data. The final reconstructed LST data still have some underestimated values. Third, the spatial scale inconsistency between in situ measured LST data and MODIS LST pixels with a 1 km resolution was not considered in the verification. The in situ measured LST data were required within a very small area while satellite measurements are the integrated responses over a heterogeneous surface. It may cause large uncertainty when the scale inconsistency is ignored in heterogeneous areas.

Conclusions
The quality of TIR-based LST data were greatly affected by frequent cloud cover, which severely restricted their applications in many related scientific fields. Therefore, to explore a practical method for the reconstruction of missing data in TIR-based LSTs is of great significance. In this study, an advanced method was adopted to reconstruct the missing information resulting from cloud cover in MODIS Aqua-Terra daytime and nighttime LST products. The instantaneous solar radiation and soil-related predictors obtained from CLDAS data were used as the linking variables to establish the relationship with instantaneous cloud-free LST data. By coupling with other predictors, we developed an XGBoost-based linking model to fit the complex relationship between cloud-free LST and the predictors, and then the constructed linking model was applied to the areas with cloud-covered pixels to generate spatially continuous MODIS LST estimates during summer time (June to August) in 2017 and 2018 across China. The SG filter method was adopted to further improve the quality of the reconstructed LST products.
The reconstructed LST products were validated and compared with daily in situ measured LST observations and CLDAS LST data. The validation results based on CMA stations showed that the nighttime and daily average LSTs agreed well with in situ data, while the daytime LST estimates showed poorer correlation due to the measurement error for platinum resistance temperature sensor. Evaluation for 6 different land-cover types led to similar conclusions. Besides, the validation results based on in situ LST measurements at five sites in the Heihe River Basin showed that the RMSE values ranged from 3.91 K to 5.53 K for the cloud-free LST and from 4.42 K to 4.97 K for the LST under cloud conditions. The accuracy of the reconstructed LST products meets the requirement comparing with other studies. The comparison with CLDAS LST data also showed similar LST spatial patterns, despite some difference in LST magnitudes.
The variable importance analysis revealed that the soil-related predictors (i.e., soil surface temperature and soil moisture) and the DEM variable made great contributions to the XGBoost model for both daytime and nighttime observation data, which could be used as a reference for the selection of predictors in the reconstruction of TIR-based LST data.