Estimating Daily Maximum and Minimum Land Air Surface Temperature Using MODIS Land Surface Temperature Data and Ground Truth Data in Northern Vietnam

This study aims to evaluate quantitatively the land surface temperature (LST) derived from MODIS (Moderate Resolution Imaging Spectroradiometer) MOD11A1 and MYD11A1 Collection 5 products for daily land air surface temperature (Ta) estimation over a mountainous region in northern Vietnam. The main objective is to estimate maximum and minimum Ta (Ta-max and Ta-min) using both TERRA and AQUA MODIS LST products (daytime and nighttime) and auxiliary data, solving the discontinuity problem of ground measurements. There exist no studies about Vietnam that have integrated both TERRA and AQUA LST of daytime and nighttime for Ta estimation (using four MODIS LST datasets). In addition, to find out which variables are the most effective to describe the differences between LST and Ta, we have tested several popular methods, such as: the Pearson correlation coefficient, stepwise, Bayesian information criterion (BIC), adjusted R-squared and the principal component analysis (PCA) of 14 variables (including: LST products (four variables), NDVI, elevation, latitude, longitude, day length in hours, Julian day and four variables of the view zenith angle), and then, we applied nine models for Ta-max estimation and nine models for Ta-min estimation. The results showed that the differences between MODIS LST and ground truth temperature derived from 15 climate stations are time and regional topography dependent. The best results for Ta-max and Ta-min estimation were achieved when we combined both LST daytime and nighttime of TERRA and AQUA and data from the topography analysis.


Introduction
Land air surface temperature (T a , also called "air temperature" or "near surface temperature") data are usually collected as point data from weather station locations, typically at 2 m above the land surface. It is an important parameter in a wide range of fields, such as agriculture, e.g., crop evapotranspiration [1], crop yield prediction [2,3], hydrology [4,5], ecology, environment and climate change [6,7]. Generally, T a values obtained from weather stations have a very high accuracy and temporal resolution [8], but do not capture information for a whole region and may therefore be unsuitable for spatial modelling applications [9][10][11].
In order to obtain T a information for a region, researchers have proposed various methods of interpolation based on known weather station sites [12][13][14]. These interpolation methods' accuracy is highly dependent on the weather station network density, as well as the scale of spatial and temporal variability [15,16]. Furthermore, station geometry and topography (elevation) change also affects the accuracy of interpolation, especially in regions with a wide range of elevation [17,18]. However, the spatial distribution of weather stations is often limited in developing countries. Our study area of Vietnam has 170 surface meteorological observing stations, including 97 synoptic and 26 international exchange stations [19], which is obviously inadequate for a country with an area of 331,688 km 2 in which about 40% is mountainous, 40% hill and the remaining 20% lowland. Therefore, interpolation techniques may not be suitable for Vietnam.
Fortunately, remote sensing data provide a promising solution to overcome the limitation of interpolation techniques in mountainous areas and a sparsity of weather station areas. The successful launch of the Advanced Very High-Resolution Radiometer (AVHRR) in 1981 and the Moderate Resolution Imaging Spectroradiometer (MODIS) on board TERRA (December 1999) and AQUA (May 2002) has driven researchers to study new satellite-based methods, as a hot topic in recent years [16,[20][21][22][23][24][25].
In recent years, there have been more and more studies employing land surface temperature (LST) obtained from remotely-sensed images for T a estimation because of high spatial and temporal resolution, free availability and easy access. Particularly, MODIS on board TERRA and AQUA can provide daily LST data with high temporal (four times per day, TERRA LST daytime, TERRA LST nighttime, AQUA LST daytime, AQUA LST nighttime, which overpass local time at around 10:30 a.m., 10:30 p.m., 1:30 a.m. and 1:30 p.m., respectively) and very high spatial resolution (1 km) are widely applied. The difference between LST and T a is strongly influenced by the surface characteristics and atmospheric conditions [26,27]. In some regions, the difference between LST and T a is high [17,28]. However, researchers from all over the world state that there is a strong linear correlation between MODIS LST and T a over many regions, e.g., in Africa [15], in Portugal [29], over the U.S. [30,31] and in Southeast Asia [32,33]. The detailed information of this difference, as well as the possible causes of this difference are still limited and need to be studied.
Some researchers [29,31,34] reviewed the types of commonly-used methods for T a estimation based on LST. There are three main distinct types of methods: The first type is the temperature-vegetation index method (TVX), which is based on the assumption that in an absolutely thick canopy, the temperature at the top of the canopy is the same as within the canopy [35]; and there is a strong negative correlation between LST and the vegetation index, such as NDVI [9,10,[36][37][38][39]. However, in some cases, this method is not satisfying due to the assumption that it often does not fit to the reality or the effect of seasonality, land cover type or soil moisture [15,40].
The second type includes the surface energy-balance-based methods. The sum of in-coming net radiation and anthropogenic heat fluxes is considered equal to the sum of the surface's sensible and latent heat fluxes [41]. The major drawback of these methods is that they require large amounts of information often not provided by remote sensing [24,25].
The last type is using statistical methods that are based on regression techniques. These methods include various levels of complexity, from basic approaches that only use LST for T a estimation [16,25] to advanced approaches that use more than one independent variable, such as elevation, NDVI, land cover, distance to water body, solar zenith angle, day length in hours, latitude and altitude [15,29,32,[42][43][44]. One of the biggest advantages of this method is that the systematic regional errors in satellite data can be reduced [45].
The most recently popular studies of T a estimation using statistic approaches are shown in Table 1. However, most of these studies have only used LST daytime and LST nighttime solely for T a maximum (T a-max ) and T a minimum (T a-min ) estimation, respectively. In a recent study [31], both LST nighttime and daytime were used for T a-max and for T a-min estimation. However, this study was only applied for the growing season (May-September) from 2008-2012, and the elevation of the study site (the Corn Belt region of the Midwestern U.S.) ranging from 87-666 m and mostly covered by crops has a small vegetation index range. There are several researchers who have studied the effect of the time of observation on the relationship between LST and T a , but the conclusions are quite different in various time and geographical regions of the study areas. For instance, [25] found that the overpass time of TERRA and AQUA has little impact on the accuracy of T a estimation in Mississippi State. Zhang et al. [39] concluded that for daily T a estimation, TERRA LST and AQUA LST give the same results. Benali et al. [29] showed that the use of both AQUA LST daytime and LST nighttime gives better accuracy of T a-max and T a-min estimation, respectively, than TERRAs in Portugal. In contrast, [38] stated that TERRA LST daytime and TERRA nighttime were better than AQUA LST daytime and nighttime for T a estimation. These differences could be understood as not only the time of observation, but also geographical location affecting the relationship between LST products and T a and, therefore, affecting the accuracy estimation of T a based on LST products. Because of all of these different conclusions about the relationship between T a and LST of TERRA and AQUA, in this study, we analyzed the relationship between T a with both TERRA LST and AQUA LST products.
In addition, as far as we know, there are no studies over Vietnam that have integrated both TERRA and AQUA LST of daytime and nighttime for T a estimation (using all four MODIS LST datasets).
The main objective of this research was to estimate daily T a (maximum and minimum) using both TERRA and AQUA MODIS LST products (daytime and nighttime) and auxiliary data, solving the discontinuity problem of ground measurements. In addition, to find out which variables, among 14 predefined variables, are the most effective to describe the relationship between LST and T a , we have tested several methods, such as: the Pearson correlation coefficient, forward selection, backward elimination, stepwise, Bayesian information criterion (BIC), adjusted R-squared and the principal component analysis (PCA) of 14 variables (including: LST products (four variables), NDVI, elevation, latitude, longitude, day length in hours, Julian day and four variables of the view zenith angle). Finally, we applied nine models for T a-max and nine models for T a-min estimation.

Study Area
The study area is located in northern Vietnam and covers more than 37,000 km 2 , comprising the provinces of Hoa Binh, Ha Noi, Vinh Phuc, Thai Nguyen, Yen Bai, Phu Tho and Son La (see Figure 1). The study area extends from 20 • 18 N-22 • 40 N and from 103 • 12 E-106 • 18 E. The elevation ranges from sea level to over 3000 m.
The northern Vietnam study site was chosen for the following reasons: (1) Good distribution of meteorological station network in comparison to southern Vietnam.
(2) Wide range of elevation (from approximately sea level to more than 3000 m).
(3) The spatial heterogeneity of land cover. (4) There is no study about T a estimation using all 4 MODIS LST datasets in northern Vietnam.
In this area, the topography is quite complex, increasing from southeast to northwest, and mostly divided into two regions: lowland and highland. The low part includes: Hoabinh, Hanoi, Phutho, Vinhphuc and Thainguyen provinces. The high part includes 2 large provinces: a part of Yenbai and Sonla. (3) The spatial heterogeneity of land cover. (4) There is no study about Ta estimation using all 4 MODIS LST datasets in northern Vietnam.
In this area, the topography is quite complex, increasing from southeast to northwest, and mostly divided into two regions: lowland and highland. The low part includes: Hoabinh, Hanoi, Phutho, Vinhphuc and Thainguyen provinces. The high part includes 2 large provinces: a part of Yenbai and Sonla. The climate and weather in this area also vary depending on the elevation and type of landscape. The humidity is high, with the average ranging around 84% a year.
In late October and early November 2008, there was torrential rain in northern and central Vietnam that triggered severe floods in this region [49].
The location and elevation of meteorological stations are shown in Table 2.  The climate and weather in this area also vary depending on the elevation and type of landscape. The humidity is high, with the average ranging around 84% a year.
In late October and early November 2008, there was torrential rain in northern and central Vietnam that triggered severe floods in this region [49].
The location and elevation of meteorological stations are shown in Table 2. Two MODIS LST products (h27v06, Collection 5, from 2003-2013, over northern Vietnam) were used in this study: MOD11A1 daily land surface temperature and emissivity from the TERRA satellites and MYD11A1 daily land surface temperature and emissivity from the AQUA satellites. There exist 4 LST data records per day, two from the TERRA satellites and two others from the AQUA satellites, which pass over the study site (local solar time) mostly around 10:30 a.m., p.m. and 1:45 a.m. and p.m., respectively. These times are relatively close to the times of maximum and minimum T a daily data. In total, there were more than 8000 images (in HDF format, from 1 January 2003-31 December 2013) downloaded from the U.S. Geological Survey [50]. The MODIS LST is generated using a split-window algorithm [51] with two thermal infrared bands, 31 (10.78-11.28 µm) and 32 (11.77-12.27 µm). These two products (MOD11A1 and MYD11A1) have been validated, and the accuracy was reported better than 1 K under clear sky conditions [52].

Elevation
In this study area, the elevation is quite complex; it ranges from sea level to above 3000 m (see Figure 1). Generally, higher elevations are associated with lower temperatures. Elevation data were obtained from ASTER Global DEM. This data are available from the U.S. Geological Survey (USGS) with a spatial resolution of 30 m. These elevation data were resized to 1-km resolution using the nearest neighbor resampling type in order to be associated with other data resolutions, such as MODIS LST or NDVI data.

Vegetation Based on NDVI
Several studies have shown that the vegetation cover on the surface can affect LST values [53]. Vegetation cover (NDVI) data are provided every 16 days at 1-kilometer spatial resolution from the MODIS satellite, including MOD13A2 (TERRA 16-days period starting Day 001) and MYD13A2 (AQUA 16-day period starting Day 009). There are some studies that use both MOD13A2 and MYD13A2 in order to composite 8-day period data by averaging these two products [18]. Miura and Yoshioka [54] stated that MOD13A2 and MYD13A2 could be interchangeable. We assumed that NDVI did not change within 16 days. Therefore, in this study, to correct the influence of vegetation, we used NDVI 16 days data at 1-km resolution obtained from the MOD13A2 product.

Meteorological Data
Daily maximum and minimum air temperature data have been collected from 15 meteorological stations in northern Vietnam (see Figure 1), from 2003-2013. The data were obtained from the Vietnam Institute of Meteorology, Hydrology and ENvironment (IMHEN).

Auxiliary Data
Day length (Dlth) is the total time that any portion of the Sun is above the horizon. Typically (for low and mid-latitude locations), this will be the elapsed time beginning at sunrise and ending at sunset [55]. In this study, the day length (in hours) was downloaded from the Astronomical Applications Department of the U.S. Naval Observatory website [55].
The Julian day (Jd) was extracted from NASA server [56]. Both Julian day and day length are proxies for the fraction of solar energy absorption during the day and emitted energy during the night, influencing the diurnal amplitude of T a throughout the year. Jang et al. [57] showed that the Julian day was a more significant parameter than altitude or the solar zenith angle in the inter-seasonal estimation of T a . In order to understand the effect of viewing angle in the temporal variations in LST, especially in rugged regions, the view zenith angle (VZA) of daytime and nighttime should be taken into consideration [58]. In this study, there are four types of VZA of TERRA daytime, TERRA nighttime, AQUA daytime, AQUA nighttime (VZA td , VZA tn , VZA ad , VZA an , respectively). These data were collected from MOD11A1 and MYD11A1 products. We also took the effects of the latitude, longitude and elevation of each station, collected from IMHEN, into consideration to estimate T a accuracy.
All key terms and their explanations are summarized in Table 3. Table 3. Description of the key terminology used in this study.

Used Terms Description
Daily minimum T a estimation T a-max-es (

Preprocessing Data
There were 8036 files of MODIS HDF format (MOD11A1 and MYD11A1, Version 5) from 1 January 2003-31 December 2013 that were downloaded for this study. Because all MODIS data were downloaded from USGS in HDF format (Hierarchical Data Format), with each file containing 12 data layers [58], we used the MODIS Reprojection Tool to extract the corresponding bands (LST_Day_1km, LST_Night_1km and view zenith angle) from MOD11A1 and MYD11A1. Next, we used the layer stacking tool in Envi5.0 to stack images for time series analysis purposes, before clipping the images to fit to our study area using ArcGIS 10.1.

LST Data Retrieval at Weather Stations
LST data under clear sky conditions at weather stations were retrieved by the following steps: (1) MODIS LST day and night was extracted from MOD11A1 and MYD11A1 data for the pixel in which the meteorological stations are located. It should be noted that in the MODIS LST (v005) product, LST values derived from a single clear-sky MODIS observation are selected from MOD11_L2 files if the viewing zenith angles are small or larger zenith angles have LST errors smaller than 2 K [58]. This means that MOD11A1 and MYD11A1 LST were not produced for pixels that are classed as cloudy or cloud contamination remaining.
(2) However, [58] also noted that pixels that are lightly or modestly cloud contaminated are not removed by using cloud-removing masks if the contamination is smaller than the normal temporal variations in clear-sky LST. Typically, these are thin and subpixel clouds that are difficult to detect by the algorithm. In this case, we applied two steps to remove those types of errors. First, we simply filter and remove all unrealistic LST data that had values greater than 100 • C and below −50 • C. Second, we calculated the difference between T a-max versus LST daytime and T a-min versus LST nighttime. Then, we applied a statistic outlier removal based on these differences' histograms to detect and remove data with unusually large differences.
(3) All valid LST day and nighttime data of TERRA and AQUA were statistical compared with T a-max and T a-min , respectively.

Estimation of Land Air Surface Temperature Using MODIS LST and Auxiliary Data
The method for this study was based on a multivariate linear regression analysis. This method was chosen because it could be applied for sparse meteorological station networks like northern Vietnam. Although some interpolation methods may give a higher accuracy result, they are not possible for regions with poorly-distributed station networks [59]; or physical models [41] that require an unreliable amount of input data.

Variable Selection
Zaksek and Schroedter-Homscheidt [34] stated that T a is driven more by LST than by direct solar radiation, meaning that LST is the most important variable for T a estimation; in reference to previous studies [15,29,31,42,57], also in consideration of all of the available data that potentially have an effect on the accuracy of T a estimation of the study area. We have collected 14 variables, including: LST products (4 variables), NDVI, elevation (DEM), latitude, longitude, day length in hours and Julian day, and 4 variables of the view zenith angle as the potential variables (candidate variables) for T a estimation.
Besides the 4 LST variables, we chose NDVI because it influences the land surface vegetation properties. Elevation, latitude and longitude were chosen for capturing the variability of climatic conditions between different regions. Day length in hours was chosen because we considered that it would affect the received solar radiation. Following [45,57], we chose Julian day because it reflects seasonal variation in air temperature. The view zenith angle was considered as a factor affecting the accuracy of LST data.
Since we have 14 potential variables (predictors), there are 2 14 (= 16,384) possible models that could be applied for T a-max and T a-min estimation. Clearly, we cannot apply all of these possible models and test them one by one. Consequently, to save the time of data processing and to find the smallest model that best fits the data, we selected the best subset of variables among the 14 available variables.
In this study, we used several popular methods for variable selection, such as: (1)

Pearson Correlation Coefficient
To select the variables for the multivariate linear regression models, we calculated the Pearson's correlation coefficient (r) of all related variables. According to [11], r close to ±1 indicates a strong correlation; r = 0 means that there is no correlation; 0.25 ≤ r ≤ 0.75, means there is a moderate degree of correlation and 0 ≤ r < 0.25, a weak correlation. Based on this, we took only variables with r greater than 0.25 for T a estimation.
The results in Table 4 show that we should use 7 variables for T a-min estimated models (4 LST, Ele, Long and Dlth); and 7 other variables for T a-max estimated models (4 LST, Ele, Dlth and Jd).
Variable Selection Based on Adjusted R-Squared (R 2 ) and BIC The best model will be the one with the largest value of R 2 : Minimizing the BIC is intended to give the best model. N is the number of observations; k is the number of predictors (variables); r 2 is the coefficient of determination; and SSE is the sum of squared error.
In Figure 2, colored parts (green for R 2 and orange for BIC) indicate that a variable was included in the model and white color that a variable was not included in the model. All of these models include the intercept. Minimizing the BIC is intended to give the best model. N is the number of observations; k is the number of predictors (variables); r 2 is the coefficient of determination; and SSE is the sum of squared error.
In Figure 2, colored parts (green for R 2 and orange for BIC) indicate that a variable was included in the model and white color that a variable was not included in the model. All of these models include the intercept.  All variables are described in Table 3.
Looking at the y-axis of Figure 2, it can be clearly seen that, in both cases of Ta-max and Ta-min, the top three models have the same results of R 2 and BIC.
In this case, we chose all top 3 models for Ta estimation. These model forms are as follows: Ta-min = a × LSTtd + b × LSTtn + c × LSTad + d × LSTan + e × Dlth + f × Jd + g.
Variable selection using forward, backward and stepwise Another popular method that is used for variable selection is the stepwise method. In this method, there are three types of operative steps: forward, backward and stepwise.  All variables are described in Table 3.
Looking at the y-axis of Figure 2, it can be clearly seen that, in both cases of T a-max and T a-min , the top three models have the same results of R 2 and BIC.
In this case, we chose all top 3 models for T a estimation. These model forms are as follows: T a-max = a × LST td + b × LST tn + c × LST ad + d × LST an + e × Ele + f × Long + g, T a-min = a × LST td + b × LST tn + c × LST ad + d × LST an + e × Dlth + f × Lat + g × Jd +h × Ele + i, T a-min = a × LST td + b × LST tn + c × LST ad + d × LST an + e × Dlth + f × Lat + g × Jd + h, T a-min = a × LST td + b × LST tn + c × LST ad + d × LST an + e × Dlth + f × Jd + g.
Variable Selection Using Forward, Backward and Stepwise Another popular method that is used for variable selection is the stepwise method. In this method, there are three types of operative steps: forward, backward and stepwise.
Forward selection starts with no variable in the model (intercept only model). In the next steps, the most significant variables are added to the model one by one. The process stops when all of the variables not in the model have a p-value greater than 0.15.
Backward elimination starts with the model including all variables. In the next step, the least significant variable will be removed. The procedure continues until all of the variables in the model have p-values less than or equal to 0.15.
The stepwise method adds or removes a variable in each step, depending on its p-value. This process continues until all variables within the model have a p-value ≤0.15, and all of the variables that were not in the model have a p-value >0.15.
Based on the results of forward, backward and stepwise selection, the T a-max estimated model should include 12 variables (using all variables except VZA tn and VZA an ); the T a-min estimated model should include 10 variables (except long, VZA an , VZA td , VZA tn ).
The T a-max and T a-min estimated models were as follows: Variable Selection-Based Principal Component Analysis (PCA) The result of the PCA showed that both T a-max and T a-min estimation models should include 4 LST, elevation, day length in hours and Julian day.
From all of the analyses above (Pearson correlation, R 2 , BIC, stepwise and PCA), it is indicated that 4 LSTs (LST td , LST tn , LST ad , LST an ) play an important role in T a-max and T a-min estimation, because they are always included in the top models of all above-mentioned methods.

From Equations
We evaluated the models using the coefficient of determination (r 2 ), the root mean square error (RMSE) and the mean absolute error (MAE). Because RMSE is very sensitive to outliers, hence MAE was chosen as an additional measure of the model quality.
where T a,i is the observed land air surface temperature from weather stations and T es,i is the corresponding land air surface temperature estimated using linear regression analysis methods.

Model Calibration and Validation
For calibration and validation purposes, all data were randomly divided into two parts: calibration and validation datasets using a 70%/30% proportion, respectively. To calculate the model coefficients a-m (see Table 5) for the models, we applied a least squares fitting analysis to the calibration dataset. The constant coefficients are listed in Appendixs A and B. To validate the models, we applied the constant coefficients, which were calculated based on the calibration dataset, to the validate dataset.

The Relationship between MODIS LST and T a
In order to analyze the relations between T a and all MODIS LST data, we calculated the coefficient of determination (r 2 ), root mean square error (RMSE) and mean absolute error (MAE) of each type of MODIS LST (TERRA daytime, TERRA nighttime, AQUA daytime and AQUA nighttime) solely with T a measured from weather stations. Figure 3a shows that LST nighttime of both TERRA and AQUA showed a stronger correlation with T a (both T a-max and T a-min ) than for LST daytime. This can be explained by the fact that during the night, solar radiation does not affect the thermal infrared signal. Looking at the overpasses time of the satellite, the overpass time of AQUA daytime (around 1:45 p.m.) is later than the overpass time of TERRA (10:30 a.m.), as more solar radiation had been received by the land surface. As a result, LST AQUA daytime is hotter than LST TERRA daytime. Similarly, the overpass time of TERRA LST nighttime is around 10.30 p.m., three hours earlier than that of AQUA (1.45 a.m.), and the LST nighttime of TERRA is slightly higher than that of AQUA. This could be reflected by the cooling process of the land's surface at night.
It also showed that LST nighttime of TERRA seems to perform better than LST nighttime of AQUA for T a-max estimation. In contrast, LST nighttime of AQUA seemed better than LST nighttime of TERRA for T a-min estimation.
In order to analyze which of the four LST variables is the best suited T a estimation, we used the adjusted R-squared and BIC criteria (Figure 3b). Among the resulting models, we chose the one that only contained one variable. This variable is as such the best to estimate T a .
It can be clearly seen that, based on the adjusted R-squared and BIC criteria (Figure 3b), the LST nighttime of TERRA and of AQUA were the best variables for T a-max and T a-min estimation, respectively. This result was consistent with the result shown in Figure 3a.
This result also indicates that the overpass time of satellites and the time of T a-max , T a-min recorded had no key influence on the relationship between T a and LST. This is consistent with other research [8,25].

Ta-Max Estimation
The parameters of models for T a-max estimation were determined when we applied Models 1-9 to the calibration dataset. This result is shown in Appendix A. Figure 4 show that the results of the nine models were consistent with the processing in Section 2.5. Model 1 showed the lowest result, followed by Models 2-4. Models 5-8 presented similar high results. However, Model 5 used only six variables; Models 6, 7 and 8 used 7, 8 and 12 variables, respectively. Based on this point, Model 5 would be chosen as the best model of T a-max estimation.
Combining LST TERRA nighttime and LST AQUA daytime (Model 2), the result of T a-max estimation was significantly better than Model 1 using LST TERRA nighttime solely. However, when we added LST TERRA daytime (Model 3, three variables), LST TERRA daytime and LST AQUA nighttime (Model 4, four variables), the performance of these models was not significantly improved in comparison to Model 2 (using two variables).
Comparing Model 4 and Model 5, it can be seen that, by combining four LSTs and elevation (Ele), longitude (Long) into the model, the result was much higher than using four LSTs only. This indicates that elevation (Ele), as well as the location (Long) of the weather station play an important role in T a-max estimation.

Ta-Max Estimation
The parameters of models for Ta-max estimation were determined when we applied Models 1-9 to the calibration dataset. This result is shown in Appendix A. Figure 4 show that the results of the nine models were consistent with the processing in Section 2.5. Model 1 showed the lowest result, followed by Models 2-4. Models 5-8 presented similar high results. However, Model 5 used only six variables; Models 6, 7 and 8 used 7, 8 and 12 variables, respectively. Based on this point, Model 5 would be chosen as the best model of Ta-max estimation.
Combining LST TERRA nighttime and LST AQUA daytime (Model 2), the result of Ta-max estimation was significantly better than Model 1 using LST TERRA nighttime solely. However, when we added LST TERRA daytime (Model 3, three variables), LST TERRA daytime and LST AQUA nighttime (Model 4, four variables), the performance of these models was not significantly improved in comparison to Model 2 (using two variables).
Comparing Model 4 and Model 5, it can be seen that, by combining four LSTs and elevation (Ele), longitude (Long) into the model, the result was much higher than using four LSTs only. This indicates that elevation (Ele), as well as the location (Long) of the weather station play an important role in Ta-max estimation.

Ta-Min Estimation
The parameters of models for T a-min estimation were determined when we applied Models 10-18 to the calibration dataset. This result is shown in Appendix B.
In general, Figure 5 shows that all models gave similar results of T a-min estimation. From the simple model with two variables to the complex model (Model 17) using 10 variables, the results are similar. Figure 5 also shows that T a-min estimation can reach a very high accuracy (r 2 = 0.85, RMSE = 2.31, MAE = 1.75) when using only one variable: LST an . However, it was difficult to increase the accuracy of T a-min estimated even with 10 variables or all 14 variables (see Figure 5).

Ta-Min Estimation
The parameters of models for Ta-min estimation were determined when we applied Models 10-18 to the calibration dataset. This result is shown in Appendix B.
In general, Figure 5 shows that all models gave similar results of Ta-min estimation. From the simple model with two variables to the complex model (Model 17) using 10 variables, the results are similar. Figure 5 also shows that Ta-min estimation can reach a very high accuracy (r 2 = 0.85, RMSE = 2.31, MAE = 1.75) when using only one variable: LSTan. However, it was difficult to increase the accuracy of Ta-min estimated even with 10 variables or all 14 variables (see Figure 5).

Performance of the Best Model
In order to test the effect of weather station location, as well as the seasonal change or any other factor related to station characteristics, we applied the best model to all datasets. The results are shown in Figure 6.
Looking at Figures 4-6, it can be clearly seen that there are similar results of Model 5 and Model 15 when we applied them to the validation dataset or to all datasets. This consistent result indicates that there is no significant factor that could affect the accuracy of Ta estimation when Model 5 and Model 15 were applied for Ta-max and Ta-min estimation, respectively.

Performance of the Best Model
In order to test the effect of weather station location, as well as the seasonal change or any other factor related to station characteristics, we applied the best model to all datasets. The results are shown in Figure 6.
Looking at Figures 4-6, it can be clearly seen that there are similar results of Model 5 and Model 15 when we applied them to the validation dataset or to all datasets. This consistent result indicates that there is no significant factor that could affect the accuracy of T a estimation when Model 5 and Model 15 were applied for T a-max and T a-min estimation, respectively.

MODIS LST Products for Ta Estimation
There are two MODIS sensors, TERRA and AQUA, which provides LST data four times daily (i.e., LSTtd, LSTtn, LSTad, LSTan). Most previous studies have used LST daytime for Ta-max and LST nighttime for Ta-min estimation. Some studies used only TERRA MODIS LST for Ta estimation [45,46]. There is a handful studies using LST daytime for Ta-min and LST nighttime for Ta-max estimation [29,42].
As far as we know, there have also been several studies using both LST daytime and nighttime of both MODIS sensors on TERRA and AQUA for their Ta estimated models [15,33,42,60]. However, the purpose of using four times daily LST was for filling the missing LST value. In our study, it is required that all four LST data have to be available.
Zhang et al. [42] combine TERRA LST daytime and nighttime (two variables), AQUA LST daytime and nighttime (also two variables) for Ta estimation and concluded that: (i) nighttime LST was better than daytime for deriving daily Ta; and (ii) incorporating daytime and nighttime LST significantly improved the estimation of Ta, as compared to using LST nighttime or daytime solely. Our results were consistent with this study [42]; however, our analysis could provide an even deeper insight into the matter. Model 2 and Figure 3b show that the combination of TERRA LST nighttime and AQUA LST daytime is even better than that of TERRA LST daytime and TERRA LST nighttime.

MODIS LST Products for T a Estimation
There are two MODIS sensors, TERRA and AQUA, which provides LST data four times daily (i.e., LST td , LST tn , LST ad , LST an ). Most previous studies have used LST daytime for T a-max and LST nighttime for T a-min estimation. Some studies used only TERRA MODIS LST for T a estimation [45,46]. There is a handful studies using LST daytime for T a-min and LST nighttime for T a-max estimation [29,42].
As far as we know, there have also been several studies using both LST daytime and nighttime of both MODIS sensors on TERRA and AQUA for their T a estimated models [15,33,42,60]. However, the purpose of using four times daily LST was for filling the missing LST value. In our study, it is required that all four LST data have to be available.
Zhang et al. [42] combine TERRA LST daytime and nighttime (two variables), AQUA LST daytime and nighttime (also two variables) for T a estimation and concluded that: (i) nighttime LST was better than daytime for deriving daily T a ; and (ii) incorporating daytime and nighttime LST significantly improved the estimation of T a , as compared to using LST nighttime or daytime solely. Our results were consistent with this study [42]; however, our analysis could provide an even deeper insight into the matter. Model 2 and Figure 3b show that the combination of TERRA LST nighttime and AQUA LST daytime is even better than that of TERRA LST daytime and TERRA LST nighttime.
For T a-max estimation, Figure 4 shows that if only LST products are used (without auxiliary data) the best result was achieved when we combined all four LST data (model 4). However, the Model 4 result (using four LSTs) was just slightly better than Model 2 (using only LST tn and LST ad ).
Similarly, the best result of T a-min estimation was achieved when we combined all four LSTs (Model 13). The combination of LST tn and LST an (Model 11) made the result of T a-min estimation better than using LST an solely (Model 10).
As we can see from Figure 2, T a-min has a stronger correlation with LST nighttime of both TERRA and AQUA than T a-max . However, as the final result, T a estimation and T a-max estimation show a much higher accuracy than T a-min .

Effect of Seasonal on the Accuracy of T a Estimation
To examine the effect of seasons on the relationship between LST and T a (T a-max and T a-min ), we separated the data into four seasons: spring (March, April, May), summer (June, July, August), autumn (September, October, November) and winter (December, January, February).
Looking at the nine upper panels (T a-max estimation) of Figure 7, it can be seen that the accuracy of T a-max estimation was significantly improved through different models. Models 1, 2, 3 and 4 show a low accuracy of T a-max estimated in summer. Models 5, 6, 7 and 8 show a similar accuracy of T a-max estimation. Model 9 again shows a low accuracy in summer. This means that, if we estimate T a-max using LST data without other auxiliary data, the results would not be accurate. This can be explained by, as in the summer, the land surface receiving more solar radiation, and as such, the relation between T a and LST becomes more complex.
In general, the changing season had no significant effect on the accuracy of T a estimation if we use all four LST data and other variables for linear regression.

Effect of View Zenith Angle on the Accuracy of T a Estimation
As mentioned in Section 2.4, LST data were collected at smaller viewing zenith angles or LST retrievals at larger zenith angles, but with LST errors smaller than at least 2 K.
In order to check the effect of VZA on the relationship between MODIS LST and T a , we divided the data into three groups based on the range of VZA: bl40 (0 • ≤ VZA ≤ 40 • ), f41t90 (41 • ≤ VZA ≤ 90 • ) and f91t130 (91 • ≤ VZA ≤ 130 • ). Using all linear regression models in Table 5 for T a-max and T a-min estimation, we then calculated the difference between T a estimated and T a measured (difference = T a measured − T a estimated). These differences are shown in Figure 8. Although the differences might vary from Models 1-9 of T a-max estimation and Models 10-18 of T a-min estimation, it was similar between the three groups through all models. This indicates that the effect of VZA on the result of T a estimation was not significant.

Effect of Station Elevation on Accuracy
To test the effect of weather station elevation on the estimation results, we divided the data into three regions: bl200m, f200t500m and ov500m. The bl200m region includes stations that have an elevation below 200 m; the f200t500m region includes stations having an elevation from 200-500 m; and ov500m region includes stations that have an elevation higher than 500 m. It should be noted that this division is just one option to test the effect of station altitude that might affect the results. It could be divided into more parts for more detail (the more parts are divided, the more detail on the differences could be achieved). We calculated the difference between T a estimated and T a measured of these three groups of elevation (difference = T a measured − T a estimated). These differences are shown in Figure 9.
As we can see in Figure 9, for T a-max , the differences of three groups were variations from Models 1-9. The most significant was found in Model 1 and Model 9. Models 5, 6, 7 and 8 showed a similar difference between three groups. This indicates that we can reduce the effect of the station's elevation by using Model 5, 6, 7 or 8.
Looking at Table 5, we found that all of Models 5, 6, 7 and 8 include longitude and elevation variable in the models. Model 9 also has the elevation variable in it, but Figure 9 shows that there is a difference between the three elevation groups. It could be concluded that the longitude variable affects the T a-max estimation.
Regarding T a-min estimation, Figure 9 shows a similar result of all models (from Models 10-18) and a similar performance between the three elevation groups. This indicates that the station's elevation has no effect on T a-min estimation.
Looking at Table 5, we found that all of Models 5, 6, 7 and 8 include longitude and elevation variable in the models. Model 9 also has the elevation variable in it, but Figure 9 shows that there is a difference between the three elevation groups. It could be concluded that the longitude variable affects the Ta-max estimation.
Regarding Ta-min estimation, Figure 9 shows a similar result of all models (from Models 10-18) and a similar performance between the three elevation groups. This indicates that the station's elevation has no effect on Ta-min estimation.

For Ta-max Estimation
When comparing the results of Model 4 (using four LST data) versus Models 5-9 (using four LST data and auxiliary data; see Figure 4), it can be clearly seen that the results were significantly improved. The coefficient of determination (r 2 ) was increased from 0.88-0.93; RMSE and MAE were decreased from 1.91 down to 1.43 and 1.50 down to 1.08 °C respectively.
The performances of Models 5, 6 and 7 were similar and had very high accuracy. As shown in the variable selection section (Section 2.5.1), these models came from the top three models of adjusted R-squared and BIC criteria.

For T a-max Estimation
When comparing the results of Model 4 (using four LST data) versus Models 5-9 (using four LST data and auxiliary data; see Figure 4), it can be clearly seen that the results were significantly improved. The coefficient of determination (r 2 ) was increased from 0.88-0.93; RMSE and MAE were decreased from 1.91 down to 1.43 and 1.50 down to 1.08 • C respectively.
The performances of Models 5, 6 and 7 were similar and had very high accuracy. As shown in the variable selection section (Section 2.5.1), these models came from the top three models of adjusted R-squared and BIC criteria. Model 8, which was chosen from stepwise regression, showed a similar result of accuracy. However, this model used up to 12 variables.
Model 9, which was chosen based on PCA analysis, also showed a high accuracy result with r 2 = 0.88, RMSE = 1.86, and MAE = 1.45; however, this result was not as high as the results of Models 5, 6, 7, and 8.
Looking at Table 5, it can be clearly seen that Models 5-8 always include elevation and longitude variables. It would be expected that because Model 9 did not use the longitude variable, therefore the accuracy was not as good as Models 5, 6, 7 and 8.
This indicates that elevation and longitude are the most important variables for T a-max estimation. This result is consistent with adjusted R-squared, BIC criteria and stepwise analysis.

For T a-min Estimation
Looking at Figure 5, it can be clearly seen that the results of all nine models (Models 10-18) were not significantly different. In other words, the accuracy of the model was not increased from the simplest (Model 10, using one variable) to the most complex model (Model 17 with 10 variables).
In comparison to previous studies (see Table 1), we achieved better results (similar r 2 , but smaller RMSE and MAE) of T a-max estimation due to the combination of four LST data and auxiliary variables. However, this combination just made a slight improvement for the T a-min estimation (comparing the result of Model 10, versus . In a further study, a better method for increasing the accuracy of T a-min estimation should be examined. Considering the coefficient of determination (r 2 ), the accuracy (RMSE, MAE) and the number of variables used per model, we would regard Model 5 (for T a-max estimation) and Model 15 (for T a-min estimation) as the best models.

Conclusions
In this study, we have analyzed and discussed the relationship between T a-max , T a-min and four LST products (LST td , LST tn , LST ad , LST an ). The simple method of multiple linear regression analysis was used, and a high accuracy was achieved with r 2 = 0.93, RMSE = 1.43, MAE = 1.08 and r 2 = 0.88, RMSE = 2.08, MAE = 1.60, for T a-max and T a-min , respectively.
When estimating T a using one LST datum solely, T a-min showed a better result than T a-max (Model 1 versus Model 10 in Figures 4 and 5). Multiple linear regressions always give better results than simple linear.
An interesting result is that when we directly compared LST data versus T a , LST nighttime showed a stronger correlation with both T a-max and T a-min than LST daytime; T a-min had a better correlation with LST data than T a-max (see Figure 3a). However, the results of modeling showed that T a-max can be estimated with better results (higher r 2 and lower RMSE, MAE) than T a-min when adding auxiliary variables into the models. It could be concluded that in T a estimation, it is not possible to see the relationship between T a and LST from a directed comparison, because there are other factors that also affect that relationship.
Several model analyses indicate that MODIS LST represents the most important variables for T a estimation. However, to achieve the best results, other variables, such as day length (in hours), Julian day, longitude, latitude and elevation, should be taken into consideration and put into the models.