**1. Introduction**

The land surface temperature (LST), generally defined as the radiative skin temperature of the ground, is closely related to the radiative budget and energy fluxes between the atmosphere and the ground [1–4]. LST plays an important role in the estimation of climate models, environmental models and evapotranspiration models, as well as the calculation of drought indices, soil moisture contents and mortality rates [5–14].

Compared with LST measurements at ground stations, satellite remote sensing observations have the advantages of easy acquisition and complete spatial coverage over large areas. Typical LST products include Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Moderate Resolution Imaging Spectroradiometer (MODIS) and Meteosat Second Generation Spinning

Enhanced Visible and Infrared Imager (MSG-SEVIRI) datasets, with spatial resolutions of 90 m, 1 km and 3 km, respectively [15–17]. Among them, the MODIS LST product is the most widely used and best suited for our research because of its appropriate spatial resolution (1 km), high temporal resolution (four overpasses per day), wide coverage (globe), and high retrieval accuracy (approximately better than 1 K). The MODIS instruments were launched on the Sun-synchronous satellites of Terra and Aqua in December 1999 and May 2002, respectively [18]. The MODIS LST products are generated with bands 31 and 32 of MODIS's 36 spectral using the split window algorithm [2]. The latest version C6 MODIS LST products have different spatial resolutions of 1 km, 6 km, and 0.05◦ and different temporal resolutions of daily, eight days and monthly. The MOD11A1 and MYD11A1 are 1 km daily Level 3 products in those MODIS LST products. The transit time of Terra corresponding to MOD11A1 is about 10:30 (22:30), while the transit time of Aqua corresponding to MYD11A1 is about 13:30 (1:30). They are both processed into sinusoidal projection and stored in tiles containing 1200 rows and 1200 columns. The quality of MODIS products is continuing to improve, from more than 2 K in the previous versions to less than 2 K (within ±1 K in most cases) in the C6 version [4]. MODIS LST products have been widely used in LST research [19–21]. However, the MODIS LST product can only provide usable values under clear-sky conditions, and its spatial integrity is thus affected by clouds or other atmospheric disturbances. Taking China as an example, more than half of the pixels per day have no observations on average, and these gaps seriously hinder the application of the MODIS LST product.

Several gap-filling methods have been developed to reconstruct LSTs under cloudy conditions to obtain spatiotemporally-continuous LST products. In general, these methods can be divided into two main groups: clear-sky LST [19–30] and cloudy-sky LST [31–37] methods. Clear-sky LST represents the retrieved LST assuming no cloud effects, whereas cloudy-sky LST represents the actual LST of the reconstruction considering cloud effects. Usually, clear-sky LST is slightly higher than cloudy-sky LST. The methods for reconstructing cloudy-sky LST, mostly based on surface energy balances, often use passive microwave remote sensing data or require ground station measurements or shortwave radiation products. Nonetheless, microwave data have a coarse spatial resolution and an accuracy that needs improvement. Moreover, ground station measurements or shortwave radiation products with a high spatial and temporal resolution are difficult to obtain. This study focused on reconstructing the clear-sky LST, first, because improving the accuracy of clear-sky LST is conducive to further determining the cloudy-sky LST better, and second, because clear-sky LST can be directly applied to research fields such as numerical weather prediction [38], the identification of diurnal patterns of urban heat islands [39], and calculation of the Temperature-Vegetation Dryness Index (TVDI) or Temperature-Vegetation-soil Moisture Dryness Index (TVMDI) [40,41].

The methods for clear-sky LST may be divided into four categories, according to the underlying principles: considering temporal correlation, considering spatial correlation, considering auxiliary information, and the hybrid method. Details of each category are as follows: (1) LST has a temporal correlation because, for the same pixel at different times, the surface properties are the same and only different weather factors, such as solar radiation and wind speed, cause LST differences. Therefore, the first category reconstructs LST based on the temporal correlation using temporal interpolation methods or methods that employ correlations at different times [22,23]; (2) LST also has a spatial correlation because different pixels at the same time have the same weather factors and different surface properties (such as elevation and land cover), but only the surface properties cause LST differences. The second category thus reconstructs LST based on the spatial correlation using spatial interpolation methods [19,24]; (3) in addition, LST is affected by related factors such as elevation and NDVI. The third category thus estimates the missing LST using the empirical relationship between LST and the auxiliary information, which has a similar spatiotemporal resolution to LST and a better spatial coverage integrity than LST [20,21]; (4) finally, the fourth category is hybrid methods that combine two or three of the above methods, such as spatiotemporal gap-filling methods or spatial interpolation methods that consider auxiliary information [25–30]. In general, the hybrid approach is the most promising. Considering only temporal correlation is not suitable for regions with high spatial heterogeneity. If only

spatial correlation is taken into account, the results will be inaccurate for areas that have large weather changes in a short period of time. If only auxiliary information is considered, the accuracy of regression and the uncertainty of auxiliary information will affect the final results. In previous studies, there have been relatively few methods suitable for LST reconstruction in large areas. In a region as large as China, where climate change is complex, spatial heterogeneity is high, and auxiliary information has considerable uncertainty, a method is needed that can comprehensively and reasonably consider time correlation, spatial correlation, and auxiliary information, and the uncertainty of auxiliary information should also be considered.

Bayesian maximum entropy (BME) is a spatiotemporal statistical method proposed by Christakos that can provide a systematic and rigorous framework for incorporating hard data, soft data and other sources of information into the estimation of variables [42,43]. BME has several attractive features; it does not need to make any assumptions regarding the linearity of the estimator, the normality of the underlying probability laws, or the homogeneity of the spatial distribution. Moreover, BME is capable of considering uncertainties contained in the data. The method has been successfully applied to numerous areas, such as air pollution, soil properties, water demand and disease [44–61]. It has also achieved good results in the gap-filling of remote sensing data [62–64]. The BME method is suitable for our research because it can not only take advantage of the temporal correlation and spatial correlation of the LST, but can also explicitly consider the uncertainties of the auxiliary information.

In this study, we applied the BME-based interpolation method to reconstruct 1 km resolution daily clear-sky LST for China's landmass considering temporal correlation, spatial correlation and auxiliary information. The goals of this article are to (1) examine the feasibility of the BME method to reconstruct LST for the whole of China, (2) discuss the accuracy of the BME method for different land cover types, and (3) compare the BME method with other commonly used LST reconstruction methods.

#### **2. Materials and Methods**

#### *2.1. Study Area*

In this study, we selected the land area of 34 provinces in China as our study area. China is located on the eastern side of the Eurasian continent and the western shore of the Pacific Ocean; it spans approximately 5500 km from north to south and 5000 km from west to east. The topography across China is complicated and includes plains, plateaus, mountains, hills and basins. It varies from the Qinghai-Tibet Plateau at more than 4000 m above sea level (peaking at 8848 m) to its eastern coastline on the Pacific Ocean. China's land resources are vast, and its use types are diverse. The cultivated lands are mainly distributed in the eastern region, the forests are distributed in the south and northeast regions, the grasslands are mainly distributed in the central and southwestern regions, and the unused lands are mainly distributed in the northwestern region. China's climate is governed by monsoonal circulations, and winters with low temperatures and little rain significantly differ from summers with high temperatures and abundant rain [65].

#### *2.2. Data Acquisition and Preprocessing*

LST can be affected by many factors [20,66,67]. In view of relatively more critical factors across the Chinese scale and the convenience of data acquisition and processing, NDVI, DEM, longitude and latitude were selected as auxiliary data to regress LST. The following specific datasets were used in this study: (1) For LST, we used the MODIS/Aqua LST Daily L3 Global 1 km SIN Grid product (MYD11A1, Collection 6). LSTs observed throughout the year 2015 were used at local 1:30/13:30 overpass times, which approximate daily minimum and maximum LST values; in the later part of the article, we call them daytime LST and nighttime LST; (2) for NDVI, we selected the MODIS/Aqua 16 day 1 km Vegetation Index product (MYD13A2, Collection6), which has great spatial completeness and the spatial resolution we need. There are 23 NDVI data sets of MYD13A2 for the entire year of 2015; (3) for DEM, we used the Shuttle Radar Topography Mission (SRTM) Digital Elevation Data Version 4 at a 90 m

spatial resolution produced to provide consistent, high-quality elevation data. The original DEM data were resized using nearest-neighbours from 90 m to 1 km; (4) for land cover data, we used the MODIS Land Cover Type Yearly Global 500 m product (MCD12Q1) from 2015, which was derived using supervised classifications of MODIS Terra and Aqua reflectance data. We combined all the land types into six categories: forests, grasslands, croplands, barren, urban and snow. Evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, mixed forests, closed shrublands and open shrublands were merged into forests, whilst woody savannas, savannas and grasslands were merged into grasslands. Water was not considered in this study. The pixels were resampled to a 1 km resolution in preparation for the subsequent verification phase; (5) for longitude and latitude, we gave each grid one longitude value and one latitude value at a 1 km spatial resolution based on the WGS84 datum.

MYD11A1, MYD13A2 and MCD12Q1 were provided by the Land Processes Distributed Active Archive Center (LP DAAC) site (https://lpdaac.usgs.gov/). DEM datasets are available from the CGIAR-CSI SRTM 90 m Database site (http://srtm.csi.cgiar.org). In this study, all the above data were downloaded, reprojected, stitched and resized with Google Earth Engine (GEE, https://earthengine. google.com). We processed all the data into 3540 rows × 6166 columns to cover the study area.

#### *2.3. Method*

As shown in Figure 1, BME was the core method of this study, and data were prepared for adapting the BME procedure. The three aspects of considering auxiliary data, time correlation, and spatial correlation were used to describe how the three most important input parameters of the BME algorithm, hard data, soft data and covariance models, were constructed. Regarding the auxiliary data, for each image to be estimated, the pixel LSTs with MODIS observations were taken as the dependent variable, and the NDVI image that was taken on the date closest to the estimated image and the elevation, longitude and latitude of the corresponding pixels were taken as independent variables to perform multiple linear regression and obtain the regression coefficients. Then, via the four independent variables of the pixels to be estimated and the respective regression coefficients, the regression LSTs of all the pixels to be estimated were calculated. The regression LSTs were used as the mean value and the mean square error (MSE) between the dependent variable and the predicted value as the variance to construct the Gaussian LST distribution as the soft data. The calculation results of the regression coefficients and the regression R<sup>2</sup> of the multiple linear regressions for the daytime and nighttime are shown in Tables A1 and A2. Regarding time correlation, we subtracted the mean LSTs of 15 days (7 days before and after) from both MODIS observed LSTs used as hard data and regression LSTs used as soft data, and the resulting difference values were input into the BME model as real hard data and soft data, respectively. This is equivalent to the 15 day mean LST minus the LST of the estimated day, which can be understood as the average trend after removing the special weather factors of the day; the LST time correlation could then be taken into account by this simple calculation. After such processing, the LSTs did not need to have trends removed before being input into the BME model, and only the spatial correlation needed to be considered, which could be achieved by the covariance function that represents the spatial dependence. The specific cause analysis can be found in lines 4 to 11 of the fourth paragraph of the introduction. Regarding spatial correlation, we calculated the spatial covariance using the real hard data and soft data mentioned above and input the obtained spatial covariance function name and parameters into the BME model. The relevant parameters of the covariance model calculated in this paper are shown in Tables A3 and A4.

From the above, it is worth noting that considering the availability of data and the simplicity of method processing, this study made the following assumptions when using the BME method to reconstruct LST: (1) LST changed linearly in a short time (time correlation was considered by subtracting the 15 day mean LST from the LST of the day to be estimated); (2) for the day to be reconstructed, one omnidirectional covariance model of that day can be used in the whole study area; (3) NDVI of each day can be represented by the NDVI data of its adjacent 7 days (the time resolution of the NDVI

data was 15 days, so the time interval between the selected NDVI and the estimated image on any day was 0 to 7 days). The main BME conceptual core and framework and the explanations for its use in this research, are shown below.

**Figure 1.** Flowchart describing the land surface temperature (LST) reconstruction model using the Bayesian maximum entropy (BME) method.

#### 200 2.3.1. BME Epistemic Paradigm and Conceptual Core

ீ

201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 The BME approach belongs to modern geostatistics, which provide insights into spatiotemporal variables. The epistemic paradigm of BME distinguishes between three main stages of knowledge acquisition, interpretation, and processing, as follows: (1) The prior stage. Spatiotemporal analysis and mapping always starts with a basic set of assumptions and the general knowledge base G. G refers to the background knowledge and the justified beliefs relative to the overall mapping situation; (2) the meta-prior stage. The specific knowledge base S is considered, including hard and soft data. S refers to a particular occurrence or state of affairs at a particular location and at a particular time. Hard data are considered accurate or have a high degree of confidence. Soft data are uncertain observations expressed in terms of interval values, probability statements, empirical charts, and others. That is to say, soft data can have varying levels of uncertainty and may be derived from the direct calculation of the probabilities or the indirect estimation from accumulated experience; (3) the integration or posterior stage. Information from (1) and (2) is processed by means of logical rules to produce the required spatiotemporal map. Therefore, the conceptual core of the BME method is that it aims at informativeness (in terms of prior information relative to the general knowledge G), as well as cogency (in terms of posterior probability relative to the specific knowledge S). BME combines the maximum entropy theory with operational Bayesian statistics to construct its scientific mathematical framework and to implement the above conceptual heart. In general, BME is used to acquire various knowledge bases and to order these bases in an appropriate manner so that, when taken together, they form a realistic picture of the phenomenon of interest.

#### 220 2.3.2. BME Framework

198 199

221 222 223 BME has a rigorous cognitive system and a mathematical reasoning framework. The complete theoretical basis, mathematical formulas and specific derivation processes can be found in reference [43]. The main BME formulas and steps involved in this study are as follows:

121

(χ, χௗ௧) <sup>=</sup> ீ൫χ൯ = ୣ୶୮൫∑ <sup>ఓ</sup>ഀ൫ೌ൯ഀ൫ೌ<sup>൯</sup>

ಿ

ഀసభ ൯ ௗೌ ୣ୶୮൫∑ ఓഀ൫ೌ൯ഀ൫ೌ൯ ಿ

ഀసభ ൯

$$f\_{\mathbb{G}}(\chi\_{\mathbb{k}}, \chi\_{\text{data}}) = f\_{\mathbb{G}}(\chi\_{\text{map}}) = \frac{\exp\Big(\sum\_{a=1}^{N} \mu\_{a}(p\_{\text{map}}) g\_{a}(\chi\_{\text{map}})\Big)}{\int d\chi\_{\text{map}} \exp\Big(\sum\_{a=1}^{N} \mu\_{a}(p\_{\text{map}}) g\_{a}(\chi\_{\text{map}})\Big)} \tag{1}$$

where χ*<sup>k</sup>* denotes the LST of the estimated pixel, χ*data* = (χ*hard*, χ*so f t*), χ*hard* represents hard data, and χ*so f t* represents soft data. In this study, MODIS LST observations were used as hard data and LST Gaussian distributions obtained by multivariate linear regression were used as soft data. The regression process is described above. *fG*(χ*<sup>k</sup>* , χ*data*) denotes the prior pdf of the map χ*map* = (χ*<sup>k</sup>* , χ*data*) given the general knowledge base G. µα(*pmap*) represents Lagrange multipliers. *g*α(χ*map*) is a set of known functions of χ*map*. In practical applications, prior knowledge usually includes the first-order statistical moment (mean trend) and second-order statistical moment (covariance). The mean trend was not adopted as the first-order statistical moment in this study; rather, the LST difference on the observation day minus the 15 day mean was used. This was done to take time correlation into account and to remove LST instability caused by different weather factors, which resulted in a more stable LST distribution. The second-order statistical moments employed in this study were the spatial covariance functions derived by the difference values calculated above. Such a priori knowledge in this study could consider both the temporal correlation dominated by weather factors and the spatial correlation dominated by surface properties.

$$f\_{\mathcal{K}}(\chi\_k) = f\_{\mathcal{G}}(\chi\_k|\chi\_{\text{hard}}, \chi\_{\text{soft}}) = f\_{\mathcal{G}}(\chi\_k|\chi\_{\text{data}}) = f\_{\mathcal{G}}(\chi\_k, \chi\_{\text{data}}) / f\_{\mathcal{G}}(\chi\_{\text{data}}) \tag{2}$$

In Equation (2), *f<sup>K</sup>* denotes the posterior pdf of the map χ*<sup>k</sup>* , given the total knowledge base K comprised of general knowledge G and specific knowledge S, including hard and soft data. The general knowledge, hard and soft data used in this study were described earlier.

$$\chi\_{k\prime}\text{ mean} = \int \chi\_k f \chi(\chi\_k) d\chi\_k \tag{3}$$

We used the BME mean value (χ*<sup>k</sup>* , mean) as the final estimated LST. The BME mean value could be calculated from the posterior PDF since we sought to penalize large errors more than smaller ones.

### 2.3.3. BME Implementation

We used the BMElib algorithm package for BME algorithm implementation in MATLAB [68]. The calculation details for each estimated day are as follows. Firstly, the mean value and mean square error of the soft data of each prediction point were input into the probaGaussian.m function of the software package to obtain the soft data information that meets the requirements of the subsequent input. Secondly, the values and position coordinates of hard and soft data were entered into the covario.m and corefit.m functions of the software package to obtain the covariance function name and parameters. Finally, the main function was used to calculate the final result. Information such as hard data, soft data, and covariance was input into the BMEprobaMoments.m function. In addition, the maximum effective distance was set to 15 km, the maximum hard data point was set to 20 points, and the maximum soft data point was set to 3 points.

#### **3. Results**

### *3.1. Spatial Patterns of the Reconstructed LSTs*

We selected the 15th day of each month in 2015 to conduct the method experiment and obtained the spatial distribution results of 12 images for both the daytime and nighttime (Figures 2 and 3). In the daytime, an average of 43% of the pixels of the MODIS LST products in the study area had LST observations, and in the nighttime, the value was 51%. That is, there was a missing rate of nearly one-half before filling gaps. The missing LST could be 100% filled using the BME method to generate a complete spatial distribution (Table 1).


**Table 1.** Availability of Moderate Resolution Imaging Spectroradiometer (MODIS) observed LST and reconstructed LST for the daytime and nighttime from 15 January to 15 December in 2015.

In general, the entire study area showed strong spatial heterogeneity that varied in different seasons for both the day and night (Figures 2 and 3).

In winter, the lowest LST for the daytime occurred in northeast China, followed by the Qinghai-Tibet Plateau, whereas the lowest LST for the nighttime appeared in the Qinghai-Tibet Plateau, followed by northeast China. In summer, the lowest LST during the daytime simultaneously occurred in the Qinghai-Tibet Plateau and northeast China, whereas during the nighttime, the LST of the Tibetan Plateau was significantly lower than that in northeast China (Figures 2 and 3). The LST of the Qinghai-Tibet Plateau in southwest China was obviously low due to its high topography, and the LST of northeast China affected by Siberian cold air in winter was also low.

281 **Figure 2.** Spatial distribution of reconstructed daytime LST from 15 January to 15 December in 2015.

**Figure 3.** Spatial distribution of reconstructed nighttime LST from 15 January to 15 December in 2015.

284 285 286 287 288 289 290 In summer, the highest daytime LST was evident in northwest and central Inner Mongolia, whereas the highest nighttime LST was widely distributed in the south-central and southeastern regions, except for the Qinghai-Tibet Plateau. In winter, the highest daytime LST was distributed in the northwest and southeast regions, but the highest nighttime LST was distributed in the southeast coastal areas (Figures 2 and 3). The LST was usually higher in the northwest and Inner Mongolia due to the large number of deserts. LST was also higher in the southern region because of its low latitude and more solar radiation on the ground.

#### 291 292 *3.2. Accuracy Assessment*

282 283

293 294 295 296 297 298 299 300 301 302 303 304 Since the study aimed at clear-sky LST reconstruction, it was not necessary to employ ground observation points for verification. First, this is because clear-sky LST is the theoretical value that is assumed not to be affected by clouds, while the ground observed value is the real value that is affected by clouds, so they cannot be directly compared. Secondly, this is because the acquisition time of the ground station is difficult to coincide with that of MODIS products, and there are also scale effects between points and surfaces. The verification method in this paper is for clear-sky LST. We selected some points with MODIS LST observations to cover them, reconstructed the LST of the covered points with the BME method, and then compared the reconstructed LST values with the known observations of MODIS LST. The verification points must have MODIS LST observations as references for the reconstructed LSTs. Therefore, we selected the points where the MODIS LST observations existed on the 15th of each month in 2015 as the verification points, which also helped to show the accuracy change of points with the same positions over time (Figure 4). There were 10,971 test pixels for the daytime (green points in Figure 4), including 330 forest pixels, 2683 grassland pixels, 537 cropland pixels, 7376 barren pixels, 38 urban pixels and 7 snow and ice pixels. There were 14,376 test pixels for the nighttime (blue points in Figure 4), including 218 forest pixels, 6040 grassland pixels, 425 cropland pixels, 7627 barren pixels, 425 urban pixels and 8 snow and ice pixels. Some of the verification points

305 306

were spatially continuous and they formed regions of various shapes. The maximum diameter of the regions formed by the verification pixels was close to 60 km, and the estimation accuracy was thus also of reference value for the missing LST values caused by large cloud cover.

**Figure 4.** Distribution of 10,971 verification points in the daytime and 14,376 verification points in the nighttime, and four verification points in big cities.

The daytime accuracy was mostly lower than that of the nighttime (Figures 5 and 6). The average mean absolute error (MAE) and RMSE values were 1.1 K and 1.6 K, respectively, in the daytime, whereas the values were 0.8 K and 1.2 K in the nighttime, respectively (Figure 7d–f). The accuracy in summer was generally lower than in winter, which decreased and then increased from January to December (Figures 5, 6 and 7d–f). During the daytime, the maximum RMSE was 3.0 K in July and the minimum was 0.8 K in January. During the nighttime, the maximum RMSE was 1.9 K in June and the minimum was 1.0 K in December (Figures 5 and 6). The higher RMSE in summer than in winter may be due to the higher surface heterogeneity in summer than in winter, and the higher RMSE during the day than during the night may be due to more serious cloud cover and more missing values during the day.

In general, the RMSEs of barren land were the largest, with averages of 1.6 K and 1.4 K during the day and night, respectively, whereas the RMSEs of urban areas were the smallest, with averages of 1.0 K and 0.6 K during the day and night, respectively (Figure 7a,b). During the daytime, RMSEs ranked from high to low for barren, grasslands, forests, croplands, urban and snow and ice. During the nighttime, RMSE was in the order of forests, barren, snow and ice, grasslands, croplands and urban (Figure 7c).

As shown in Figures 5 and 6, the changes between day and night of urban LST were the most obvious, and the seasonal changes of ice and snow LST were most obvious. For the daytime, urban LST was close to the average LST of different land cover types, whereas for the nighttime, urban LST was generally higher than the average LST of different land cover types. This indicates that the urban areas have a relatively strong heat island effect at night. The LST of snow and ice was lower in winter and higher in summer.

329 330

331 332

**Figure 5.** Scatter plots of reconstructed LST versus observed LST for 10,972 pixels for the daytime from 15 January to 15 December in 2015.

**Figure 6.** Scatter plots of reconstructed LST versus observed LST for 14,376 pixels for the nighttime from 15 January to 15 December in 2015.

**Figure 7.** RMSE of (**a**) daytime LST and (**b**) nighttime LST for each land cover type from 15 January to 15 December in 2015; (**c**) overall average RMSE for each land cover type in 2015; (**d**) overall average mean absolute error (MAE) from 15 January to 15 December in 2015; bias and overall average RMSE for the daytime (**e**) and nighttime (**f**) from 15 January to 15 December in 2015.

337 338 339 340 341 342 343 344 345 346 347 348 349 In addition, we selected one point in Beijing, Wuhan, Shanghai and Guangzhou to estimate the LST for 365 days in 2015 and validated the accuracy with MODIS observations. These four verification points were geographically discrete and located in four major cities of China from north to south (Figure 4). The results are shown in Figure 8. Beijing had fewer than 200 days with LST observations, whereas in the other three cities, the number was less than 100. The maximum number of consecutive missing days was 30 days in Beijing, 40 days in Wuhan, 44 days in Shanghai and 54 days in Guangzhou. BME could fill 100% of the missing LSTs and reconstruct the uninterrupted LST time curve of each pixel for 365 days. As seen from the variation range of the curve, the BME method can describe the change of LST in a relatively fine manner, without smoothing out the maximum and minimum values. The R<sup>2</sup> values of the four urban test sites were all greater than 99%, and, except for the RMSE in Wuhan of 0.6 K, the RMSE in the other cities in the day or night was less than 0.5 K. Therefore, the single point test accuracy of the BME method was very high in large cities. The time distribution of LST demonstrated that there were obvious temperature differences between the day and night in the four cities. The four seasons changed most obviously in Beijing because it is located in a typical

353 ate.

north temperate semi-humid continental monsoon climate zone, whereas Guangzhou had the smallest difference between the four seasons because of the Marine subtropical monsoon climate. 350 351 352

**Figure 8.** Temporal pattern of observed LST and reconstructed LST for the (**a**) daytime in Beijing, (**b**) nighttime in Beijing, (**c**) daytime in Wuhan, (**d**) nighttime in Wuhan, (**e**) daytime in Shanghai, (**f**) nighttime in Shanghai, (**g**) daytime in Guangzhou, and (**h**) nighttime in Guangzhou.

356

354

374 375

375

376

376 377

#### *3.3. Factors that Influence Accuracy* 357

Figures 9 and 10 illustrate the influence of different factors on the LST estimation accuracy represented by RMSE. For the daytime, the Pearson correlation coefficients between the RMSE of the reconstructed LST and the four influencing factors multiple linear regression R<sup>2</sup> , average LST, ratio of pixels with LST observations to total pixels (namely, completeness) and average NDVI were 0.32, 0.85, −0.21 and 0.86, respectively; for the nighttime, the corresponding Pearson correlation coefficients were 0.22, 0.55, −0.6 and 0.44, respectively. Therefore, the average temperature and average NDVI were strongly correlated with RMSE in the daytime, and the mean LST and completeness were moderately correlated with RMSE in the nighttime. Three rather interesting aspects emerged from the results: (1) The average LST affected the accuracy of the method, where the higher the average temperature, the larger the RMSE; (2) the completeness, or the number of missing pixels, slightly affected the accuracy of the method; (3) the accuracy R<sup>2</sup> of multiple linear regression did not affect the accuracy of the method. In this study, R<sup>2</sup> varied from 0.39 to 0.90 (from 0.39 to 0.82 during the daytime and from 0.79 to 0.90 during the nighttime). This suggests that the BME method has a great ability to consider the uncertainty of soft data. Since the BME method does not have high requirements for the accuracy of soft data, it can be applied to other large-scale regions, and there is no need to improve the regression accuracy by random forest or other regression methods. 357 358 359 360 361 362 − 363 − 364 365 366 367 368 369 370 371 372 373 358 359 360 361 362 <sup>−</sup> 363 − 364 365 366 367 368 369 370 371 372 373 374

**Figure 9.** Correlation between multiple regressive R<sup>2</sup> , mean LST and RMSE for the (**a**) daytime and (**b**) nighttime.

**Figure 10.** Correlation between the completeness of the observed LST, mean NDVI and RMSE for the (**a**) daytime and (**b**) nighttime.

#### 377 378 *3.4. Comparisons with Other Methods*

378 We compared the BME method with four other commonly used LST gap-filling methods, including Crosson's method of supplementing MYD data with the same day's MOD data [22], the time

interpolation method HANTS [23], the Kriging spatial interpolation method, and the hybrid gap-filling method proposed by Li [30]. It is worth noting that the same hard data and the same spatial covariance model of the BME method were entered into Kriging, and the only difference was that the Kriging method does not consider soft data. The RMSEs and error distributions of each method are shown in Figure 11. In general, the accuracy of each method ranked from high to low, as follows: BME > Kriging > Hybrid > HANTS > Crosson. It appeared that Crosson's method had the lowest accuracy, the BME method had the highest accuracy during the daytime, the Kriging method had the highest accuracy during the nighttime, the Hybrid method had a stable accuracy in the day and night, and HANTS had a significantly higher accuracy in the night than in the day.

**Figure 11.** Error distribution of LST using the five methods of Crosson, HANTS, Kriging, Hybrid and BME for the (**a**) daytime on 15 January, 2015; (**b**) nighttime on 15 January, 2015; (**c**) daytime on 15 July, 2015; and (**d**) nighttime On 15 July, 2015.

#### **4. Discussion**

392 393 394

395

#### *4.1. Accuracy Analysis*

The mean RMSEs in this study were 1.6 K for the daytime and 1.2 K for the nighttime, which were slightly lower than the RMSEs of 3.3 K for the daytime and 2.7 K for the nighttime in Li's study in a comparably large area [30]. The single point RMSE of approximately 0.5 K is comparable with the RMSE of approximately 2 K under cloud-free conditions in Duan's study, which selected four ground points for validation [31]. The accuracy of this method is acceptable for large areas with complex geographical and climatic conditions. In addition, this method has a high accuracy in estimating urban LST and can be applied to urban heat island research.

There were different accuracies for different land cover types, which indicates that the accuracy is affected by land cover [34]. The accuracy of barren land and forest was lower than that of urban and cropland because the terrain of barren and forest is more complex, and the spatial heterogeneity is greater. Therefore, the model accuracy can be improved by dividing various terrain regions and then adopting different covariance models for various regions.

The accuracy decreased with the increase of NDVI and average LST because a high NDVI and high temperature usually represent the summer climate in China, when the cloud cover is large and the distribution is concentrated, which results in large LST gaps. The completeness has only a slight impact on the accuracy, possibly because the accuracy is influenced not only by the number of missing pixels, but also by their maximum diameter and distribution characteristics [69].

When constructing soft data, the accuracy of multiple linear regression does not affect the accuracy; this may be due to the ability of BME to fully consider the uncertainty of soft data. The average regression R<sup>2</sup> in this study was nearly 0.6, and in the subsequent application of the BME method to LST reconstruction, when the average regression R<sup>2</sup> in the construction process of soft data is greater than 0.6, it is unnecessary to adopt more complex regression methods to improve the regression accuracy.

#### *4.2. Suggestions for Method Selection*

We can learn the characteristics of each LST reconstruction method from Figure 11 and in combination with previous studies. The accuracy of spatial interpolation models is usually higher than that of temporal interpolation models [70]. The time interpolation models have some difficulties in capturing extreme values, and their accuracy is relatively low. Spatiotemporal gap-filling methods are often unable to fill all the missing values at one time, and one usually needs to iterate several times until all the missing values are filled. Spatial interpolation methods, especially the ones that consider auxiliary information, have a high accuracy, but usually take a relatively long time for calculations [25,26]. The study area in this study was large. To balance the calculation time and accuracy, we did not select the spatiotemporal covariance model, but rather the spatial covariance model, to consider the spatial dependence characteristics and the simple 15 day mean value to consider the time dependence characteristics. The spatiotemporal covariance model can be selected in small areas [34,71]. With the development of computing power and multi-core parallelism in the future, the computing speed will become faster.

Suggestions on how to choose an appropriate method to reconstruct LST are as follows: (1) If one hopes for a short computation time and simple computation steps, one can choose the time interpolation method; (2) if a high precision and simple calculation steps are required, the spatial interpolation method is recommended; (3) if one wants to balance the calculation speed and accuracy, we suggest using the spatiotemporal gap-filling method. All of the above three methods can consider introducing auxiliary data to improve the accuracy. In addition, note that the Kriging spatial interpolation method achieved good results in this study area and that Kriging is thus a simple method worth trying. The BME method that we used is a spatial interpolation method that considers auxiliary information, and its precision was very high in this study area.

#### *4.3. Exploration of the Accuracy Improvement*

In this part, we will explore methods to improve the accuracy and model applicability based on the assumptions of the method. One assumption of this method was that LST changed linearly in a short time. We thus subtracted the 15 day mean LST from the LST of the constructed day to consider the time correlation. Doing so can also make the data closer to a normal distribution and thus replace the step of trend removal in other BME studies. We selected 18 October, 2015 as an example (there were relatively more MODIS observed LSTs available on that day), as shown in Figure 12. The green part shows that the LST of all pixels on this day presents a non-normal distribution, while the purple part exhibits an approximately normal distribution after subtraction calculation (LST–the mean LST of 15 days) is performed. This operation achieved the desired effect. However, as can be seen from Figure 8, the time

461 462

469

470

variation of LST exhibited both an overall trend and fluctuation. According to the limitation of this study's assumption, we may improve the accuracy by taking better account of the time correlation. We can do so by introducing the Annual Temperature Cycle (ATC) model, which is a general and smooth curve description of the LST annual change. Considering the time correlation with the LST of the estimated day minus the LST of the corresponding point on the ATC curve, theoretically, will be more accurate than considering the time correlation with the LST of the estimated day minus the LST of the 15 day mean value. We plan to conduct follow-up studies with regards to this.

**Figure 12.** Data distribution of the MODIS observed LST and the difference of LST minus the 15 day mean LST on 18 October, 2015.

The second assumption of this method was that one omnidirectional covariance model of the constructed day can be used in the whole study area. Here, we want to explore the directivity of the covariance model. We constructed omnidirectional and directional covariance models for 18 October, 2015 (Table 2, Figure 13). As can be seen from the parameters of directional covariance in the study area, the overall characteristic is that the directional covariance is significantly affected by longitude and latitude, and the latitude direction changes faster than the longitude direction.


**Table 2.** Omnidirectional and directional covariance model parameters for 18 October, 2015.

**Figure 13.** Directional covariance models of the daytime and nighttime on 18 October, 2015.

Omnidirectional and directional covariance models were input into the method in this study to calculate the results, and 10,000 points were randomly selected during the day and night for accuracy verification and comparison. The results show that after considering the directionality of the covariance, the accuracy in the daytime is slightly improved, and the accuracy of the night is not improved (Table 3). It is suggested that directionality can be considered during the daytime in the following research.


**Table 3.** Omnidirectional and directional accuracy verification results for 18 October, 2015.

In addition, large study areas will have problems with spatial covariance models that differ in different regions. First of all, however, we have not explored how the covariance model changed in different regions of China's mainland. Moreover, if we want to fill in the LST of the whole study area at once, more detailed regional division may make the model more complicated. How to apply different covariance models in different subregions of the study area is challenging and worthy of further exploration.

The third assumption of this method was that the NDVI of each day can be represented by the NDVI data of its adjacent 7 days. We used 15 day NDVI data because it had values on all pixels to ensure that soft data can be constructed on all LST missing pixels. When the LST was reconstructed on the 15th day of each month in 2015, the nearest NDVI data were selected on 9 January, 10 February, 14 March, 15 April, 17 May, 18 June, 20 July, 21 August, 22 September, 8 October, 8 November and 11 December. There were no results to prove that the closer the reconstruction date was to the obtained date of NDVI, the higher the fitting accuracy and the final LST accuracy were. We believe it was feasible to reconstruct the LST daily with a 15 day NDVI product. Lastly, we cannot accurately calculate the impact of the uncertainty of NDVI datasets on the uncertainty of the final results, which is a limitation of this method.

### *4.4. Recommendations for Future Studies*

Due to the characteristics of the BME method, we cannot accurately determine the uncertainty of the results. At present, we can only obtain the conclusion that at a 1 km spatial resolution, the accuracy of reconstructing the daily LST of China's landmass with the BME method is acceptable.

For daytime or nighttime LST reconstruction on a certain day, one covariance model can be adopted in the whole research area on that day, which can achieve a reasonable accuracy. If one wants to use this model later, we suggest that it be directly used in the same study area and at the same spatiotemporal resolutions. If other areas or other resolutions are studied, an accuracy analysis should be conducted first to see whether it can meet the actual requirements. The soft data should also be reconstructed according to the geographic and data characteristics of the other study areas. One can try to use the four independent variables in this paper or may introduce other auxiliary data, such as soil moisture and temperature. In the future, we can try to introduce the ATC model, divide different subregions and adopt different covariance models to improve the accuracy of the BME method presented in this study.

When using the BMElib package, it is important to be aware of some parameter settings. The maximum effective distance can be set to a value similar to the range in the covariance model. In order to balance the calculation accuracy and calculation time, we recommend that the maximum number of hard data points should not exceed 50 and the maximum number of soft data points should not exceed 5.

Although clear-sky LST can be applied in some research, the real surface LST is still needed in many fields. This study calculated clear-sky LST, and if one wants to obtain cloudy-sky LST from clear-sky LST, one can refer to Zeng's method [69]. Adding microwave and ground observation data can also be considered.

#### **5. Conclusions**

The MODIS LST product has many missing values over wide areas, which hinders its practical application. In this study, we reconstructed the seamless 1 km resolution daily clear-sky LST for China's landmass based on the BME method, considering spatiotemporal correlation and taking auxiliary data as soft data. The average RMSE was 1.6 K for the daytime and 1.2 K for the nighttime, with the mean absolute error (MAE) of 1.1 K for the daytime and 0.8 K for the nighttime, and the corresponding R<sup>2</sup> of 0.92 for the daytime and 0.98 for the nighttime.

This method has the following advantages: (1) It simultaneously considers spatiotemporal correlation and auxiliary data and has a high accuracy in a large area. It has the ability to capture extreme values; (2) the data in this method are easy to obtain and process; (3) simple linear regression is used to construct soft data, and there is no need to adopt more complex regression methods to improve the regression accuracy, as long as the average regression R<sup>2</sup> is greater than 0.6; (4) even if the diameter of the missing area is large or the continuous missing time is long, this method does not need multiple step-by-step calculations to gradually fill in the missing pixels, and can estimate all the missing pixels at one time.

There are also some limitations for this method: (1) This method is not applicable when an accuracy of less than 1 K across the entire Chinese landmass is required; (2) when using the method in other study areas and spatiotemporal scales, it is necessary to first consider whether the hypothesis of LST linearity change in a short time and one omnidirectional covariance model can be applied to the entire study area are valid; (3) the method cannot quantitatively calculate the influence of the uncertainty of NDVI and DEM data on the uncertainty of the results; (4) the clear-sky LST should be converted to cloudy-sky LST if the real LST is required.

The results of this study provide a data basis for daily LST analysis and subsequent relevant studies in large areas of China. For the method, its high accuracy and great ability to capture extreme values in urban areas can help improve urban heat island research. It can also be applied to the reconstruction of missing LST values of other years, other regions and other spatial resolutions (such as Landsat), as well as the estimation of missing values of other geographical attributes.

**Author Contributions:** Conceptualization, Y.Z. and Y.C.; data curation, Y.Z. and Y.L.; formal analysis, Y.Z. and Y.C.; funding acquisition, Y.C. and J.L.; investigation, Y.Z.; methodology, Y.Z.; project administration, Y.C.; software, Y.Z. and Y.L.; supervision, Y.C. and J.L.; validation, Y.Z. and H.X.; visualization, Y.Z.; writing original draft, Y.Z.; writing—review and editing Y.Z., Y.C., Y.L., H.X. and J.L.

**Funding:** This work was supported by the National Key R&D Program on monitoring, early warning and prevention of major natural disasters under Grant 2017YFC1502406; the Natural Science Foundation of China under Grants 41571342, 41771448 and 51579135; the Beijing Natural Science Foundation under Grant 8192025; and in part by the Beijing Laboratory of Water Resources Security.

**Acknowledgments:** We are very grateful to Patrick Bogaert (Université Catholique de Louvain; Belgium) and Marc Serre (University of North Carolina at Chapel Hill; USA) for providing us with the MATLAB algorithm package BMElib and to Alexander Kolovos (SpaceTimeWorks, LLC; USA), Hwa-Lung Yu (National Taiwan University, Taipei; Taiwan), Steve Warmerdam and Boris Dev (San Diego State University, CA, USA) for providing us with the MATLAB visualization algorithm package SEKS-GUI. We would also like to thank the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) of NASA for MODIS data and the Consultative Group of International Agricultural Research-Consortium for Spatial Information (CGIAR-CSI) for SRTM DEM data. We would like to thank the anonymous reviewers for their insightful comments and substantial help in improving this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** Regression R<sup>2</sup> and regression coefficients of multiple linear regression for the daytime.


**Table A2.** Regression R<sup>2</sup> and regression coefficients of multiple linear regression for the nighttime.


**Table A3.** Names and parameters of the spatial covariance model for the daytime.



**Table A4.** Names and parameters of the spatial covariance model for the nighttime.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
