Next Article in Journal
GPS/BDS Medium/Long-Range RTK Constrained with Tropospheric Delay Parameters from NWP Model
Previous Article in Journal
Tree Species Classification of the UNESCO Man and the Biosphere Karkonoski National Park (Poland) Using Artificial Neural Networks and APEX Hyperspectral Images
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Reconstruction of MODIS Land Surface Temperature Products Based on Multi-Temporal Information

1
Key Laboratory of Remote Sensing of Gansu Province, Heihe Remote Sensing Experimental Research Station, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China
2
CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100049, China
3
Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100101, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(7), 1112; https://doi.org/10.3390/rs10071112
Submission received: 5 June 2018 / Revised: 4 July 2018 / Accepted: 10 July 2018 / Published: 12 July 2018

Abstract

:
Land surface temperature (LST) products derived from the moderate resolution imaging spectroradiometer (MODIS) sensor are one of the most important data sources used to research land surface energy and water balance at regional and global scales. However, MODIS data are severely contaminated by cloud cover, which limits the applications of LST products. In this paper, based on the spatio-temporal autocorrelation of land surface variables, a reconstruction algorithm depending on the correlations between spatial pixels in multiple time phases from available MODIS LST data is developed to reconstruct clear-sky LST values for missing pixels. Considering the impacts of correlation and bias between predictors and reconstructed data on the modeling error, the known data in the reconstructed time phase are combined with the data temporally nearest to them as predictor variables to establish their temporal relationships with the reconstructed data. The reconstructed results are validated by a series of evaluation indices. The average correlation coefficient between the reconstructed results and ground-based observations is 0.87, showing high temporal change accuracy. The difference in Moran’s I, representing spatial structure characteristics between the known and reconstructed data, is 0.03 on average, indicating a slight loss of spatial accuracy. The average reconstruction rate is approximately 87.0%. The modeling error, as part of the reconstruction error, is only 1.40 K on average and accounts for 5.0% of the total error. If the product and modeling errors are removed, the residual error represents approximately 3.5 K and 5.6 K of the annual mean difference between the cloudy and cloudless LST at night and during the day, respectively. In addition, different reconstruction cases are demonstrated using various predictor data, including many combinations of multi-temporal MODIS LST data, the microwave brightness temperature, and the combination of the normalized difference vegetation index and terrain data. Comparisons among cases show that the known MODIS LST data are more reliable as predictor variables and that the data combination advocated in this paper is optimal.

Graphical Abstract

1. Introduction

Land surface temperature (LST) plays an important role in land-atmosphere interactions. Therefore, LST data are widely required for studies on climate change, the water cycle, and energy budgets [1,2,3]. With the development of remote sensing technology, the LST data obtained by transforming satellite-based thermal infrared observations are more valuable [4] and can better meet regional and global research requirements than ground-based observations. Due to high spatio-temporal resolution, the moderate resolution imaging spectroradiometer (MODIS) sensor has become one of the most widely used sources of satellite-derived LST data [5]. However, daily MODIS LST products are strongly influenced by clouds, and only clear-sky pixel values are available, which leads to limited research applications. To overcome this restriction, several MODIS LST reconstruction algorithms have been developed.
To date, most of studies on the reconstruction of MODIS LST products have been focused on recovering clear-sky LST values for missing pixels. There are three categories of statistical methods using different information, including spatial, temporal, and spatio-temporal information. Spatial methods are the most common, wherein the relationship in space between the auxiliary data with spatial continuity and the available LST data is adopted, and the missing LST data are filled in by sharing the relationship model assuming that the known and unknown LST data have the same statistical relationship with the auxiliary data [6]. Neteler [7] used the digital elevation model (DEM) as a covariate to generate reconstructed LST maps through linear regression when the LST map consisted of over 10% valid pixels; otherwise, 16-day MODIS LST data were used as the target variable. Fan, et al. [8] introduced land cover, normalized difference vegetation index, and MODIS band 7 as auxiliary data to fill the missing pixels using three methods: linear regression, regression tree analysis, and artificial neural networks. Ke, et al. [9] adopted latitude, longitude, DEM, and normalized difference vegetation index (NDVI) as predictor variables to establish a regression model used to capture large-scale spatial variability; they then applied the Kriging method to generate small-scale spatial information to supplement the regression prediction [10]. For this category of methods, modeling accuracies are determined by the correlation between LST and its covariates. However, the LST variable is synthetically affected by many environmental factors, and limited kinds of predictor data cannot accurately represent the spatial heterogeneity of LST. The second category of methods is based on temporal information. Xu and Shen [11] used the harmonic analysis of time series (HANTS) algorithm [12] to fit available time series of LST data at each pixel and then recovered missing pixels using fitted values of the temporal model. Na, et al. [13] used another temporal method, called Savitzky-Golay filter [14], to reconstruct MODIS LST data under clear-sky conditions. However, the temporal filter method smooths LST diurnal variation and can represent only seasonal LST changes. Thus, a temporal method such as HANTS or Savitzky-Golay is not generally suitable to reconstruct daily MODIS LST products. The third category consists of spatio-temporal methods performed by establishing a regression relationship between similar pixels and missing pixels. Yu, et al. [15] selected pixels whose environmental factor vector (composed of DEM, slope, aspect, NDVI, and solar radiation) had the nearest Euclidean distance to the reconstructed pixel as the most similar pixel. However, the consistency in direction between environmental factor vectors of available and reconstructed pixels was ignored. It is possible that both vectors have almost the same magnitude but opposite directions, indicating a low degree of similarity between two pixels. Zeng, et al. [16] classified multi-temporal LST data and then considered several pixels in the class which the reconstructed pixels belong to as similar pixels. The more clearly land surface types are distinguished, the more accurate the regression relationship; however, the use of too many classes can lead to a lack of available similar pixels. Sun, et al. [17] defined similar pixels according to spatial distance. Pixels closer to the reconstructed pixel made a greater contribution to the reconstruction process.
The above methods can reconstruct cloud-free LST values but not cloudy ones; they are nonetheless valuable for long-term trend analysis and further estimation of real LST values under clouds. Few studies on recovering pixels under cloudy conditions have been performed based on model-based and statistical methods. Due to the difficultly of acquiring additional information, e.g., radiation and flux data at large spatial scales, model-based reconstruction methods have been given little attention. Jin [18] proposed a neighboring pixel method based on the surface energy balance theory assuming that the difference in LST values between a cloudy pixel and its neighboring clear-sky pixel is caused by differences in energy fluxes, e.g., net radiation and sensible and latent heat. Lu, et al. [19] improved the neighboring pixel method by introducing satellite observations. Zhang, et al. [20] used a one-dimensional heat transfer model to reconstruct cloudy pixels. In statistical methods, the solution for recovering cloudy MODIS LST pixels is to introduce microwave temperature brightness (TB) data due to penetrating clouds. Based on the strong correlation between infrared LST and TB data, TB data, particularly in the high-frequency range, have been widely used to estimate land surface emissivity [21,22] and LST [23,24] by employing infrared LST as auxiliary data. However, due to the difference in spatial scale, the 25–40 km microwave pixels cannot be directly filled into gaps of MODIS LST data. Kou, et al. [25] reconstructed missing MODIS LST pixels under cloudy conditions by merging microwave data using a Bayesian Maximum Entropy method. However, a complex land surface environment, e.g., with variations in terrain and vegetation, may reduce the correlation between LST and TB and thereby increase the reconstruction error.
In this paper, a reconstruction algorithm based on spatio-temporal information is used to recover daily MODIS LST data. The available LST data at different spatial pixels and times are used as predictor variables and are optimally selected by analyzing the correlation and bias between multi-temporal MODIS LST data. In addition, reconstructed results employing different kinds of predictor data, e.g., NDVI, DEM, and TB, are compared with each other.

2. Materials

2.1. Study Area

This work was performed in the Babao River Basin, which is situated in the eastern branch of the upper reach of the Heihe River on the northeastern margin of Tibetan Plateau in Northwestern China (Figure 1). The Babao River Basin is a cloudy area, and the cloudy time generally exceeds half of a year due to a combination of westerlies, East Asia monsoon, and Tibetan Plateau monsoon. The stronger the atmospheric circulation becomes, the more there is cloud cover. In addition, the characteristics of the study area increase the probability of clouds. The elevation of the Babao River Basin is 3604 m on average and ranges from 2640 m to 5000 m; the annual rainfall is approximately 400 mm, which drives the formation of orographic clouds by forcing humid air to rise. The high vegetation cover and soil water content in the study area increase the surface evapotranspiration under the strong solar radiation, especially in summer afternoons, and the enhanced convective motion of the air promotes the formation of convective clouds, causing the cloud cover to be higher during the day than the night. The Babao River Basin is an ideal study area to research the reconstruction of MODIS LST products.

2.2. Data

MODIS Aqua LST products are taken as the experimental data to be reconstructed. Figure 2 shows the spatial and temporal distribution of pixels requiring reconstruction. The missing Aqua LST data at night are substantially fewer and more homogeneous in space than those during the day; the highest intensity appears from April to October. The Aqua LST data are severely affected by clouds during the day, showing a high missing data rate throughout the year. There is a significant correlation between the number of absent data and elevation, indicating that a higher elevation leads to a larger number of missing data.
A large amount of remote sensing data is employed to perform the reconstruction algorithm (see Table 1). All available pixel values derived from MODIS Aqua and Terra LST products (10.5067/MODIS/MOD11A1.006 and 10.5067/MODIS/MYD11A1.006) are potential predictor data for filling missing pixels. To compare reconstruction results based on different data sources, two groups of predictor variables are used to reconstruct MODIS Aqua LST data. One group is multi-frequency and multi-polarization TB data derived from the Advanced Microwave Scanning Radiometer 2 (AMSR2) sensor (http://gcom-w1.jaxa.jp), and the other group is NDVI (10.5067/MODIS/MYD13A2.006) and terrain data, e.g., slope and aspect calculated using DEM data derived from shuttle radar topography mission (SRTM, https://lta.cr.usgs.gov/SRTM). As shown in Table 2, the upwelling and downwelling longwave radiation data observed by six automatic meteorological stations (10.3972/hiwater.248.2015.db, 10.3972/hiwater.249.2015.db, 10.3972/hiwater.255.2015.db, 10.3972/hiwater.254.2015.db, 10.3972/hiwater.253.2015.db, and 10.3972/hiwater.252.2015.db) are applied to estimate the ground-based LST through combination with MODIS emissivity observations (10.5067/MODIS/MYD11B1.006). The calculated ground-based LST data are employed to validate the reconstruction results. All ground-based and remotely sensed observations from January 2014 to December 2014 are adopted. NDVI = normalized difference vegetation index; AMSR2 = Advanced Microwave Scanning Radiometer 2; TB = temperature brightness.

2.3. Data Reprocessing

2.3.1. Resampling

The spatial reconstruction method employing terrain data as predictor data requires a consistent spatial scale between the predictor and target variables. First, the slope and aspect values are calculated using an algorithm incorporating the DEM values of the calculated cell’s eight neighbors [26], and terrain data, including DEM, slope, and aspect, are then resampled at the same spatial scale as the MODIS LST data by averaging all 90 m pixels in the 1 km grids.

2.3.2. Processing Outliers

Outlier data may increase prediction uncertainties. Generally, it is recommended to leave the data in the 95% confidence interval, which cannot directly be applied to temporal data due to the existence of trends. It is feasible to find time series outlier values by identifying the unreasonable bias between the data and their trend. The criterion to select available data that are not affected by clouds is defined as follows:
μ 1.96 σ ε b i a s μ + 1.96 σ  
where ε b i a s is calculated by removing the temporal trend from data and μ and σ are the mean and standard deviation of ε b i a s , respectively. The data with eligible ε b i a s are retained.

2.3.3. LST Retrieval from Ground-Based Measurements

Ground-based radiation observations cannot be used to evaluate reconstructed results and need to be converted into LST values. Based on thermal radiative transfer theory, the ground-based LST T i n s i t u can be calculated using in-situ longwave radiation observations [27]:
T i n s i t u = [ L i n s i t u ( 1 ε b ) L i n s i t u ε b δ ] 1 4  
where L i n s i t u and L i n s i t u indicate ground-based upwelling and downwelling longwave radiation, respectively. δ is the Stefan-Boltzmann constant ( 5.67 × 10 8   W m 2 K 4 ) and ε b is the broadband emissivity, which can be estimated by MODIS narrowband emissivity products [28,29]:
ε b = 0.2122 ε 29 + 0.3859 ε 31 + 0.4029 ε 32  
where ε 29 , ε 31 , and ε 32 are the narrow emissivities of MODIS bands 29, 31, and 32, respectively.

3. Methodology

3.1. MODIS LST Reconstruction Algorithm

In consideration of the spatio-temporal autocorrelation characteristics of land surface variables, the MODIS LST data that are not affected by clouds as the predictor variables are more reliable than other data sources, e.g., TB and NDVI, for predicting missing MODIS LST pixels. MODIS LST data are observed at 01:30, 10:30, 13:30, and 22:30 local time. The missing MODIS LST pixel values can be estimated based on the spatio-temporal correlation among LST values at different spatial pixels and times. Due to the strong temporal variability, suppose that the temporal correlation range is less than one day. Then, the MODIS LST pixels without values at the jth pixel at local time t on the dth day are filled using the equation as follows:
T j d , t = i = 1 n ( d , t 1 ) w i t 1 T i d , t 1 + i = 1 n ( d , t 2 ) w i t 2 T i d , t 2 + i = 1 n ( d , t 3 ) w i t 3 T i d , t 3 + i = 1 n ( d , t 4 ) w i t 4 T i d , t 4  
where T i d , t k is the available MODIS LST value situated at the ith pixel and observed at local time t k on the dth day. In eight directions with equal intervals from 0° to 360°, the nearest T i d , t k to T j d , t in the spatial distance measurement is selected to maintain the spatial continuity between the reconstructed and known LST data. t k ( k = 1 , 2 , 3 , 4 ) represents four MODIS LST data observation times. n ( d , t k ) indicates the number of data T i d , t k . w i t k , the weighting coefficient, which is assigned to T i d , t k , can be estimated using the following equation:
T n ( d , t k ) m , t k ω n ( d , t k ) t k = T j m , t  
where ω n ( d , t k ) t k is the weighting coefficient vector and each element is matched with w i t k in Equation (4). T j m , t is a temporal vector, including the known LST data at the jth pixel at time t on the mth ( m = 1 , 2 , 3 , ; d m ) day, represented as [ T j 1 , t , T j 2 , t , , T j m , t ] . The observation matrix T n ( d , t k ) m , t k is constructed using LST time series data [ T n ( d , t k ) 1 , t , T n ( d , t k ) 2 , t , , T n ( d , t k ) m , t k ] and is written as follows:
[ T 1 1 , t 1 T n ( d , t 1 ) 1 , t 1 T 1 1 , t 2 T n ( d , t 2 ) 1 , t 2 T 1 1 , t 3 T n ( d , t 3 ) 1 , t 3 T 1 1 , t 4 T n ( d , t 4 ) 1 , t 4 T 1 2 , t 1 T n ( d , t 1 ) 2 , t 1 T 1 2 , t 2 T n ( d , t 2 ) 2 , t 2 T 1 2 , t 3 T n ( d , t 3 ) 2 , t 3 T 1 2 , t 4 T n ( d , t 4 ) 2 , t 4                                                    T 1 m , t 1 T n ( d , t 1 ) m , t 1 T 1 m , t 2 T n ( d , t 2 ) m , t 2 T 1 m , t 3 T n ( d , t 3 ) m , t 3 T 1 m , t 4 T n ( d , t 4 ) m , t 4 ]  
However, if the correlation among column vectors in the observation matrix is high, especially when both sets of sequence data are spatially close, invalid weighting coefficients that lead to over-fitting will be estimated. This problem can be addressed by performing ridge regression [30,31]:
ω ^ n ( d , t k ) t k = [ ( T n ( d , t k ) m , t k ) T n ( d , t k ) m , t k + λ I ] 1 ( T n ( d , t k ) m , t k ) T j m , t  
where ω ^ n ( d , t k ) t k is the estimator of the weighting coefficient vector. λ is a regularization parameter. Herein, we set λ to 0.1.

3.2. Evaluation Index

The performance of the reconstruction algorithm is considered from four aspects: (1) the level of uncertainty in the reconstruction process; (2) whether the reconstructed data is consistent with the temporal changes with the known data; (3) whether the reconstructed data maintains similar spatial structure characteristics to the known data; and (4) the value of the reconstructed rate.
(1)
Error Analysis:
When the ground-based LST data are used as the ground truth to evaluate reconstructed data, the sources of reconstruction error can be decomposed as follows:
ε r e = ε p r o d u c t + ε c l o u d y s k y + ε m o d e l i n g  
where the reconstruction error ε r e can be estimated by calculating the differences between the reconstructed data and the ground-based LST data:
ε r e = E ( T r e T i n - s i t u ) 2  
where E ( · ) represents the mean and T r e and T i n - s i t u are the reconstructed value and ground-based validation data calculated by Equation (2), respectively.
ε p r o d u c t indicates the accuracy of the MODIS LST product, which can be obtained by evaluating the known LST data:
ε p r o d u c t = E ( T k n o w n T i n - s i t u ) 2  
where T k n o w n is the LST data that are not affected by clouds.
ε m o d e l i n g , the modeling error representing the modeling reliability, can be obtained from the residual errors of the model described by Equation (5):
ε m o d e l = E ( T j m , t T n ( d , t k ) m , t k ω ^ n ( d , t k ) t k ) 2  
The missing LST data are reconstructed in clear sky days, but the ground-based LST values under clouds are employed as validation data. Thus, the difference ε c l o u d y s k y between the reconstructed data T r e and true LST data under the clouds is also an important component of ε r e . Suppose that all errors ε are uncorrelated and have a zero mean and that ε c l o u d y s k y can be written as
ε c l o u d y s k y = D ( ε r e ) D ( ε p r o d u c t ) D ( ε m o d e l i n g )  
where D ( · ) represents variance.
(2)
Temporal Consistency:
The stochastic characteristics of reconstruction errors can decrease the temporal consistency between the reconstructed data and the ground-based LST data. The correlation coefficient is typically used to evaluate the prediction accuracy of temporal changes:
r r e = C o v ( T r e , T i n - s i t u ) D ( T r e ) D ( T i n - s i t u )  
where C o v ( · ) and D ( · ) indicate covariance and variance, respectively.
(3)
Spatial Structure:
In addition to temporal variability, spatial heterogeneity is a significant inherent characteristic of land surface variables. Thus, it is necessary to evaluate the spatial structure of reconstructed LST data. Moran’s I is generally used to quantifiably evaluate the spatial autocorrelation of the variable:
I = n i = 1 n j = 1 n w i j i = 1 n j = 1 n w i j ( T i T ¯ ) ( T j T ¯ ) i = 1 n ( T i T ¯ ) 2  
where T i is the attribution value for feature i . w i j represents the spatial weight between features i and j and is estimated using the inverse distance between T i and T j . n is equal to the total number of features.
(4)
Reconstruction Rate
The reconstruction rate is an index to evaluate algorithm performance. The reconstruction rate is defined as follows:
R R = N r e N m i s s i n g × 100 %  
where N m i s s i n g is the number of the missing pixels and N r e is the number of reconstructed pixels.

4. Results and Discussion

4.1. Selection of Predictor Data

If the correlation coefficient between both temporal variables is equal to 1, when using one as a predictor variable to estimate the other, the prediction error will be 0. However, a complete correlation is almost impossible, and then temporal bias between the predictor and target variables affects prediction accuracy. Due to the use of multi-temporal LST data as predictor data, bias between both sets of temporal data observed at different local times is obvious. Figure 3 shows the influence of the combination of correlation coefficient and absolute bias on the prediction accuracy, which is recorded using a temporal MODIS LST observation at one pixel to estimate those at other pixels via Equation (5). With the decrease of the correlation coefficient and the growth of bias between LST data, the prediction errors gradually increase. Thus, it is necessary to optimally select multi-temporal LST data as predictor variables to decrease prediction uncertainty.
In MODIS LST data reconstruction, Aqua LST products are taken as the experiment data to perform the reconstruction algorithm. Based on the above analysis, the suitable predictor data need to be selected for the reconstruction of Aqua LST products. Figure 4 and Figure 5 show analyses of correlation and bias, respectively, between the LST data at the reconstructed time and themselves and the LST data at other times. The highest average correlation and smallest average bias occur between the Aqua LST data and themselves, showing an average correlation coefficient of 0.98 at night and 0.86 during the day and an average bias of 2.2 K at night and 6.0 K during the day. Thus, the Aqua LST data that are not affected by clouds are the best predictor data. According to Figure 3, if the available Aqua data are used to predict the missing pixels, the average prediction errors are 1.6 K at night and 3.9 K during the day.
However, the reconstruction rate is low if only using the available Aqua LST data at the reconstructed time. Without substantially increasing the prediction errors, the multi-temporal LST data can also be taken as predictor data. For the nighttime Aqua data, the correlation and difference with nighttime Terra data show the following average correlation coefficient and bias, representing 0.96 and 2.7 K, respectively, which leads to an average prediction error of 2.3 K. For daytime Aqua data, the average correlation coefficients with MODIS LST data at three other times are approximately 0.82, but the daytime Terra data have the smaller bias of 6.7 K, and the average prediction error is approximately 4.3 K.
Through above analysis, in addition to the Aqua data that are not contaminated by clouds, the nighttime Terra data one day ahead may be introduced as predictor data to reconstruct nighttime Aqua data, and daytime Terra data may be used to reconstruct daytime Aqua data.

4.2. Reconstruction of MODIS Aqua LST Data

Figure 6 shows the monthly reconstruction rates of MODIS Aqua LST data. The average percentage of reconstructed pixels is 92.8% at night and 81.7% during the day. The average reconstruction rates in spring and winter are higher than those in summer and autumn. Depending on the reconstruction algorithm, the missing MODIS LST data can be completely filled as long as the predictor data exist. Thus, a high missing rate does not imply a low reconstruction rate (e.g., missing pixels accounting for more than 90% of daytime Aqua LST data in December are recovered to a degree of 96.7%). Even if MODIS LST data are completely affected by clouds, the absence of pixel values is still estimated by introducing other available temporal LST data. MODIS Aqua LST products have a number of data with 100% loss, specifically, 67 at night and 121 during the day, and 42 nighttime and 66 daytime empty data are reconstructed, indicating that multi-temporal data can be used to increase the reconstruction rate.

4.3. Spatial Characteristics of Reconstructed Data

The examples of the reconstructed Aqua LST data with different missing rates in space are shown in Figure 7. Strong spatial continuity generally exists between the known and reconstructed data, and the renewed nighttime data show better continuity. Near the northeast boundary at coordinate (100.4°, 38.2°) and around coordinate (100.2°, 38.1°), the reconstructed daytime data exhibit spatial randomness caused by the stronger temporal heterogeneity of the daytime data (see Figure 5) and severe data loss (see Figure 2) representing a lack of prior knowledge of the target variable T j m , t in Equation (5) at those two locations.
The LST variable exhibits spatial autocorrelation [32], which can be quantified by Moran’s I index. The closer the value of Moran’s I is to 1, the stronger the spatial autocorrelation. To detect the ability of the reconstruction algorithm to recover spatial structure, the differences in the Moran’s I index of monthly data between the reconstructed and known data were calculated, as shown in Figure 8. The average absolute differences are 0.02 at night and 0.04 during the day, meaning that the spatial characteristics of reconstructed nighttime data are closer to the known data. The average spatial autocorrelations of the known nighttime and daytime data are separately quantified as 0.93 and 0.86, respectively, using Moran’s I index, indicating that the LST variable has stronger randomness during the day than at night (especially the daytime data in June, July, and August, with an average Moran’s I value of 0.77). Spatial randomness increases the difficulty of recovering spatial characteristics. In addition to the weak spatial autocorrelation, severe data loss is a significant factor resulting in missing data that cannot be accurately reconstructed in terms of spatial structure because it is challenging to capture spatial randomness when employing few data.

4.4. Validation Using Ground-Based Observations

The upwelling and downwelling longwave radiation observations at six ground-based stations (see Table 2) are transformed into LST values using Equation (2) to validate the reconstructed results. The comparisons between ground-based observations and reconstructed Aqua LST data are shown in Figure 9 and Figure 10, for the nighttime and daytime, respectively. The missing data affected by clouds are estimated under clear sky conditions, but the ground-based observations used to validate the reconstructed data represent values under clouds. Thus, the fitted line of the reconstructed nighttime data (red dotted line) is located below the fitted line of known data (blue dotted line) due to the higher LST under clouds than under clear sky caused by the downwelling longwave radiation of clouds at night. In contrast, the reconstructed daytime data values are higher than the ground-based LST data because the solar radiation reaching the ground is reduced by clouds, as indicated by the fitted line of the reconstructed data lying above the fitted line of the known data.
The detailed accuracy information of reconstructed data is shown in Table 3. The goal of the algorithm is to reconstruct cloud-free LST data; thus, the reconstruction error ε r e calculated using cloudy ground-based LST data cannot exactly represent the achievement of the goal. As shown in Equation (8), ε r e comprises three parts: ε p r o d u c t , ε c l o u d y s k y , and ε m o d e l i n g . The unchangeable product error ε p r o d u c t makes a significant contribution to ε r e , and the average percentage of ε p r o d u c t to ε r e is 54.3% at night and 44.2% during the day. The modeling error ε m o d e l i n g can be taken as an indirect measure of the accuracy of reconstructed LST data under clear sky. The smaller ε m o d e l i n g is, the better ε c l o u d y s k y expresses the differences in the LST values between cloudy and cloud-free conditions. In this case, the average ε m o d e l i n g is approximately 1.4 K and accounts for only 5.0% and 4.0% of ε r e at night and during the day, respectively. Thus, ε c l o u d y s k y can indicate the average temporal differences between cloudy and cloud-free LST values, which is reliable information to further recover LST values under clouds. However, it is more challenging to recover daytime cloudy LST values due to the stronger spatial heterogeneity of ε c l o u d y s k y during daytime, with considerable differences existing among the six stations.
In addition, the correlation coefficient r r e calculated using Equation (13) is adapted to measure the temporal consistency between ground-based observations and reconstructed data. r r e is close to the reference value r k n o w n representing the correlation coefficient between ground-based observations and the known Aqua LST data. The average difference between r k n o w n and r r e is 0.05 at night and 0.03 during the day, indicating that the reconstructed data can accurately maintain temporal consistency with the ground validation data. In addition to the reconstruction accuracy, the reconstruction rate is an important index of the performance of the reconstruction algorithm. The average missing percentages of time series at the pixels where the six stations are located are 47.7% at night and 76.4% during the day, which corresponds to Figure 2. The average reconstruction rates of MODIS Aqua at these six pixels exceed 73.0%, and the average reconstruction rate of nighttime data exceeds 90%.

4.5. Comparisons among Reconstruction Results Based on Different Predictor Data

To demonstrate the influence of different predictor data on reconstruction accuracy, five cases are designed using different data sources as predictor variables: the available MODIS LST data at the reconstructed time (case A), the optimal combination of multi-temporal MODIS LST data analyzed in Section 4.1 (case B), the combination of all temporal MODIS LST data (case C), TB data (case D), and the combination of NDVI and Terrain data (case E). Each case is performed by replacing temporal data vectors in matrix T n ( d , t k ) m , t k of Equation (6) using its own predictor data except in the NDVI-terrain-based case, which is implemented by establishing the matrix composed of spatial data vectors at the reconstructed time. The evaluation indices of the five cases are shown in Figure 11. For case D, adopting TB data, only the accuracies of the temporal changes and spatial structure of the reconstructed daytime data are acceptable. Other indices are not ideal because the high correlation proven by previous research between LST and TB is not shown in this study area, which may be caused by the complex terrain. For case E, the correlation coefficient is high depending on the temporal consistency between NDVI and LST, but the combination of NDVI and terrain data cannot accurately represent the spatial heterogeneity of LST leading to a large modeling error and Moran’s I values far from the reference values. In addition, the reconstruction rates in case E are also unsatisfactory. For cases employing MODIS LST data as predictor variables, the differences in the correlation coefficient and spatial structure among reconstruction results are slight. Case A, adopting the known MODIS LST data at the reconstructed time as the predictor data, has the lowest modeling errors due to this case having the highest temporal correlation and lowest temporal bias between the predictor and reconstructed data (see Figure 4 and Figure 5), but it has low reconstruction rates. To improve the reconstruction rate, cases B and C introduce multi-temporal LST data. For case B, in addition to the known data at the reconstructed time, the predictor data at the nearest time to reconstructed data are adopted without seriously increasing the modeling error, and the accuracies on the correlation coefficient and spatial structure are consistent with those of case A. Case C employs all temporal LST data, increasing the reconstruction rate. However, relative to cases A and B, the reconstruction accuracies decrease, especially for the reconstructed daytime data. Overall, case B advocated in this paper provides a favorable balance between reconstruction accuracy and reconstruction rate.

5. Conclusions

MODIS LST products are one of the most important remote sensing LST observations. However, due to clouds, remotely sensed LST data contain a large number of missing values. It is very important to reconstruct MODIS LST data to perform large-scale and long-term LST application studies.
In this paper, the reconstruction goal is to obtain cloud-free LST values, and the recovery accuracy is decided by the correlation between predictor and reconstructed data. Relative to other predictor data, e.g., BT, NDVI and terrain data, the available MODIS LST data that are not affected by clouds have a higher correlation with the reconstructed data due to the autocorrelation characteristics of the land surface variable. Thus, using the available MODIS LST data as predictor data is the most reliable approach.
Considering the high correlation between time series MODIS LST data at different spatial pixels and times, a spatio-temporal algorithm employing multi-temporal MODIS LST data is adopted to reconstruct MODIS LST data. When the available LST data at the reconstructed time are taken as the predictor data, the algorithm shows satisfactory performance except for the reconstruction rate. When the extra LST data that are nearest in time to the reconstructed data are introduced, the reconstruction rate increases without a severe decrease in reconstruction accuracy; this approach is thus advocated as the optimal case. Although the reconstruction rate can be further improved by adopting all temporal LST data, the overall reconstruction accuracies significantly decrease due to the low correlation and high bias of the time series between nighttime and daytime data.
Although the real LST under clouds cannot be recovered based on the reconstruction algorithm, the differences between the cloudy and cloud-free LST can be estimated by the error decomposition as long as the modeling accuracy is high enough, which can provide important evidence to further estimate ground truth values. However, the differences between the LST under clouds and clear sky show stronger spatial heterogeneity during the day than at night, which increases the challenge of estimating real LST data under clouds.
BT data have the potential to obtain cloudy LST values due to microwaves penetrating through cloud cover and the strong correlation between LST and BT. In this paper, however, the high correlation does not appear, which may be caused by the complex terrain in the study area. In future work, we plan to use the reconstruction algorithm to obtain cloudy LST data by merging multiple sources of remotely sensed microwave data in flat areas.

Author Contributions

J.K. and R.J. designed the work, implemented the concept, and analyzed the results. J.K. wrote the manuscript. J.T., X.L., and Y.Z. provided suggestions.

Funding

This work was jointly supported by the Strategic Priority Research Program of the Chinese Academy of Science (Grant No. XDA19070104), the National Natural Science Foundation of China (Grant Nos. 41531174, 41701419 and 41701420), and the Foundation for Excellent Youth Scholars of NIEER, CAS.

Acknowledgments

MODIS data were retrieved from the online Data Pool, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota (https://e4ftl01.cr.usgs.gov). TB data was provided by JAXA (http://gcom-w1.jaxa.jp). DEM data is provided by USGS (https://lta.cr.usgs.gov/SRTM).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Corbari, C.; Mancini, M. Calibration and validation of a distributed energy–water balance model using satellite data of land surface temperature and ground discharge measurements. J. Hydrometeorol. 2014, 15, 376–392. [Google Scholar] [CrossRef]
  2. Dai, A.; Fyfe, J.C.; Xie, S.-P.; Dai, X. Decadal modulation of global surface temperature by internal climate variability. Nat. Clim. Chang. 2015, 5, 555–559. [Google Scholar] [CrossRef]
  3. Sun, Z.; Wang, Q.; Batkhishig, O.; Ouyang, Z. Relationship between evapotranspiration and land surface temperature under energy-and water-limited conditions in dry and cold climates. Adv. Meteorol. 2016, 2016, 1835487. [Google Scholar] [CrossRef]
  4. Li, Z.; Tang, B.; Wu, H.; Ren, H.; Yan, G.; Wan, Z.; Trigo, I.F.; Sobrino, J.A. Satellite-derived land surface temperature: Current status and perspectives. Remote Sens. Environ. 2013, 131, 14–37. [Google Scholar] [CrossRef] [Green Version]
  5. Wan, Z. New refinements and validation of the collection-6 modis land-surface temperature/emissivity product. Remote Sens. Environ. 2014, 140, 36–45. [Google Scholar] [CrossRef]
  6. Guillemot, C.; Le Meur, O. Image inpainting: Overview and recent advances. IEEE Signal Proc. Mag. 2014, 31, 127–144. [Google Scholar] [CrossRef]
  7. Neteler, M. Estimating daily land surface temperatures in mountainous environments by reconstructed modis lst data. Remote Sens. 2010, 2, 333. [Google Scholar] [CrossRef]
  8. Fan, X.; Liu, H.; Liu, G.; Li, S. Reconstruction of modis land-surface temperature in a flat terrain and fragmented landscape. Int. J. Remote Sens. 2014, 35, 7857–7877. [Google Scholar] [CrossRef]
  9. Ke, L.; Ding, X.; Song, C. Reconstruction of time-series modis lst in central qinghai-tibet plateau using geostatistical approach. IEEE Geosci. Remote Sens. Lett. 2013, 10, 1602–1606. [Google Scholar] [CrossRef]
  10. Goovaerts, P. Geostatistics for Natural Resources Evaluation; Oxford University Press: Oxford, UK, 1997. [Google Scholar]
  11. Xu, Y.; Shen, Y. Reconstruction of the land surface temperature time series using harmonic analysis. Comput. Geosci. 2013, 61, 126–132. [Google Scholar] [CrossRef]
  12. Menenti, M.; Azzali, S.; Verhoef, W.; van Swol, R. Mapping agroecological zones and time lag in vegetation growth by means of fourier analysis of time series of ndvi images. Adv. Space Res. 1993, 13, 233–237. [Google Scholar] [CrossRef]
  13. Na, F.; Gaodi, X.; Wenhua, L.; Yajing, Z.; Changshun, Z.; Na, L. Mapping air temperature in the lancang river basin using the reconstructed modis lst data. J. Res. Ecol. 2014, 5, 253–262. [Google Scholar] [CrossRef]
  14. Savitzky, A.; Golay, M.J. Smoothing and differentiation of data by simplified least squares procedures. Anal. Chem. 1964, 36, 1627–1639. [Google Scholar] [CrossRef]
  15. Yu, W.; Nan, Z.; Wang, Z.; Chen, H.; Wu, T.; Zhao, L. An effective interpolation method for modis land surface temperature on the qinghai–tibet plateau. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2015, 8, 4539–4550. [Google Scholar] [CrossRef]
  16. Zeng, C.; Shen, H.; Zhong, M.; Zhang, L.; Wu, P. Reconstructing modis lst based on multitemporal classification and robust regression. IEEE Geosci. Remote Sens. Lett. 2015, 12, 512–516. [Google Scholar] [CrossRef]
  17. Sun, L.; Chen, Z.; Gao, F.; Anderson, M.; Song, L.; Wang, L.; Hu, B.; Yang, Y. Reconstructing daily clear-sky land surface temperature for cloudy regions from modis data. Comput. Geosci. 2017, 105, 10–20. [Google Scholar] [CrossRef]
  18. Jin, M. Interpolation of surface radiative temperature measured from polar orbiting satellites to a diurnal cycle: 2. Cloudy-pixel treatment. J. Geophys. Res. 2000, 105, 4061. [Google Scholar] [CrossRef]
  19. Lu, L.; Venus, V.; Skidmore, A.; Wang, T.; Luo, G. Estimating land-surface temperature under clouds using msg/seviri observations. Int. J. Appl. Earth Obs. Geoinf. 2011, 13, 265–276. [Google Scholar] [CrossRef]
  20. Zhang, X.; Pang, J.; Li, L. Estimation of land surface temperature under cloudy skies using combined diurnal solar radiation and surface temperature evolution. Remote Sens. 2015, 7, 905–921. [Google Scholar] [CrossRef]
  21. Prakash, S.; Norouzi, H.; Azarderakhsh, M.; Blake, R.; Prigent, C.; Khanbilvardi, R. Estimation of consistent global microwave land surface emissivity from amsr-e and amsr2 observations. J. Appl. Meteorol. Clim. 2018, 57, 907–919. [Google Scholar] [CrossRef]
  22. Prakash, S.; Norouzi, H.; Azarderakhsh, M.; Blake, R.; Tesfagiorgis, K. Global land surface emissivity estimation from amsr2 observations. IEEE Geosci. Remote Sens. Lett. 2016, 13, 1270–1274. [Google Scholar] [CrossRef]
  23. Prigent, C.; Jimenez, C.; Aires, F. Toward “all weather” long record, and real-time land surface temperature retrievals from microwave satellite observations. J. Geophys. Res. Atmos. 2016, 121, 5699–5717. [Google Scholar] [CrossRef]
  24. André, C.; Ottlé, C.; Royer, A.; Maignan, F. Land surface temperature retrieval over circumpolar arctic using ssm/i–ssmis and modis data. Remote Sens. Environ. 2015, 162, 1–10. [Google Scholar] [CrossRef]
  25. Kou, X.; Jiang, L.; Bo, Y.; Yan, S.; Chai, L. Estimation of land surface temperature through blending modis and amsr-e data with the bayesian maximum entropy method. Remote Sens. 2016, 8, 105. [Google Scholar] [CrossRef]
  26. Burrough, P.A.; McDonell, R.A. Principles of Geographical Information Systems; Oxford University Press: New York, NY, USA, 1998. [Google Scholar]
  27. Liang, S. Quantitative Remote Sensing of Land Surfaces; John Wiley & Sons: Hoboken, NJ, USA, 2005; Volume 30. [Google Scholar]
  28. Wang, K.; Liang, S. Evaluation of aster and modis land surface temperature and emissivity products using long-term surface longwave radiation observations at surfrad sites. Remote Sens. Environ. 2009, 113, 1556–1565. [Google Scholar] [CrossRef]
  29. Ogawa, K.; Schmugge, T.; Jacob, F.; French, A. Estimation of land surface window (8–12 μm) emissivity from multi-spectral thermal infrared remote sensing—A case study in a part of sahara desert. Geophys. Res. Lett. 2003, 30, 1067. [Google Scholar] [CrossRef]
  30. Kang, J.; Jin, R.; Li, X.; Zhang, Y.; Zhu, Z. Spatial upscaling of sparse soil moisture observations based on ridge regression. Remote Sens. 2018, 10, 192. [Google Scholar] [CrossRef]
  31. Hoerl, A.E.; Kennard, R.W. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics. 1970, 12, 55–67. [Google Scholar] [CrossRef]
  32. Yu, W.; Ma, M.; Wang, X.; Geng, L.; Tan, J.; Shi, J. Evaluation of modis lst products using longwave radiation ground measurements in the northern arid region of china. Remote Sens. 2014, 6, 11494. [Google Scholar] [CrossRef]
Figure 1. Study area and locations of ground-based stations.
Figure 1. Study area and locations of ground-based stations.
Remotesensing 10 01112 g001
Figure 2. Spatial and temporal distribution of missing moderate resolution imaging spectroradiometer (MODIS) Aqua land surface temperature (LST) data.
Figure 2. Spatial and temporal distribution of missing moderate resolution imaging spectroradiometer (MODIS) Aqua land surface temperature (LST) data.
Remotesensing 10 01112 g002
Figure 3. Influence of bias and correlation between the temporal predictor and target LST data on the prediction accuracy.
Figure 3. Influence of bias and correlation between the temporal predictor and target LST data on the prediction accuracy.
Remotesensing 10 01112 g003
Figure 4. Correlation at different directions and distances between MODIS Aqua LST data at the reconstructed time and themselves and MODIS LST data at other times.
Figure 4. Correlation at different directions and distances between MODIS Aqua LST data at the reconstructed time and themselves and MODIS LST data at other times.
Remotesensing 10 01112 g004
Figure 5. Average temporal absolute bias between one pixel from the MODIS Aqua LST data at the reconstructed time and all pixels from themselves and MODIS LST data at other times.
Figure 5. Average temporal absolute bias between one pixel from the MODIS Aqua LST data at the reconstructed time and all pixels from themselves and MODIS LST data at other times.
Remotesensing 10 01112 g005
Figure 6. Monthly reconstruction rates of MODIS Aqua LST data.
Figure 6. Monthly reconstruction rates of MODIS Aqua LST data.
Remotesensing 10 01112 g006
Figure 7. Examples of reconstructed MODIS Aqua LST data with different missing rates in space. The reconstructed data are within the red polygon.
Figure 7. Examples of reconstructed MODIS Aqua LST data with different missing rates in space. The reconstructed data are within the red polygon.
Remotesensing 10 01112 g007
Figure 8. The differences in monthly Moran’s I index between the known and reconstructed data.
Figure 8. The differences in monthly Moran’s I index between the known and reconstructed data.
Remotesensing 10 01112 g008
Figure 9. Scatter diagrams between the reconstructed data and ground-based observations at night.
Figure 9. Scatter diagrams between the reconstructed data and ground-based observations at night.
Remotesensing 10 01112 g009
Figure 10. Scatter diagrams between the reconstructed data and ground-based observations during the day.
Figure 10. Scatter diagrams between the reconstructed data and ground-based observations during the day.
Remotesensing 10 01112 g010
Figure 11. Comparisons among five reconstruction cases employing different predictor data.
Figure 11. Comparisons among five reconstruction cases employing different predictor data.
Remotesensing 10 01112 g011
Table 1. Remote sensing data information.
Table 1. Remote sensing data information.
Data SourceVariablePassing Time (Local Time)Spatial Resolution Temporal ResolutionDateVersion
MODIS TerraLST (MOD11A1)10:30/22:301 kmDailyJanuary 2014–December 20146.0
MODIS AquaLST (MYD11A1)01:30/13:301 km
NDVI (MYD13A2)
Emissivity (MYD11B1)5 km
AMSR236 GHZ TB (V/H)01:30/13:3010 km2.2
89 GHZ TB (V/H)
SRTMElevation--90 m----3.0
Table 2. Ground-based station information.
Table 2. Ground-based station information.
StationLongitudeLatitudeElevationLand CoverDateTemporal ResolutionVariable
ARC100.464°38.047°3033 mAlpine meadowJanuary 2014–December 201410 minUpwelling and downwelling longwave radiation
ARS100.411°37.984°3536 mAlpine meadow
JYL101.116°37.838°3750 mAlpine meadow
HZS100.192°38.225°2612 mCropland
HCG100.731°38.003°3137 mAlpine meadow
EBZ100.915°37.949°3294 mAlpine meadow
Table 3. Evaluations of the reconstructed data.
Table 3. Evaluations of the reconstructed data.
Station ε r e ε p r o d u c t ε m o d e l i n g ε c l o u d y s k y r k n o w n r r e N m i s s i n g N r e
Aqua Night (01:30)
ARC6.404.151.484.650.950.87164149
ARS4.934.071.352.430.920.88170155
JYL5.794.721.343.080.920.88182176
HZS4.442.761.073.320.970.92184169
HCG5.784.311.303.620.930.88165150
EBZ6.044.491.443.780.930.88172157
Aqua Day (13:30)
ARC6.714.551.294.760.920.92265206
ARS5.943.851.294.340.890.88301244
JYL8.816.901.905.130.850.83324269
HZS6.804.001.715.220.940.80232170
HCG10.666.671.388.200.870.85265201
EBZ8.335.391.456.180.870.86287226

Share and Cite

MDPI and ACS Style

Kang, J.; Tan, J.; Jin, R.; Li, X.; Zhang, Y. Reconstruction of MODIS Land Surface Temperature Products Based on Multi-Temporal Information. Remote Sens. 2018, 10, 1112. https://doi.org/10.3390/rs10071112

AMA Style

Kang J, Tan J, Jin R, Li X, Zhang Y. Reconstruction of MODIS Land Surface Temperature Products Based on Multi-Temporal Information. Remote Sensing. 2018; 10(7):1112. https://doi.org/10.3390/rs10071112

Chicago/Turabian Style

Kang, Jian, Junlei Tan, Rui Jin, Xin Li, and Yang Zhang. 2018. "Reconstruction of MODIS Land Surface Temperature Products Based on Multi-Temporal Information" Remote Sensing 10, no. 7: 1112. https://doi.org/10.3390/rs10071112

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop