Reconstructing One Kilometre Resolution Daily Clear-Sky LST for China’s Landmass Using the BME Method

Abstract

The land surface temperature (LST) is a key parameter used to characterize the interaction between land and the atmosphere. Therefore, obtaining highly accurate, spatially consistent and temporally continuous LSTs in large areas is the basis of many studies. The Moderate Resolution Imaging Spectroradiometer (MODIS) LST product is commonly used to achieve this. However, it has many missing values caused by clouds and other factors. The current gap-filling methods need to be improved when applied to large areas. In this study, we used the Bayesian maximum entropy (BME) method, which considers spatial and temporal correlation, and takes multiple regression results of the Normalized Difference Vegetation Index (NDVI), Digital Elevation Model (DEM), longitude and latitude as soft data to reconstruct space-complete daily clear-sky LSTs with a 1 km resolution for China’s landmass in 2015. The average Root Mean Square Error (RMSE) of this method was 1.6 K for the daytime and 1.2 K for the nighttime when we simultaneously covered more than 10,000 verification points, including blocks that were continuous in space, and the average RMSE of a single discrete verification point for 365 days was 0.4 K for the daytime and 0.3 K for the nighttime when we covered four discrete points. Urban and snow land cover types have a higher accuracy than forests and grasslands, and the accuracy is higher in winter than in summer. The high accuracy and great ability of this method to capture extreme values in urban areas can help improve urban heat island research. This method can also be extended to other study areas, other time periods, and the estimation of other geographical attribute values. How to effectively convert clear-sky LST into real LST requires further research.

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,2,3,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,6,7,8,9,10,11,12,13,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,16,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,20,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,20,21,22,23,24,25,26,27,28,29,30] and cloudy-sky LST [31,32,33,34,35,36,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,26,27,28,29,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,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61]. It has also achieved good results in the gap-filling of remote sensing data [62,63,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² of the multiple linear regressions for the daytime and nighttime are shown in

Table A1. Regression R² and regression coefficients of multiple linear regression for the daytime.

Date	Regression R2	Intercept	Coefficient of NDVI	Coefficient of DEM	Coefficient of Latitude	Coefficient of Longitude
15 January, 2015	0.69	366.07	−2.10	−0.004	−1.07	−0.39
15 February, 2015	0.58	371.31	−5.23	−0.004	−1.26	−0.31
15 March, 2015	0.50	381.25	−24.29	−0.005	−1.52	−0.17
15 April, 2015	0.40	369.12	−4.05	−0.005	−0.18	−0.50
15 May, 2015	0.59	388.10	−26.26	−0.007	−0.63	−0.39
15 June, 2015	0.59	343.54	−23.48	−0.007	−0.10	−0.13
15 July, 2015	0.82	373.64	−21.95	−0.008	−0.16	−0.35
15 August, 2015	0.59	320.09	−28.35	−0.005	−0.07	0.10
15 September, 2015	0.39	317.19	−14.30	−0.003	−0.02	−0.03
15 October, 2015	0.50	349.92	−16.69	−0.004	−0.52	−0.19
15 November, 2015	0.67	370.45	−16.05	−0.004	−1.38	−0.21
15 December, 2015	0.59	358.56	−8.76	−0.003	−1.65	−0.14

and

Table A2. Regression R² and regression coefficients of multiple linear regression for the nighttime.

Date	Regression R2	Intercept	Coefficient of NDVI	Coefficient of DEM	Coefficient of Latitude	Coefficient of Longitude
15 January, 2015	0.87	337.62	−18.39	−0.005	−0.97	−0.18
15 February, 2015	0.89	336.19	−15.67	−0.006	−0.98	−0.13
15 March, 2015	0.85	330.99	−5.02	−0.005	−1.13	−0.06
15 April, 2015	0.87	338.49	−9.03	−0.005	−0.58	−0.25
15 May, 2015	0.83	350.82	−3.16	−0.007	−0.66	−0.32
15 June, 2015	0.86	300.98	−0.42	−0.005	−0.30	0.03
15 July, 2015	0.90	326.60	−2.70	−0.006	−0.23	−0.21
15 August, 2015	0.85	308.03	−2.69	−0.005	−0.55	0.07
15 September, 2015	0.87	295.45	−2.60	−0.004	−0.34	0.07
15 October, 2015	0.87	310.69	−2.56	−0.005	−0.57	−0.03
15 November, 2015	0.85	331.74	−1.86	−0.005	−1.25	−0.06
15 December, 2015	0.79	324.41	−12.13	−0.005	−0.98	−0.09

. 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

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

Date	Model Name	Nugget	Partial Sill	Range (km)
15 January, 2015	exponential	0.34	0.66	4.03
15 February, 2015	exponential	0.67	0.32	11.59
15 March, 2015	exponential	0.38	0.62	5.62
15 April, 2015	spherical	0.50	0.48	5.29
15 May, 2015	exponential	0.40	0.60	7.62
15 June, 2015	exponential	0.57	0.43	9.25
15 July, 2015	spherical	0.69	0.31	10.13
15 August, 2015	gaussian	0.41	0.57	9.78
15 September, 2015	spherical	0.44	0.55	14.39
15 October, 2015	exponential	0.32	0.68	13.78
15 November, 2015	spherical	0.52	0.46	13.28
15 December, 2015	exponential	0.68	0.32	18.56

and

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

Date	Model Name	Nugget	Partial Sill	Range (km)
15 January, 2015	exponential	0.34	0.66	6.34
15 February, 2015	spherical	0.35	0.65	12.08
15 March, 2015	spherical	0.57	0.43	10.17
15 April, 2015	spherical	0.56	0.43	9.57
15 May, 2015	spherical	0.39	0.61	17.30
15 June, 2015	exponential	0.45	0.55	3.56
15 July, 2015	exponential	0.53	0.47	10.54
15 August, 2015	gaussian	0.31	0.67	10.75
15 September, 2015	exponential	0.39	0.61	9.57
15 October, 2015	spherical	0.29	0.71	14.79
15 November, 2015	exponential	0.50	0.50	12.71
15 December, 2015	exponential	0.30	0.70	9.28

.

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.

2.3.1. BME Epistemic Paradigm and Conceptual Core

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.

2.3.2. BME Framework

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:

$$f_G\left({\mathsf{\chi }}_k,{\mathsf{\chi }}_{data}\right)=f_G\left({\mathsf{\chi }}_{map}\right)=\frac{\exp \left({\sum }_{\alpha =1}^N{\mu }_{\alpha }\left(p_{map}\right)g_{\alpha }\left({\mathsf{\chi }}_{map}\right)\right)}{\int d{\mathsf{\chi }}_{map}\exp \left({\sum }_{\alpha =1}^N{\mu }_{\alpha }\left(p_{map}\right)g_{\alpha }\left({\mathsf{\chi }}_{map}\right)\right)},$$

(1)

where ${\mathsf{\chi }}_k$ denotes the LST of the estimated pixel, ${\mathsf{\chi }}_{data}$ = (${\mathsf{\chi }}_{hard}$, ${\mathsf{\chi }}_{soft}$), ${\mathsf{\chi }}_{hard}$ represents hard data, and ${\mathsf{\chi }}_{soft}$ 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. $f_G$(${\mathsf{\chi }}_k$, ${\mathsf{\chi }}_{data}$) denotes the prior pdf of the map ${\mathsf{\chi }}_{map}$ = (${\mathsf{\chi }}_k$, ${\mathsf{\chi }}_{data}$) given the general knowledge base G. ${\mu }_{\alpha }$($p_{map}$) represents Lagrange multipliers. $g_{\alpha }$(${\mathsf{\chi }}_{map}$) is a set of known functions of ${\mathsf{\chi }}_{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_K\left({\mathsf{\chi }}_k\right)=f_G\left({\mathsf{\chi }}_k{|\mathsf{\chi }}_{hard},{\mathsf{\chi }}_{soft}\right)=f_G\left({\mathsf{\chi }}_k{|\mathsf{\chi }}_{data}\right)=f_G\left({\mathsf{\chi }}_k,{\mathsf{\chi }}_{data}\right)/f_G\left({\mathsf{\chi }}_{data}\right)$$

(2)

In Equation (2), $f_K$ denotes the posterior pdf of the map ${\mathsf{\chi }}_k$, 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.

$${\mathsf{\chi }}_k, \mathrm{mean}=\int {\mathsf{\chi }}_k{f}_K\left({\mathsf{\chi }}_k\right)d{\mathsf{\chi }}_k$$

(3)

We used the BME mean value (${\mathsf{\chi }}_k, \mathrm{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 (Figure 2 and Figure 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. 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.

Date		Observed	Reconstructed		Observed	Reconstructed
15 January, 2015	Daytime	39.8%	100%	Nighttime	54.6%	100%
15 February, 2015		40.6%	100%		53.2%	100%
15 March, 2015		40.1%	100%		45.2%	100%
15 April, 2015		52.8%	100%		58.2%	100%
15 May, 2015		39.6%	100%		46.2%	100%
15 June, 2015		33.6%	100%		30.1%	100%
15 July, 2015		39.1%	100%		48.1%	100%
15 August, 2015		36.7%	100%		42.2%	100%
15 September, 2015		48.6%	100%		56.8%	100%
15 October, 2015		63.0%	100%		74.8%	100%
15 November, 2015		46.2%	100%		48.0%	100%
15 December, 2015		34.8%	100%		52.7%	100%
Average		42.9%	100%		50.9%	100%

).

In general, the entire study area showed strong spatial heterogeneity that varied in different seasons for both the day and night (Figure 2 and Figure 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 (Figure 2 and Figure 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.

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 (Figure 2 and Figure 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.

3.2. Accuracy Assessment

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 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.

The daytime accuracy was mostly lower than that of the nighttime (Figure 5 and Figure 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 (Figure 5, Figure 6 and Figure 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 (Figure 5 and Figure 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 Figure 5 and Figure 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.

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² 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 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.

3.3. Factors that Influence Accuracy

Figure 9 and Figure 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², 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² of multiple linear regression did not affect the accuracy of the method. In this study, R² 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.

3.4. Comparisons with Other Methods

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.

4. Discussion

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² in this study was nearly 0.6, and in the subsequent application of the BME method to LST reconstruction, when the average regression R² 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 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.

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. Omnidirectional and directional covariance model parameters for 18 October, 2015.

Time (omnidirectional)	Model name	Nugget	Partial Sill	Range (km)
Daytime	spherical	0.20	0.76	12.94
Nighttime	spherical	0.20	0.79	15.59
Time (directional)	The angle between the principal axis and the horizontal axis (°)		Principal/secondary axes (km/km)
Daytime	0		20.68/7.69
Nighttime	3.52		18.39/10.06

, 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.

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. Omnidirectional and directional accuracy verification results for 18 October, 2015.

Time	Model Name	MAE (K)	RMSE (K)	R2
Daytime	omnidirectional	0.79	1.32	0.98
	directional	0.78	1.30	0.98
Nighttime	omnidirectional	0.45	0.71	0.99
	directional	0.45	0.71	0.99

). It is suggested that directionality can be considered during the daytime in the following research.

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² 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² 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.