Estimation of 30 m land surface temperatures over the entire Tibetan Plateau based on Landsat-7 ETM+ data and machine learning methods

Land surface temperature (LST) is an important parameter in land surface processes. Improving the accuracy of LST retrieval over the entire Tibetan Plateau (TP) using satellite images with high spatial resolution is an important and essential issue for studies of climate change on the TP. In this study, a random forest regression (RFR) model based on di ﬀ erent land cover types and an improved generalized single-channel (SC) algorithm based on linear regression (LR) were proposed. Plateau-scale LST products with a 30 m spatial resolution from 2006 to 2017 were derived by 109,978 Landsat 7 Enhanced Thematic Mapper Plus images and the application of the Google Earth Engine. Validation between LST results obtained from di ﬀ erent algorithms and in situ measurements from Tibetan observation and research platform showed that the root mean square errors of the LST results retrieved by the RFR and LR models were 1.890 and 2.767 K, respectively, which were smaller than that of the MODIS product (3.625 K) and the original SC method (5.836 K)


Introduction
The Tibetan Plateau (TP) is known as the highest plateau in the world, with an average elevation of over 4,000 m. The thermal and dynamical effects of the TP have important impacts on the Asian monsoon system and atmospheric circulation of Asia and even the globe (Yao et al. 2012;. Additionally, the TP has the most widely distributed elevational permafrost at mid-latitudes and globally due to its high average altitude. Land surface temperature (LST) plays an important role in the energy exchange and water cycle processes between the land surface and atmosphere, and it is critical for the evaluation of climate change impacts on the ecosphere, hydrosphere and particularly permafrost in the cryosphere (Becker and Li 1995;Hachem, Duguay, and Allard 2012;Li et al. 2013;Luo et al. 2018;Sekertekin and Bonafoni 2020;Wan et al. 2004;. Many previous studies have indicated that the TP has been warming rapidly over the past decades. The linear rates of air temperature increase in the permafrost region over the TP are approximately three times the global warming rate (Harris 2010;Niu et al. 1978;Qiu 2008). It is critical to monitor air temperature changes over the permafrost area. However, air temperature data are either not spatially representative (from meteorological stations) or too coarse (reanalysis) for regional warming studies. Different from the ground surface temperature (GST) measured at 0-5 cm under the land surface, LST has been found to be better correlated with air temperature 1-3 m above the ground. High-resolution LST also offers enormous potential for monitoring air temperature changes over the TP.
In situ measurements, satellite remote sensing and numerical modelling are the three primary methods used to derive the land surface characteristic parameters (Liu and Chen 2000;Luo, Jin, and Bense 2019). Generally, the results from in situ measurements are considered to be the most reliable. However, due to the complex underlying surface status and harsh environment of the TP, meteorological stations are both sparse and unevenly distributed (Song et al. 2014). The lack of in situ measurements across large areas of the TP means that this data source alone cannot meet the demands of land-atmosphere interaction studies. Satellite remote sensing provides a more efficient way to retrieve LST data at larger spatial scales and longer temporal scales (Ma et al. 2003;Oku and Ishikawa 2004). LST datasets over the entire TP and even globally have been obtained in previous studies (Ma et al. 2011;Salama et al. 2012;Wan 2008;). However, these LST datasets were primarily obtained from satellite data with a low to medium spatial resolution. Unfortunately, results with coarse spatial resolutions are insufficient when investigating the details of heterogeneous pixels (Liu, Hiyama, and Yamaguchi 2006).
The Landsat series (Landsat 5, 7, 8 and 9) have the potential to provide LST information with a high spatial resolution. A number of LST retrieval algorithms for Landsat series have been proposed in previous studies Jiménez-Muñoz et al. 2009;Jiménez-Muñoz and Sobrino 2003;Parastatidis et al. 2017;Qin, Karnieli, and Berliner 2001;Vlassova et al. 2014). However, most previous studies on LST retrievals with a high spatial resolution only focused on the mesoscale due to limitations in both computing capability and data size (Ge et al. 2019;Ma et al. 2006;Zhao et al. 2020). LST data with a fine spatial resolution over the entire TP are particularly lacking. The TP has a higher altitude with different underlying surface types from other plain areas. Because in situ measurements in the TP have been seldom used in retrieval algorithm development, those algorithms based on the measurements in the plain area are not suitable for LST retrieval in the TP. Therefore, we propose two algorithms for retrieving the plateau LST from Landsat 7 Enhanced Thematic Mapper Plus (ETM+) data. The generalized singlechannel (SC) algorithm was developed using a linear regression (LR) model and 10-year in situ measurements over the TP. LST datasets covering 12 years with a 30 m spatial resolution for the entire TP area were obtained with the optimized algorithm. The retrieved results can describe a finer spatiotemporal distribution of the LST over the entire TP. Additionally, a new algorithm to estimate LST based on a random forest (RF) model was also proposed in this study. Using the cloud computing platform Google Earth Engine (GEE), the 12-year LST results at 30 m resolution were obtained efficiently within 9 h.
The goals of this study are (1) to develop efficient algorithms that can more effectively determine the LST over the TP using machine learning models; and (2) to construct a time series of land surface temperature with a fine spatial resolution for the entire TP. In this study, LST with the normalized difference vegetation index (NDVI) and land surface emissivity (LSE) with a 30 m spatial resolution over the entire TP were retrieved based on a total of 109,978 Landsat ETM+ images from 2006 to 2017. LST data of approximately 2 × 10 11 pixels at a spatial resolution of 30 m between 20°N-50°N and 60°E-110°E were processed within 9 h using the improved SC algorithm with the GEE platform. The data and methodology are described in Section 2. The retrieved LST results were validated with in situ measurements, as described in Section 3. A discussion of various topics is presented in Section 4, and concluding remarks are given in Section 5.

Satellite data
Landsat 7 was launched in 1999 and carries an enhanced thematic mapper plus (ETM+) sensor. ETM+ has six reflective bands and one thermal infrared band (Sekertekin and Bonafoni 2020) (Table 1). The Landsat 7 ETM+ data used in this study were provided by the United States Geological Survey (USGS) and are available in GEE, which provides access to multiple ETM+ datasets, including raw scenes (level-1 product), surface reflectance (SR) data and top of atmosphere reflectance (TOA) data. GEE is an open-source platform provided by Google to a wide range of professional scientists who focus on atmospheric science, environmental science and geophysics to process massive amounts of data. It provides a programming front-end for users to write code in JavaScript and Python (code.google.earthengine.com), and provides a variety of interfaces for satellite remote sensing data, including MODIS products, Landsat products, and NCAR products. (Gorelick et al. 2017). The USGS also provides Level-2 LST products retrieved with a single-channel (SC) algorithm that was jointly developed by the Rochester Institute of Technology (RIT) and National Aeronautics and Space Administration (NASA) Jet Propulsion Laboratory (JPL). In addition, the MODIS/Terra Land Surface Temperature/ Emissivity Daily L3 Global LST product (MOD11A1) is one of the most commonly used daily LST products and is produced using a generalized split-window algorithm. In addition to the retrieved LST in this study, the USGS Collection-2 Level-2 Tier-1 LST products of Landsat 7 and 8 and MOD11A1 were validated with in situ measurements for comparison.

Reanalysis data and in situ measurements
As an important variable in the LST retrieval algorithm, atmospheric total column water vapour (TCWV) is required in addition to satellite radiation data Jiménez-Muñoz and Sobrino 2003;Jiménez-Muñoz and Sobrino 2010;Sekertekin and Bonafoni 2020). GEE only provides water vapour reanalysis data from the National Centres for Environmental Prediction and the National Centre for Atmospheric Research (NCEP/NCAR), and its 6-hour temporal resolution is too low. Datasets from the European Centre for Medium-Range Weather Forecasts Reanalysis 5 (ERA5) can provide TCWV data with hourly temporal resolution that are more consistent with satellite transit times. The ERA5 water vapour data were validated using atmospheric radiosonde data from the Japan International Cooperation Agency (JICA) China-Japan climate disaster mitigation project (JICA/Tibet project) over the TP. There are six radiosonde sites used to provide the measured TCWV data, including Changdu, Dingri, Gaize, Lhasa, Linzhi, and Naqu in Section 3.1. The RMSE between the sounding data and ERA5 is calculated by Equation (1): where ω JICA is the TCWV from the JICA sounding experimental data and ω ERA5 is the TCWV from the ERA5 reanalysis data. The spatial resolution of the ERA5 data used in this study is 0.5°× 0.5°. Although the water vapour data of ERA5 are not provided by GEE currently, GEE provides a convenient way for users to upload grid data. The upload interface provided by GEE was used in this study to upload 12 years (2006-2017) of gridded ERA5 water vapour reanalysis data.
The in situ LST used in this study is obtained from measurements at NADORS, NAMORS, MAWORS, BJ, SETORS, QOMS, NPAM (MS3478) and D105 on the TP (Figure 1) in the Tibetan observation and research platform (TORP) (Ma et al. , 2014. The temporal resolution of these in situ measurements is one hour, and the basic information of these stations is shown in Table 2. The in situ land surface temperature, LST in−situ , was directly computed as: where LWU and LWD are the measured upwards and downwards longwave radiances, respectively; e b is the broadband LSE (BBE); and s is the Stefan-Boltzmann constant, which equals 5.67 × 10 −8 W · m −2 · K −4 . The BBE was calculated by Equation (3) (Cheng et al. 2017): where r i is the surface reflectance of ETM+ band i and e b is the BBE. The LST data from in situ measurements were linearly interpolated to match the satellite transit times. There were 1,592 LST records from TORP that were divided into two parts: the training set consisted of 1,273 records, which accounted for 80% of the total records and was used to train the model; the remaining 319 LST records, which accounted for 20% of the total records, were used to test the performance of the RF model. These two datasets were completely randomly divided.

Cloud identification and gap filling of Landsat SLC-off images
The accuracy of calculating LST from optical remote sensing data is often influenced by clouds; it is thus typically necessary to mask cloudy pixels in satellite images. Cloud identification was based on a quality assessment band (BQA) to ensure that pixels with cloud cover greater than 0 were masked out. The SLC of Landsat 7 failed at the end of May 2003, and the images acquired by the Landsat 7 ETM+ after that time had a problem that manifested as missing stripes (Markham et al. 2004). To fill these image gaps, a gap-filling method depending on adjacent nonempty linear interpolation was developed in this study ( Figure 2). For images with gaps, an image containing the median value of each pixel around the same days between the previous year and the following year was obtained. Then, this image was used to establish a linear relationship between the missing pixel and surrounding 21*21 no-null pixels. Values of missing stripes can be obtained from this linear relationship. This gap-filling method was applied to the red band (B3), near infrared band (B4) and thermal infrared band (B6) of ETM+. As shown in Figure 3, after applying the gap-filling algorithm, the problem of missing stripes was effectively resolved. In addition, the masked pixels with cloud cover greater than 0 were also compensated for via this interpolation method to a certain extent. The spatial distributions of gaps filled with NDVIs (Figure 3(b)) agree with the general spatial NDVI distributions shown in Figure 3(a). Figure 3(b) also shows that for pixels not located in the gaps, the NDVIs maintained their original values. Therefore, the missing data were effectively addressed after applying the gap-filling algorithm. The performance of the gap-filling method is shown in the discussion section.

Estimation of LST
The SC algorithm has been widely used to retrieve LST based on different satellite sensors, such as Thematic Mapper (TM)  Guo, and Wu 2014). The process of retrieving LST from the raw data by the SC algorithm is shown in Figure 4. The use of the generalized SC method requires high-quality atmospheric transmittance data to estimate the key parameters in the radiative transfer equation. Equation (4) is the primary formula for the SC algorithm.
where LST SC is the LST obtained from the SC algorithm; e is the LSE; C 1 , C 2 and C 3 are three parameters called atmospheric functions (AFs), which are related to the atmospheric water vapour content; L S is the radiation (in W · m −2 · sr −1 · mm −1 ) measured by the satellite sensor; and d and g are parameters that depend on Planck's function.
The formulae for calculating d and g are as follows: where c 1 and c 2 are Planck's constants (c 1 = 1.19104 × 10 8 W · mm 4 · m −2 · sr −1 , c 2 = 1.43877 × 10 4 mm · K); l is the wavelength of the TIR band, which equals 11.45 mm for the Landsat 7 ETM+ sensor; and T b is the brightness temperature (in units of K) of the satellite sensor. The function that defines the relationship between T b and L S is as follows: where K 1 and K 2 are the calibration constants of the satellite TIR band and are considered to be 666.09 and 1282.71 for the Landsat 7 TIR band, respectively (Chander, Markham, and Helder 2009;Coll et al. 2010). The AFs in (1) are related to the atmospheric water vapour content v. They can be calculated from atmospheric total column water vapour by Equation (8) (Jiménez-Muñoz and Sobrino 2003;Sobrino, Jiménez-Muñoz, and Paolini 2004): a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 The atmospheric total column water vapour (TCWV, in units of g · cm −2 ), v, in Equation (8) is the only external variable of the SC algorithm in addition to the satellite radiation. As a key variable in the SC algorithm, water vapour has important effects on TIR observations, which shows atmospheric absorption to TIR radiation. v can be obtained from the water vapour datasets.
The coefficients a ij in Equation (8) can be determined by simulations and vary according to different TCWV datasets (Jiménez-Muñoz and Sobrino 2003). In this study, the ERA5 reanalysis was used to obtain hourly water vapour data for the entire TP. However, there are currently no suitable AF coefficients a ij for ERA5 data. In this study, AFs were developed by the LR model, and new coefficients with higher precision were proposed. We input the AFs into Equation (4) to obtain the quadratic relationship and fitted it with in situ LST by linear regression to obtain these coefficients (Table 3)

Estimation of LSE
LSE is an important land surface parameter that represents the ability of the land surface to transform heat into radiation (Peres and Dacamara 2005;Sobrino et al. 2008). Several methods in previous studies have been proposed to retrieve the LSE (Sobrino et al. 2008;Valor 1996), including the classification-based emissivity method (CBEM) (Peres and Dacamara 2005), NDVI-based emissivity method (NBEM) (Valor 1996) and two temperature methods (Watson 1992). An NDVI-based method called NDVI THM was used to calculate the LSE (e) in this study. The land surface is considered to be a mix of bare soil and complete vegetation cover in this model and is divided into three cases (Sobrino et al. 2008): (1) the pixels with NDVI < 0.2 are considered as soil with sparse vegetation and even bare soil; (2) the pixels with 0.5 ≥ NDVI≥0.2 are considered as a mixture of bare soil and vegetation; and (3) the pixels with NDVI > 0.5 are considered as fully vegetated. For each pixel in the ETM+ imagery, the LSE values are: where r Red is the reflectance of the red band and P V is the vegetation coverage, which can be calculated by Equation (10): where NDVI S is the NDVI value of barren land and NDVI V is the NDVI value of vegetation, which are set to 0.2 and 0.5, respectively.

Random forest
The random forest (RF) algorithm proposed by Breiman in 2001 is a nonlinear ensemble machine learning (ML) method that was developed based on the bagging algorithm. RF is a machine learning method that has been widely used in hydrological assessment and meteorological forecasting in recent years. RF is composed of multiple decision trees, and the final output of the model is jointly determined by each decision tree in the forest. With regression problems, the mean of the output of each decision tree is the final result. RF can effectively model nonlinear relationships between predictors X i and response variables Y. In this study, RF was considered a regression model (random forest regression, RFR) to estimate LST. To obtain the spatial distribution of ETM+ LST results by training the in situ LST at a point scale, we proposed a three-step methodology: (1) generation of the RF model; (2) validation of the models obtained to ensure the required precision; and (3) estimation of the LST spatial distribution using satellite and reanalysis data. RFR models the relationship between remote sensing indices and LST simulation by a set of decision rules, and assesses the contribution of individual predictors to the model automatically and outputs the importance score. However, the selection of predictors X i must refer to physical correlations with LST retrieval, such as water vapour and brightness temperature (Hutengs and Vohland 2016). In this study, the variables that describe the major drivers of LST are considered predictors X i in the RFR model (Table 4). In previous studies on LST estimation by RFR, types of land cover (LCT) were seldom regarded as predictors. Therefore, in this study, we distinguished different RFR models according to the underlying cover types to characterize the LST estimation process of different underlying surfaces. The LCT data were obtained from Finer Resolution Observation and Monitoring Global Land Cover (FROM-GLC) (Gong et al. 2013). The spatial resolution of the LCT data is 30 m, which is consistent with the spatial resolution of ETM+ data. For each type of land surface, model-trained LST results are obtained from input variables as follows: where subscript i indicates different land cover types.

Statistical indices for ground validation
For the assessment of the performances of the RF model, three statistical indices, including RMSE (root-mean-square error), rRMSE (relative root-mean-square error) and R 2 (coefficient of determination), were applied in this study: where N is the amount of data; x i and y i are the observed LST and LST values computed through different models and LST products, respectively; and x and y are the averages of x i and y i , respectively. In terms of these metrics, the model is denoted as a perfect fit when R 2 = 1, RMSE and rRMSE = 0.

Validation of ERA5 Using JICA sounding data
The ERA5 water vapour data were validated using atmospheric radiosonde experimental data from the Changdu, Dingri, Gaize, Lhasa, Linzhi and Naqu stations in the JICA/Tibet project. As shown in Figure 5, the measured values and ERA5 data were in good agreement with each other. The RMSE between v JICA and v ERA5 is 0.328 g · cm −2 , and the determination coefficient (R 2 ) is 0.880.  images for QOMS, 18 images for SETORS, and 16 images for NADORS. In this study, we used in situ measurements from eight stations on the TP over eleven years (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) to validate the retrieved LST results. The RMSE, rRMSE and R 2 were calculated by Equations (12a)-(12c), respectively. A model with smaller RMSE, rRMSE and higher R 2 can be considered more reliable. As shown in Figure 6(a), the validation set of Landsat 7 ETM+ RFR-based LST (red triangles) was in good agreement with the in situ measurements with an RMSE of 1.890 K and an rRMSE of 0.164. The validation set of LR-based LST (blue points) has an RMSE of 2.767 K and an rRMSE of 0.241, which is less than that of the LST results obtained by the original SC algorithm (green crosses) on the same days with an RMSE of 5.836 K and an rRMSE of 0.501, respectively. The accuracy of LST retrieved from RFR and LR was superior to LST results from the original SC algorithm, which also indicates that the original coefficients of AFs proposed by Jiménez-Muñoz and Sobrino in 2003 are not suitable for the TP area. In addition, the original SC algorithm overestimated higher temperatures, which was improved by the RFR and LR models. Furthermore, we validated the LST results from the MOD11A1 and USGS Level-2 products on the same days as ETM+ data. The spatial resolution of the Landsat 7 ETM+ retrieval results was 30 m, while that of the MODIS LST product was 1000 m. In this study, LST retrieved from ETM+ was upscaled to 1 km to match the spatial resolution of the MODIS LST product. As shown in Figure 6(b), the RMSEs (rRMSEs) between the MOD11A1, USGS LST product and in situ measurements were 3.625 K (0.315) and 4.493 K (0.386), respectively.

Retrieval
To validate the portability of the model, we applied the same RFR and LR models to the B10 band of the Landsat 8 Thermal Infrared Sensor (TIRS). As shown in Figure 6(c), due to the different wavelengths of TIRS B10 (10.60-11.19 µm) and ETM+ B6 (10.40-12.50 µm), the RMSEs between the LST results and in situ measurements increased to 3.415 and 3.011 K for the LR and RF models, respectively. However, they are still lower than those of the USGS and MOD11A1 LST products, as shown in Figure 6(d). Additionally, when compared with the USGS official LST product, the rRMSE Figure 5. Validation of ERA5 TCWV data against JICA densified radiosonde data over the TP. The radiosonde experiment was carried out at 01:00, 07:00, 13:00 and 19:00 (BST) at Changdu, Dingri, Gaize, Lhasa, Linzhi, and Naqu from February to July in 2008. It should be noted that the observation time is not continuous.
of RFR-based and LR-based LST decreased from 0.407 to 0.290 and 0.256, respectively. These assessments indicate that the RFR model and the SC algorithm improved by LR can reduce the biases in the retrieval of LST on the TP. The performance of machine learning is critical for improving the accuracy of LST retrieval.
Because the spatial resolutions of the MODIS LST product (MOD11A1, with a 1,000 m spatial resolution) and the LST retrieved from the Landsat 7 ETM+ data (30 m spatial resolution) are different, one pixel of the MODIS LST may contain multiple pixels of ETM+ data. As shown in Figure 7, the ETM+ results have higher resolution and are more consistent with the land cover types. The LST results of ETM+ (Figure 7(b)) show that the LST of the grassland around the SETORS station is higher than that of the forest area. These cannot be seen from the MODIS LST with a coarser resolution. The LST results that were obtained by the improved algorithm and ETM+ data can thus describe the LST distribution more precisely. An LST with high spatial resolution can provide more accurate data for studies of land-atmosphere interactions and the energy and water cycles on the TP. LST from the USGS Landsat 7 Level 2 product (blue points) and MOD11A1 V6 daily LST product (red triangles) on the same day as ETM+ results, and (c) LST retrieved from Landsat 8 TIRS Band 10. The retrieved LSTs were obtained from two methods: RFR (red triangles) and LR (blue points), (d) LST from the USGS Landsat 8 Level 2 product (blue points) and MOD11A1 V6 daily LST product (red triangles) on the same day as the TIRS results against in situ measurements. LST results from ETM+ and TIRS have been upscaled to a 1 km spatial resolution.

Spatiotemporal distribution of land surface temperatures over the Tibetan Plateau
According to the comparison of different LST algorithms and in situ LST, the RFR model can provide more accurate LST results while it is cost-computing when compared to the LR model. Considering the limitations of the computing abilities of the GEE platform, we provide a GEE repository with the code required to derive LST by the LR model to obtain the optimized annual average LST in the study area. With the source code, the GEE platform can be used to calculate land surface characteristic parameters (e.g. LST, LSE, and NDVI) with a high spatial resolution efficiently. GEE also allows users to directly retrieve the land surface parameters of their study area by simply modifying the geometry of the region of interest (ROI) variable in the code. Additionally, with recent improvements in computing ability, the RFR model is easily reproducible. The original RFR model files have also been provided and can be used to obtain higher-precision LST data.
Retrieved from a total of 109,978 ETM+ images, the annual average LST result from 2006 to 2017 with a 30 m spatial resolution for the entire TP is shown in Figure 8. The LST on the TP is markedly lower than that in the surrounding area, primarily due to the higher elevation of the TP. The Qaidam Basin has a higher LST than the primary part of the plateau, primarily because it has a lower altitude than a desert cover type.

Discussion
This study presents a method to retrieve the LST over the entire TP from the Landsat 7 ETM+ TIR band with GEE. The generalized SC algorithm was modified and improved by the LR model to obtain LST with a fine spatial resolution of 30 m. The quality of the LST retrievals was also assessed by comparisons with in situ measurements from eight stations over the TP and MOD11A1 LST products. Results showed that LST retrieved from the improved SC algorithm was more stable and accurate.
While training the RFR model, the accuracy of the final model improves as the size of the training set increases. A small training set will reduce the stability of the model, making the model unable to distinguish different situations. However, due to the complex underlying surface status and harsh environment of the TP, the stations with radiation measurement capabilities on the TP are not distributed as densely as those in the plain areas. The data from the 8 stations used in this study (Table 2) are all the data that we could obtain currently. As shown in Figure 9, the surface types of the TP are primarily divided into grassland, forest and bare land, which together occupy 95.79% of the area. The 8 stations used in this study covered these three primary underlying surface types, which can represent the performance of areas far away from selected stations to a certain extent.
For LST validation, errors may come from both the retrieval process of satellite LST and calculation of in situ measured LST. In this study, data from Landsat ETM+ and TIRS were upscaled to 1 km to match the spatial resolution of the MODIS LST product, which may bring uncertainties. A more reliable method is to establish multiple observation stations within a MODIS pixel and use these stations to make a comprehensive evaluation of MODIS product. However, due to the complex underlying surface status and harsh environment of the TP, meteorological stations are both sparse and unevenly distributed, which is not feasible for the time being. Instead, we analysed the homogeneity of the land cover within 1 km radius around the TORP stations. As shown in Figure 10, except for SETORS and NAMORS, nearly all stations are surrounded by only two types of land cover: bare soil and grassland. As shown in Table 5, the land cover types of the stations (Bareland for NADORS; Grassland for BJ, MAWORS, NAMORS, NPAM and D105) account for   more than 85% of the area within 1 km radius around most stations except for the SETORS and QOMS stations. Therefore, the land cover around the stations is homogeneous to a certain extent. However, different microtopographies and soil moisture conditions may result in uncertainties that are not easy to quantify.
The uncertainties in the retrieved LST may also come from the ERA5 TCWV in AFs (Equation (8)). When water vapour is greater than 3 g · cm −2 , the error in the algorithm in the retrieved results become large Jiménez-Muñoz and Sobrino 2010). In this study, the RMSE between the retrieved LST and the in situ measurements is approximately 4 K for regions with TCWV greater than 3 g · cm −2 , and the RMSE of LST increases by 0.5 K when the RMSE of TCWV is 0.35 g · cm −2 . Additionally, the gap filling method used to compensate for the image gaps caused by Landsat 7 SLC-off may also introduce errors. To validate the accuracy of the gap-filling method, LST results from 2008 to 2018 at BJ station were selected as an example. As shown in Figure 11, the RMSE between gap-filled LST and in situ measurements is 3.034 K, which is larger than that of LST without gaps (2.091 K).

Conclusion
In this study, we proposed RFR models accounting for different land cover types that were trained by in situ measurements. Compared with the LST results from the SC algorithm, the LR model and MOD11A1, the results estimated by the trained RFR models are the most accurate, with an RMSE of 1.890 K. Another novel aspect of this study is the use of the GEE platform. We provide a GEE repository with the improved SC algorithm to retrieve the LST for the TP area. Using Google's online computing capability, it was possible to derive high-resolution land surface characteristic parameters at a large scale in a relatively short time. The primary conclusions of this study are summarized as follows: (1) The coefficients of AFs in the SC algorithm that are more suitable for the TP area were obtained by the LR model based on in situ measurements from eight stations. The cross validation between the LR LST results and the test set of in situ measurements showed that the retrieved LST from ETM+ exhibited good agreement with in situ measurements with an RMSE of 2.767 K.
(2) A strategy to estimate LST was proposed by training multiple remote sensing indices using the RFR model. Different RFR models were established according to the different LCTs to distinguish the LST characteristics of different LCTs. The RFR method achieved better accuracy in LST estimation with an RMSE (rRMSE) of 1.890 K (0.164) when compared with the MODIS LST product with an RMSE (rRMSE) of 3.625 K (0.315), USGS Landsat 7 Level-2 LST product with an RMSE of 4.493 K (0.386) and LR model with an RMSE of 2.767 K (0.241).
(3) We applied the improved SC algorithm based on the LR model to the GEE platform. A total of 109,978 ETM+ images acquired from 2006 to 2017 over the TP were processed with the improved SC algorithm. The spatial distributions of the annual average LST, NDVI and LSE with a spatial resolution of 30 m for the entire TP were obtained. Using GEE, large-scale and high-precision land surface characteristic parameters can quickly be calculated. The code provided in this study allows researchers to process data quickly and efficiently.