The study area is the Shiyang River basin in the arid region of Northwest China, which is located in the Hexi Corridor of Gansu Province and the coordinate range is between 101°40′E-104°20′E and 36°30′N-39°30′ N (Figure 1). It covers an area of 41,600 km², and the range of altitude is 1,157m–5012 m. This area is of a continental temperate arid climate, with the characteristics of aridity where precipitation is 300 mm per year and average annual pan evaporation is 2,000 mm.

Meteorological stations in Northwest China are sparse, with only five stations in and near the study area (Minqin, Wuwei, Yongchang, Wushaoling, Jintai station), which can be collected from the China Meteorological Science Data Center (http://data.cma.cn). The observation station (National Field Scientific Observation and Research Station on Efficient Water Use of Oasis Agriculture in Wuwei of Gansu Province) has a meteorological station to collect Ta data, SI-111 thermal infrared radiometers to obtain Ts data (the coordinates in Table 1). Google Earth Pro was used to verify the location of the meteorological station and ensure that Ta data corresponds totally to the site location. Then obtain the available satellite images from 2014 to 2019 and calculate and extract the Ts data corresponding to the location of the meteorological station through the GEE platform.

Google Earth Engine (GEE) is a geospatial processing platform based on cloud computing developed by Google, which promotes fast analysis by using Google’s computing infrastructure and providing a convenient platform for applications based on linking to the cloud computing engine (Becker et al., 2021). Most of the algorithms built into the GEE cloud computing platform use pixel-by-pixel calculation functions, so no matter what the area or proportion of the calculation and analysis is required, as long as the research area has available data. The platform is suitable for scientific researchers with a background in non-professional programming and can quickly realize global-scale remote sensing data processing and mining. This research uses the online JavaScript API of the GEE platform (https://earthengine.google.com/) to access and analyze the data sets used from the public catalog, without downloading images, only outputting the processing results, which improves computing efficiency.

This study uses the Landsat 8 Raw and SR (Surface Reflectance product) dataset in the GEE platform, screening the images with cloudiness not exceeding 30% from 2014 to 2019 and perform cloud removal processing. To cover the entire area of the Shiyang River basin using images of four Landsat-8 tiles 131/033,131/034 and 132/033,132/034 (Figure 1). By using the GEE platform, the remote sensing images were spliced efficiently, clipped, and parameter inversion, also the meteorological data were interpolated.

At present, there are three main methods for remote sensing to retrieve land surface temperature: atmospheric Radiative Transfer Equation (RTE) method, single-channel algorithm, and split-window algorithm. This study uses the Landsat 8 SR data set in the GEE platform, to retrieve the Ts based on the RTE method.

The expression of the radiation transfer equation of the thermal infrared radiation value (L_λ) received by the satellite sensor is (Li et al., 2013; Windahl and Beurs, 2016): L λ = [ ε B ( T s ) + ( 1 − ε ) L ↓ ] τ + L ↑ (1) where L λ is the spectral radiance value at the top of the atmosphere at the band λ ( W ⋅ m − 2 ⋅ μ m − 1 ⋅ s r − 1 ); ε is the surface specific emissivity; B (T_s) is the blackbody thermal emissivity brightness ( W ⋅ m − 2 ⋅ μ m − 1 ⋅ s r − 1 ); τ is the atmospheric thermal infrared band transmittance; L↓ is the downward radiance of the atmosphere after reflection on the ground ( W ⋅ m − 2 ⋅ μ m − 1 ⋅ s r − 1 ); L↑ is the upward radiance of the atmosphere ( W ⋅ m − 2 ⋅ μ m − 1 ⋅ s r − 1 ).

Knowledge of land surface emissivity (LSE) is necessary to apply the above methods to a Landsat image. Considering different situations, obtain the emissivity value from NDVI (Sobrino et al., 2004; Orhan and Yakar, 2016): ε = { 0.986                 NDVI ≤ NDVI S o i l 0.004 [ NDVI − NDVI S o i l NDVI V e g − NDVI S o i l ] + 0.986 NDVI S o i l < N D V I < NDVI V e g 0.990             NDVI ≥ NDVI V e g (2) where NDVI V e g = 0.7 and NDVI S o i l = 0.05 (Ma and Pu, 2020).

When using Landsat8 images, take the 10th band to provide the thermal infrared radiance value. The calculation formula of B ( T s ) is: B ( T s ) = [ L λ − L ↑ − ( 1 − ε ) τ L ↓ ] τ ε (3)

The calculation of the surface temperature (T_s) uses the Planck formula: T s = K 2 ln ( K 1 B ( T s ) + 1 ) − 273.15 (4) where Ts is surface temperature (°C); K 1 , K 2 can be obtained from the header file of the remote sensing data. For Landsat8 TIRS Band10, K 1 = 774.89 W · m − 2 ⋅ μ m − 1 ⋅ s r − 1 ,   K 2 = 1321.08 K .

It can be seen that the use of Radiative Transfer Equation method to retrieve the Ts needs to have the atmospheric profile parameters, which can be obtained by entering the shadowing time, latitude, and longitude in the website provided by NASA (http://atmcorr.gsfc.nasa.gov/).

The near-surface air temperature (Ta) has a good correlation with the surface temperature (Ts) (Benali et al., 2012; Ruiz-Álvarez et al., 2019). Ts is a physical quantity that reflects the degree of cold and heat on the surface of a ground object because it is affected by the characteristics of the underlying surface, such as vegetation coverage and dry and wet conditions. Ta is a physical quantity reflecting the degree of cold and hot air in the atmosphere. The atmosphere has strong fluidity, which is easily affected by the surrounding environment (Xu et al., 2012; Gholamnia et al., 2017). Therefore, when looking for the correlation between Ts and Ta, the influence of various factors such as ground characteristic and environment must be considered. For the estimation of Ta, Jang et al. (2004) showed that Julian day is a more significant parameter than altitude or the solar zenith angle. In addition, we have chosen variables that are often selected as predictors in air temperature modeling literature, such as surface albedo, elevation (DEM), normalized vegetation index (NDVI), and wind speed (Riddering and Queen, 2006; Cristóbal et al., 2008; Hou et al., 2013). The full name of the DEM (WGS84/EGM96) data source is “NASA SRTM Digital Elevation 30 m,” which is used directly on the GEE platform and provided by NASA JPL (Farr et al., 2007).

The broadband albedo ( α ) of the ground surface is a critical variable for many scientific applications, which is the ratio of the total radiant flux reflected by the ground surface to the incident flux (Liang et al., 2003). The calculation formula of α is: α = 0.356 α 2 + 0.130 α 4 + 0.373 α 5 + 0.085 α 6 + 0.072 α 7 − 0.0018 (5) Where α is the surface reflectance, and its value is between 0–1.0; α₂, α₄, α₅, α₆, and α₇ are the 2, 4, 5, 6, and 7 bands of Landsat 8 surface reflectance products.

As has been commented in the previous section, the LSE can be retrieved from NDVI values. The data can be used to construct NDVI according to the following equation: N D V I = N I R − R N I R + R (6) Where NIR and R are the reflection values of the near-infrared and infrared bands, respectively, which are the fifth and fourth bands of Landsat 8.

This study mainly considers the relationship between Ta and Ts, NDVI, surface albedo, DEM, wind speed and Julian day. The flow chart of the models used is shown in Figure 2, and the process can be summarized as Data preparation, processing and prediction. Data collection is the basis for establishing a reliable model, but it is also necessary to filter the input variables and try to use as few variables as possible under the condition of high model accuracy. Variables used in this study are readily available, which are closely related to the changes of Ta. Previous studies on Ta estimation lack the verification of Ts and analysis of related results. Therefore, this research explored whether it is necessary to improve the precision of Ts to achieve better results in actual application. After verification and comparison, the best model is selected from the combined methods to estimate Ta.

The spatial and temporal matching among the variables is the main issue for the reliability of the regression implementation. Table 1 shows the sources and resolutions of all data sets. On the spatial scale, the resolutions of spatial variability independent variables such as Ts, Albedo, NDVI, elevation, wind speed (resampled after interpolation) are all 30 m × 30 m, which is spatially consistent. On the time scale, Ts, Albedo, and NDVI are instantaneous data, and only wind speed is daily data.

SI-111 collected Ts data for the 5 years from 2014 to 2018. Based on Landsat8, using RTE to invert the Ts, 40 points corresponding to the position of the ground instrument were extracted. The RTE method estimates Ts from 2014 to 2017, compared with observation data and its R ² = 0.77, RMSE = 3.27°C. This study found that the RTE method overestimated the value of Ts. Especially when the surface vegetation coverage is low, the error will be more obvious (Figure 3). To resolve this uncertainty, a correction method using NDVI is proposed to improve the accuracy of the RTE method for inversion of Ts.

Correcting RTE with NDVI, the expression is: Ts R T E − N D V I = Ts R T E − 0.22 NDVI − 1.5 (11) where NDVI is Normalized Difference Vegetation Index; Ts R T E − N D V I is the surface temperature calculated using the Radiation Transfer Equation corrected by NDVI (°C);  Ts R T E is the surface temperature calculated using the original Radiation Transfer Equation (°C).

The 31 measured data from 2014 to 2017 were used as the target to revise the RTE (Figure 4A), and the 2018 data were used as verification (Figure 4B). After correction using NDVI of 2014–2017, compared with the verified data, the accuracy of the model is R ² = 0.83, RMSE = 2.09°C. And after the verification of the data in 2018, it is confirmed that Eq. 11 has improved the accuracy of RTE (before: R ² = 0.63, RMSE = 4.80°C; After: R ² = 0.76, RMSE = 2.41°C). The fitting line is close to 1:1, indicating that the estimation accuracy of the Ts is significantly improved.

Seven variables were tested for multicollinearity, and there was no variable with VIF greater than 10, indicating that the regression model can be established. As shown in Figure 5, RSS monotonically decreases with the increase of the number of independent variables, which cannot be directly judged. Adjusted R ² and Cp values tended to be the largest when the number of variables was 6 and 7. The adjusted R ² for Ta models ranged from 0.87 to 0.95, which eventually fixed at 0.95. Further addition of the independent variable did not improve the adjusted R ², indicating that 6 independent variables were considered optimal. The change of BIC with the number of variables also proves that the model with six variables has the highest accuracy.

The statistical analysis indicated that Ts, Albedo, NDVI, DEM, x7, and x8 were significant independent variables for estimating Ta. As shown in Figure 6, R ² was increased with the increase of the independent variables, so optimal combination of independent variables was not proposed in this study. However, in the absence of the wind speed, the adjusted R ² reached a maximum while Cp and BIC arrived at a minimum, indicating that wind speed does not improve the accuracy of the Ta model. Except that the wind speed does not affect the model construction, the addition of other variables can improve the fitting accuracy of the model. The date of Ta estimation is almost always sunny when the average wind speed is slightly different in space. Therefore, there is no significant correlation between wind speed and Ta variation. Moreover, the average wind speed data is interpolated from the data obtained from meteorological stations, which may be inaccurate at the regional scale. To sum up, it is sufficient to use remote sensing data, elevation, and Julian days as the independent variables for the model.

After determining the independent variables, the formula of the constructed multiple linear regression equation is as follows: T a = b 0 + b 1 T s + b 2 α + b 3 N D V I + b 4 H + b 5 x 7 + b 6 x 8 (12) where Ta is the near-surface air temperature; Ts is the surface temperature; a is the surface albedo; NDVI is Normalized Difference Vegetation Index; H is the altitude; x 7 = ( J − 200 ) 2 ;   x 8   =   x 7 × T s ; J is Julian Day (the number of the day of the year); b_0∼b₆ are undetermined coefficients.

The random forest algorithm cannot provide an estimation form similar to MLR because it cannot be parameterized. Two parameters are critical to the operation of the random forest model: the number of trees to grow (ntree) and the number of variables to be used at each node (mtry). These two parameters are selected based on the percentage of model interpretation and mean of squared residuals. The mtry value was tested from 1 to 8, and the ntree value was tested using the following 6 values: 100, 300, 500, 1,000, 1,500, and 2,000.

Figure 7 shows the changes in Mean of squared residuals and Percent variance explained after running a random forest with different ntree and mtry values. The highest Percent variance explained and the lowest Mean of squared residuals were obtained when mtry = 3 and ntree = 1,500 or 2,000. Since the higher the value of ntree, the more cost and time for calculation, so the value of ntree is set to 1,500, and the value of mtry is set to 3 to run the random forest model.

The amount of Ta-related data from 2014 to 2019 is 397. Randomly select 70% of the total amount as the training set to determine the variables and parameters, and the remaining data as the testing set. Figure 8 is a scatter plot of the simulated and observed values of the two models constructed from the training set data. The R ² and RMSE of the MLR model were 0.953 and 1.74°C, respectively. Most of the points of the RF are concentrated and close to the 1:1 line. The comparison between the results of the RF model and the MLR model shows that the random forest has higher estimation accuracy.

To understand the universality of the built model, the MLR model was verified using the 2018–2019 dataset; the RF algorithm was verified using the reserved test dataset. At the same time, it is explored whether the improvement of the estimation accuracy of Ts will affect the estimation accuracy of Ta.

The prediction result of Ta estimated through the validation data set (Figure 9) cannot reach the accuracy of the initial training set. When Ta is estimated from the validation data set of the RF, the prediction could not achieve the accuracy of the initial training set, respectively. When using NDVI to correct RTE, the R ² and RMSE of the RF model was 0.943 and 2.12°C, respectively, indicating that the estimation accuracy of using RTE-NDVI is reduced. In contrast, the simulation results of MLR on Ta are relatively stable. The two calculation methods of Ts have little effect on MLR. Therefore, the accuracy of Ts inversion will not affect the accuracy of the two methods for estimating Ta. The estimation results of all models show good performance, and the fitting line is close to 1:1.

Through the analysis of the average daily Ta and error estimated by MLR for 2018–2019, it was found that the temperature and error on the time scale have spatial variability. As shown in Figure 10, Ta had a clear trend with Julian day and season, and the estimation accuracy had an apparent seasonal trend. Sometimes, the meteorological station is covered by clouds and cannot participate in the calculation or verification, especially at the high-altitude Wushaoling station. The amount of remote sensing image data in winter is insufficient. Although the error is small, it cannot be fully verified. There was not much difference in the accuracy between spring and summer. The RMSE values in spring and summer were both 1.66°C, which had no difference in the estimation accuracy. With the increase of temperature, the error value was relatively stable and had no apparent change trend. In autumn, the RMSE was 2.35°C. At this time, the study area is in the rainy season, and the estimate of Ta will be affected by the increase in rainfall.

In 2018 and 2019, the MEF value of each station was between 0.84 and 0.97, and the overall MEF was 0.946, indicating that MLR has good performance. Analyze the error distribution of Ta estimated by the MLR model at each site. Figure 9 shows that Ta is more accurate at locations far from the mountainous area.

The MEF of the Minqin station for 2 years was close to 1, and the degree of the fitting was relatively high. The MLR model behaved differently in different altitude ranges. As the altitude increases, the performance of the MLR model gradually weakens, which may be caused by the complexity of the terrain. The RMSE has the maximum value at the highest altitude site (Wu Shaolin) (Figures 11C,D). The average RMSE of Minqin and Jintai Station are below 1.5°C. The accuracy of these two stations is better than that of other stations, possibly because the terrain and environmental conditions of the stations are not complicated. The locations of these two stations are far away from the Qilian Mountains, which are as high as 5,000 m above sea level. The estimation model had the highest relevance and accuracy at the Minqin station, which is located on the north side of the basin, and the surrounding terrain is flat and less affected by high mountains. Although Jintai Station is close to the east of Qilian Mountain, which is far from the mountain range, the estimation accuracy of Ta is also very high.

Because the data of the two adjacent columns of the Landsat satellite are of different periods, remote sensing images of the entire study area cannot be obtained on the same day, so the basin can only be divided into two parts for comparison. Chosen four images with fewer clouds, then calculated the difference between the spatial interpolation result and the MLR model estimate. Inverse Distance Weighted (IDW) assumes that the influence of variables on the surrounding area decreases as the distance to the sample location increases, which was used to interpolate Ta data from the meteorological station. The error of Ta distribution obtained by IDW in the area is related to the underlying surface. Using meteorological and remote sensing data from July 31, August 22, 2018, May 21, May 28, 2019, Ta spatial interpolation and estimation maps were obtained and compared.

It iswidely believed that Ta is related to altitude. Temperature decreased with the increase of altitude in most cases. In high altitude areas, where no meteorological station provides temperature data, the Ta of the regions was determined by data from nearby stations. Therefore, Ta in high altitude areas was often overestimated by interpolation (the red area below Figures 12A,C). Simultaneously, there are also large desert areas in the study area (the blue area in Figure 12). Since the temperature in the desert during the day will be much higher than in other zones, the interpolation method can underestimate the temperature by up to 15°C. The Ta of farmland is overestimated during spatial interpolation (the part of the red area above d in Figure 12B), which is also related to the location of meteorological stations. Most of the station is established in cities or border areas. The vegetation coverage of farmland is very high during the growth period of crops, which will affect its Ta. Therefore, when using interpolation methods such as IDW, more meteorological stations are needed to ensure the accuracy of the interpolation results, and at the same time, should pay attention to the unreliability of interpolation methods in mountains or deserts.

In areas with vegetation cover, the temperature value retrieved from remote sensing data is often a mixed value of Ts and canopy temperature (Tc). The accurate estimation of Ts still needs to be analyzed according to the vegetation situation and terrain (Zhang and Li, 2018). Comparing the Ts calculated by RTE with the measured values, it was found that RTE overestimated the Ts, especially when the vegetation coverage was low, consistent with the conclusions published by USGS in 2015 (Barsi et al., 2014; Xu and Huang, 2016). However, due to the characteristics of fewer parameters and suitable for any thermal infrared band, the improved accuracy can also be widely used. Ts is considered to be the most important independent variable the many estimated Ta models (Park, 2011; Zhang et al., 2016; Song and Wu, 2018). In this study, the RMSE of Ts directly calculated by RTE is 3–5°C, and it is 2–3°C after correction. The correction equation is proposed based on the local estimation results, and its universality still needs to be verified. The error is related to the surface emissivity or atmospheric profile parameters (Li et al., 2013; Sekertekin and Bonafoni, 2020).

The variables of this study to be selected are closely related to Ta in many studies (Cristóbal et al., 2008; Gholamnia et al., 2017; Shen et al., 2020). Although the correlation between time parameters and Ta was not strong, the adjusted x7 had a significant influence on model accuracy indicating that the adjustment using the Julian Day was useful. Studies have suggested that wind speed is significant, especially when using energy balance methods, which require wind speed to calculate aerodynamic resistance (Hou et al., 2013; Zhang et al., 2015). Through indicators comparison, it is found that when using a statistical method to build the model, the average wind speed is not helpful to improve Ta estimation accuracy (Stisen et al., 2007).

The selected variables participate in MLR and RF to estimate Ta, using the available satellite dataset. The RMSE of the two methods were both lower than 2.0°C, which indicated that they were both suitable for northwest China. Previously, Ho et al. (2014) Using the RF algorithm to estimate the maximum Ta in Vancouver, Canada, and the RMSE was 2.3°C. The optimal MLR model established by Yang et al. (2017) estimates the average Ta and the RMSE is 3.6°C. Therefore, the methods of estimating Ta in this study had relatively high accuracy. However, due to the 16-days revisit cycle of the Landsat satellite and the low temporal resolution of the images, there are limitations in monitoring daily Ta (Yoo et al., 2018).

When performing regression, the random forest cannot make predictions that exceed the range of the training dataset, which may lead to overfitting when modeling some data with specific noise. For many statistical modelers, the random forest is like a black box (Zheng et al., 2019), almost unable to control the internal operation of the model, can only try between different parameters and random allocation. The simulation results of the random forest algorithm for the existing data are very good, but when it is used in the prediction and estimation, its accuracy may drop suddenly. It is difficult to improve it because it cannot control the internal operation of the model. For this study, multiple linear regression is still recommended. The simulation accuracy of the MLR model is similar to that of validation, which is conducive to the evaluation of subsequent prediction results. From the comparison results of the models, the MLR and RF models that without improving the accuracy of Ts performed better, and the actual operation was simpler. The reason for this phenomenon may be that the NDVI involved in the Ts correction is also one of the independent variables of the Ta estimation model, which makes the internal adjustment of the model during operation can ignore the errors of Ts.

The seasonal variation of the average temperature and the distribution of the station will affect the estimation accuracy of Ta (Holden et al., 2011; Chen et al., 2015). From the time point of view, the RMSE in autumn reaches 2.24°C, which is 0.57°C higher than that in spring and summer. This is consistent with the results of Golkar et al. (2018), which believe that the estimation of Ta in spring and summer has higher accuracy, while there are more uncertainties in autumn. Therefore, the estimation model of daily average temperature is more suitable for spring and summer days. Yang et al. (2017) integrated various statistical indicators and believed that each model performed better in spring, and the estimation accuracy decreased due to the influence of rainfall and cloudy weather. The study area is located in the inland of Northwest China, and rainfall is mainly concentrated in autumn. Benali et al. (2012) believe that the cloud cover of remote sensing images is inversely proportional to the model performance, and higher cloud cover harms model performance. That is, because cloud cover could reduce the accuracy of the Ts, which affects the accuracy of estimated values of daily average temperature. Therefore, the model can be optimized to the season or month scale, and the process of climate and environmental impacts can be added to adapt to the impact of climate change and cloud cover.

The spatial distribution map of RMSE shows that the model in mountainous and plateau areas generally performs worse than in plain areas. This conclusion is consistent with the research results of Yang et al. (2017). Due to the influence of terrain differences, mountainous terrain has a broader spatial variation range than flat terrain, so mountain temperature changes are more complicated. On relatively flat terrain and vegetation, horizontal uniformity leads to stable atmospheric conditions (Blandford et al., 2008; Lin et al., 2016). Ta raster data is usually limited by site coverage, especially in mountainous areas where the density and elevation distribution of meteorological observations vary greatly (Holden et al., 2011). Shen et al. (2020) believe that in the case of large-scale use of MLR and RF, the accuracy of the estimation results is poor due to the scarcity of sites in Northwest China that can participate in model training. When MLR is estimated in the Shiyang River Basin, the MEF of the highest altitude station (3000 m) is higher than 0.83, and the RMSE = 2.10°C, the model also performs well.

The method proposed by the research can obtain high-precision Ta estimation results when using the Landsat 8 data set, and improve the spatial resolution. However, the actual situation of only six meteorological stations limits the accuracy of Ta inversion in this area. We hope that more reliable verification data can be obtained in future studies.