Dataset of daily near-surface air temperature in China from 1979 to 2018

Near-surface air temperature (T_a) is an important physical parameter that reflects climate change. Many methods are used to obtain the daily maximum (T_max), minimum (T_min), and average (T_avg) temperature, but are affected by multiple factors. To obtain daily T_a data (T_max, T_min, and T_avg) with high spatio-temporal resolution in China, we fully analyzed the advantages and disadvantages of various existing data. Different T_a reconstruction models were constructed for different weather conditions, and the data accuracy was improved by building correction equations for different regions. Finally, a dataset of daily temperature (T_max, T_min, and T_avg) in China from 1979 to 2018 was obtained with a spatial resolution of 0.1^∘. For T_max, validation using in situ data shows that the root mean square error (RMSE) ranges from 0.86 to 1.78^∘, the mean absolute error (MAE) varies from 0.63 to 1.40^∘, and the Pearson coefficient (R²) ranges from 0.96 to 0.99. For T_min, the RMSE ranges from 0.78 to 2.09^∘, the MAE varies from 0.58 to 1.61^∘, and the R² ranges from 0.95 to 0.99. For T_avg, the RMSE ranges from 0.35 to 1.00^∘, the MAE varies from 0.27 to 0.68 ^∘, and the R² ranges from 0.99 to 1.00. Furthermore, various evaluation indicators were used to analyze the temporal and spatial variation trends of T_a, and the T_avg increase was more than 0.03 ^∘C yr⁻¹, which is consistent with the general global warming trend. In summary, this dataset has high spatial resolution and high accuracy, which compensates for the temperature values (T_max, T_min, and T_avg) previously missing at high spatial resolution and provides key parameters for the study of climate change, especially high-temperature drought and low-temperature chilling damage. The dataset is publicly available at https://doi.org/10.5281/zenodo.5502275 (Fang et al., 2021a).

1 Introduction

Near-surface air temperature (T_a) is an important variable that reflects global climate change and significantly affects the cyclical conversion of energy and matter in all spheres of the earth (Gao et al., 2012, 2014). Obtaining accurate grid T_a is helpful for research on urban heat island effects, ecological environment changes, vegetation phenology development, crop yield fluctuation, and energy dynamic balance (Lin et al., 2012; Bolstad et al., 1998). In this study, T_a refers to the daily maximum (T_max), minimum (T_min), and average temperatures (T_avg) of daily near-surface air temperature, which are important input parameters for hydrological, environmental, and crop models (Han et al., 2020; He et al., 2020; Mostovoy et al., 2006; Schaer et al., 2004). These parameters can accurately reflect the frequency and extent of the occurrence and development of extreme climate events (Zhang et al., 2017; Miao et al., 2016). With the intensification of global warming, the temperature gradually rises, the number of extremely cold days and cold nights gradually decreases, and the frequency of extreme weather events also increases (Ding et al., 2006; Liao, 2020; Ryoo et al., 2010). China is a country where extreme weather events frequently occur, causing substantial economic losses (Kharin et al., 2007; Kong, 2020). Therefore, obtaining spatio-temporal changes in T_a is necessary to study extreme weather events and meteorological disasters leading to decreased agricultural yield.

T_a is affected by many factors of the Earth's system, resulting in frequent, complicated diurnal temperature fluctuations (Schwingshackl et al., 2018; Chen et al., 2014). At present, T_a is obtained mainly via three methods: T_a observed via meteorological stations, T_a estimated from land surface temperature (T_s) retrieved from remote sensing, and T_a obtained from the assimilation model. Temperatures with high temporal resolution can be obtained via measurements from meteorological stations. This detection method can avoid the influence of clouds and rain, preserving relatively good data integrity, continuity, and accuracy. However, the number of meteorological stations is limited and they are unevenly distributed, especially for mountainous regions (Mao et al., 2008; Gao et al., 2018; Zhao et al., 2020). Most meteorological stations are in sparsely populated areas far from cities and cannot accurately monitor changes in urban temperature caused by the urban heat island effect (He and Wang, 2020). Moreover, due to the aging of meteorological station equipment, the observation data may be incomplete. Although many interpolation methods, such as Kriging, cubic spline, and inverse distance weight interpolations are available, the difference in density among stations affects the interpolation accuracy (Tang et al., 2020; Berezowski et al., 2016; Tencer et al., 2011).

Satellite sensors provide global coverage and high-spatial-resolution data used to estimate T_a. The most commonly used estimation methods are the statistical regression method (Wen, 2020; Zhu et al., 2013; Zhang et al., 2015), the temperature vegetation index method (Xing et al., 2020), the energy balance method (Benali et al., 2012), the atmospheric temperature profile extrapolation method (Wen, 2020), and the machine learning method (Mao et al., 2008; Wen, 2020). Sensors are susceptible to weather phenomena, such as clouds and rain, leading to missing data or reduced quality. In addition, these methods are mostly suitable for clear-sky conditions, which need to be further expanded to establish an estimation model of T_s to T_a under different weather conditions.

Reanalysis data generated by the global assimilation model have provided many datasets of geophysical parameters, including near-surface temperature, which overcome most of the aforementioned problems caused by abnormal weather. The NCEP/NCAR reanalysis dataset was developed by the National Center for Environmental Prediction and the National Center for Atmospheric Research (January 1948–September 2021), with a temporal resolution of 6 h and a spatial resolution of 2.5^∘ (Kalnay et al., 1996). The ERA5 dataset was released by the European Center for Medium-Range Weather Forecast (ECMWF; January 1950–September 2021), with a temporal resolution of 1 h, and a spatial resolution of 0.3^∘ (Hersbach et al., 2020; Dee et al., 2011; Taszarek et al., 2021; Lei et al., 2020). The land surface modeling forcing dataset was developed by Princeton University (January 1948–December 2006), with a temporal resolution of 3 h and a spatial resolution of 1.0^∘ (Deng et al., 2010). To improve the accuracy of regional data, some researchers have developed meteorological forcing datasets for China. The representative dataset is the China Meteorological Forcing Dataset (CMFD) released by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (January 1979–December 2018), with a temporal resolution of 3 h and a spatial resolution of 0.1^∘ (He, 2010; Yang et al., 2010; Yang and He, 2019). However, the dataset does not provide daily maximum and minimum temperatures. The grid dataset of daily surface temperature in China (V2.0) was released by the China Meteorological Administration (CMA; January 1961–September 2021), with a spatial resolution of 0.5^∘. This dataset comprises the daily maximum, minimum, and average temperatures; its spatial resolution is low; and the accuracy of local areas needs improvement. Although reanalysis datasets can obtain global near-surface air temperature data, the number of T_max, T_min, and T_avg datasets with high spatial resolution and high precision is insufficient.

In this study, we aimed to obtain a long-term T_a (T_max, T_min, and T_avg) dataset with high spatial resolution in China. We first analyzed the advantages and disadvantages of various data (e.g., reanalysis, remote sensing, in situ data). Next, we constructed daily T_a models for clear- and non-clear-sky conditions. This method compensates for the deficiency that studies have estimated T_a mostly under clear-sky conditions rather than under all-sky conditions. We further improve data accuracy by building correction equations for different regions. Finally, a dataset of daily T_a (T_max, T_min, and T_avg) in China from 1979 to 2018 was obtained with a spatial resolution of 0.1^∘, and we cross-validated this dataset with existing datasets.

2 Study area

China's vast territory has significant undulations on the Earth's surface, and a wide range of climate changes. To explore the temporal and spatial characteristics of T_a, we divided China into six subregions (Fig. 1) according to climatic conditions, such as temperature and rainfall, and topographical conditions, such as elevation. (I) The northeastern region mainly includes Northeast China, located to the east of the Greater Khingan Range. This region is located in the temperate monsoon climate zone, the annual precipitation is 400–1000 mm, and cumulative temperature is between 2500 and 4000 ^∘C (Mao and Wan, 2000). (II) The North China region is located north of the Qinling-Huaihe River and south of the Inner Mongolia Plateau. This region is mostly located in the temperate monsoon climate zone, and the annual accumulated temperature is between 3000 and 4500 ^∘C (Xu et al., 2017), with hot, rainy summers and cold, dry winters. (III) The central southern region is located south of the Qinling-Huaihe River and north of the tropical monsoon climate type. This region is located in the subtropical monsoon climate zone, the annual accumulated temperature is between 4500 and 8000 ^∘C, and the precipitation is mostly between 800 and 1600 mm. (IV) The southern region is south of the Tropic of Cancer. This region is located in the tropical monsoon climate zone, the annual accumulated temperature is greater than 8000 ^∘C, the annual minimum temperature is not less than 0 ^∘C, and there is no frost year round. Annual precipitation mostly ranges from 1500 to 2000 mm. (V) The northwest region is mainly distributed in the inland areas above 40^∘ N latitude in China, located northwest of the Greater Khingan Range-Yin Shan-Ho–lan Mountains-Qilian Mountains line. This region is far from the coast, water vapor transport is limited, annual precipitation is between 300 and 500 mm, and the annual accumulated temperature is between 2000 and 3500 ^∘C. The daily and annual temperature differences are large, including those in the temperate desert, temperate grassy, and subfrigid coniferous climates. (VI) The Qinghai-Tibet Plateau region includes the Qinghai-Tibet Plateau, the Andes Mountains, Mount Everest, and other areas. This region is located in the plateau and mountainous climate zone, the annual accumulated temperature is lower than 2000^∘, the daily temperature range is large, and the annual temperature range is small. This region has strong solar radiation, sufficient sunshine, and little precipitation.

Figure 1 Scope map of the total study area and the six subregions. Black dots indicate distribution locations of meteorological stations; blue frame lines indicate the substudy area range, represented by I, II, III, IV, V, and VI.

3 Data

3.1 Reanalysis data

The reanalysis dataset contains drivers of surface elements in a large area, which can provide highly complementary information and avoid data gaps and low-quality pixels caused by abnormal weather conditions. This study primarily used the CMFD and ERA5 data as reanalysis data sources.

The CMFD is a set of meteorological forcing datasets developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (He et al., 2020; Yang et al., 2010; Yang and He, 2019). It is mainly based on the Global Land Data Assimilation System (GLDAS) as a background dataset, using empirical knowledge algorithms and combining GLDAS with measured data to obtain temperature data with a spatial resolution of 0.1^∘. The CMFD contains seven variables: 2 m air temperature, surface pressure, specific humidity, 10 m wind speed, downward shortwave radiation, downward longwave radiation, and precipitation rate. The CMFD covers January 1979 to December 2018 and provides four types of temporal resolution (3 h, daily, monthly, and yearly). The CMFD is comprehensive and has the longest time series and the highest spatial resolution in China. Studies have used the temperature data as input parameters to construct a surface air temperature model, which shows that the correlation coefficient between the CMFD temperature and the measured data is greater than 0.99 and has high consistency, and that grid data can reflect the temporal and spatial changes in regional air temperature (Zhang et al., 2019; Wang et al., 2017). The CMFD as an input element to build a surface temperature model can also significantly reduce model deviation and improve model accuracy (Chen et al., 2011). Therefore, we used the 3 h temperature of the CMFD to build the T_a Model and verified the new product with the daily temperature from the CMFD. The CMFD is available from the China National Qinghai–Tibet Plateau Science Data Center (http://data.tpdc.ac.cn/zh-hans/data/8028b944-daaa-4511-8769-965612652c49/, last access: 1 November 2020).

ERA5 is the fifth-generation product of the atmospheric reanalysis global climate data launched by the ECMWF, replacing the ERA-Interim reanalysis data which were discontinued on 31 August 2019. ERA5 data are generated based on the Cy41r2 model of the integrated forecasting system, which has benefited from the development of data assimilation, model simulation, and model physics, and is generated by assimilating many ground-monitoring, aircraft weather observation, and radio-detection data. ERA5 data are significantly better than ERA-Interim data; for example, the former has a higher spatio-temporal resolution, more vertical mode levels, and more parameter products than the latter. ERA5 provides timely, updated quality checks on the data, which is convenient for providing stable, real-time, and long-term climate information. ERA5 provides many meteorological elements, including 2 m air temperature, 2 m relative humidity, sea level pressure, sea surface temperature, and precipitation. Since the release of the ERA5 reanalysis data, many researchers have tested their applicability and accuracy. The results show that the accuracy of the ERA5 is better than that of the ERA-Interim data, and the higher spatio-temporal resolutions are conducive to the precise description of regional atmospheres. The details of these improvements are convenient for studying changes in small-scale atmospheric environments (Meng et al., 2018; Mo et al., 2021; Hillebrand et al., 2021). These data can be obtained from https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset&text=ERA5 (last access: 1 December 2020).

3.2 In situ data

The in situ data from 1979 to 2018 used in this study were employed to build a T_a model and evaluate existing datasets and new products. The measured data of meteorological stations were from the China National Meteorological Information Center (http://www.nmic.cn/site/index.html, last access: 1 November 2020), including hourly air temperature, hourly land surface temperature, maximum daily temperature (T_max), minimum daily temperature (T_min), daily average temperature (T_avg), and weather condition records. Due to the inconsistency of recorded data of meteorological conditions at many stations, some data are missing. Furthermore, there are no meteorological stations in most areas; thus, the data are used as auxiliary data.

The ground observations obtained from the China Meteorological Administration underwent uniform data processing and homogeneity testing. To further ensure the quality of the data, we checked the in situ data. First, we set a fixed threshold to eliminate the overflow value. Second, we tested the time series of station data and eliminated abnormal and missing data due to instrument damage or bad weather (Zhao et al., 2020). Finally, we checked the spatio-temporal consistency of the in situ data, deleted the meteorological stations with location migration during the study period, and maintained the temperature data of meteorological stations with a long monitoring time and stable temperature values.

3.3 Supplementary data

China's daily near-surface temperature grid dataset was released by the CMA with a spatial resolution of 0.5^∘. This grid dataset contains the daily maximum, minimum, and average temperatures in China (http://www.nmic.cn/site/index.html, last access: 11 April 2021). The CMA dataset was obtained by combining the daily temperature data monitored by meteorological stations and the digital elevation model (DEM) data generated by resampling with three-dimensional geospatial information via a thin-plate spline interpolation algorithm. The spatial resolution of the CMA data is 0.5^∘, and we used these data for cross-validation. The Moderate Resolution Imaging Spectroradiometer (MODIS) is an important sensor in the Earth Observation System program and is mounted on the Terra and Aqua satellites. Terra is a morning-orbiting satellite that passes through the Equator at approximately 10:30 local time from north to south. Aqua is an afternoon-orbiting satellite that passes through the Equator at approximately 01:30 local time from south to north. The Terra satellite has been in service since 1999, the Aqua satellite since 2002. Since 2002, surface temperature data can be obtained four times per day from MODIS data through inversion calculation. In this study, we used the MOD11A1 and MYD11A1 products, which provide daily surface temperature data on a global scale with a spatial resolution of 1 km. MODIS LST (land surface temperature) has a quality control (QC) field that indicates data quality and is encoded in binary form. MODIS data can be downloaded from the LAADS DAAC (Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center) website (https://ladsweb.modaps.eosdis.nasa.gov/search/order, last access: 1 December 2020).

In addition to the aforementioned data, DEM data were used. The Shuttle Radar Topography Mission (SRTM) DEM used in this study was a radar topographic mapping project jointly implemented by NASA and the National Imagery and Mapping Agency, which was implemented by the Space Shuttle Endeavour. Temperature data were regulated via the topographical correction of the SRTM DEM, with 90 m resolution to eliminate the influence of topographical fluctuations on air temperature. SRTM DEM data can be obtained from the Geospatial Data Cloud (http://www.gscloud.cn/search, last access: 10 February 2021).

4 Methodology

The T_max, T_min, and T_avg data were provided by meteorological stations. Other non-station locations or grid values were estimated by interpolation or indirect methods such as remote sensing. Because of the limited number of meteorological stations and their uneven distribution, it is difficult to guarantee the accuracy of T_max, T_min, and T_avg obtained through interpolation in some areas. Under rainfall and cloud-cover weather conditions, estimating the air temperature from remotely sensed surface temperature data is impossible. Even in clear-sky conditions, the formula for estimating near-surface air temperature is not universally applicable, which hinders the development of a high-precision T_a dataset to a certain extent. Therefore, to obtain a T_a dataset with high spatio-temporal resolution and long time series, it is necessary to build a reliable and robust T_a model to estimate T_max and T_min, and further improve the accuracy of T_avg. Consequently, the product could be widely used for climate change and research on extreme weather events.

Daily temperature changes are affected by many factors and are extremely sensitive to fluctuations under different weather conditions. This study used multiple methods to calculate T_a. First, the daily weather conditions were divided into clear-sky and non-clear-sky conditions. Second, based on the physical process of daily temperature changes and combined with existing reanalysis data, in situ data, and remote-sensing data, we estimated T_max and T_min under different weather conditions. To further improve the accuracy of the dataset, we constructed a modified model for each region. Details are provided in the following sections. The overall process of this study is illustrated in Fig. 2. The construction of the dataset was mainly divided into three steps: (1) the process of daily weather condition determination, (2) the process of establishing T_a models under different weather conditions, and (3) data correction.

Figure 2 Summary flowchart of T_a dataset establishment.

4.1 Strategies for division of weather conditions and T_a estimation

4.1.1 Scheme for dividing weather conditions

Different weather conditions have different rules of temperature change. To improve the estimation accuracy of the maximum and minimum temperature, we conducted specific calculations by distinguishing daily weather conditions. The quality of observation data is affected by weather, and some remote-sensing products, such as MODIS LST products, have quality control fields. Therefore, the quality control field of MODIS can be used to distinguish between clear-sky and non-clear-sky conditions. However, we have only been able to obtain MODIS observation data four times per day since 2002, which cannot cover the timespan involved in this study. Therefore, we divided the time series of this study into two periods: 1979–2001 and 2002–2018, and different methods were used for the two time series to distinguish the daily weather conditions. For the study period from 2002 to 2018, we distinguished each pixel based mainly on the MODIS quality control field. When the MODIS quality control of all four T_s corresponding to a pixel was in the clear-sky condition, the pixel was judged to be in the clear-sky condition; otherwise, it was judged to be in the non-clear-sky condition.

For the study period from 1979 to 2002, we used the in situ, CMFD, and ERA5 data to determine the daily weather condition. First, we filtered each pixel and divided it into two types: meteorological stations corresponding to pixels with and without weather condition records. For pixels with weather condition records, we used many statistical discrimination methods to analyze the impact of non-clear-sky weather phenomena on temperature fluctuations, which can facilitate the subsequent determination of pixels without weather condition records. Statistical analysis shows a significant difference in daily temperature fluctuations between clear-sky and non-clear-sky conditions, and non-clear-sky weather conditions may cause abnormal temperature fluctuations. Therefore, we converted the judgment of the weather state into the abnormal judgment of the time and frequency of the occurrence of T_max and T_min (occurrence time of T_max and T_min is hereinafter cited as H_max and H_min, respectively). Specifically, when H_max and H_min occur abnormally or the temperature change is wavy, a non-clear-sky condition is used (Zhao and Duan, 2014; Ren et al., 2011). In other cases, they are regarded as clear-sky conditions, and the position of each pixel is marked. Therefore, we had to further fill the daily time series of each pixel to determine the weather condition. In this study, we used two strategies to perfect the temperature series for distinguishing weather conditions. The specific implementation steps for determining weather conditions are shown in Fig. 3.

Figure 3 Summary flowchart for classification of weather conditions.

In the first strategy, when the pixel location had a corresponding meteorological station or when the Euclidean distance between adjacent stations was less than 0.3^∘, we filled in the gaps to improve the integrity and continuity of the time series. The time series-filling process was as follows: (1) when the temperature data at the observation sites were missing and not consecutively missing, in the case of the same spatial range, we used the average temperature of two adjacent timepoints before and after the missing value at the same site to fill in the missing value; and (2) when the observation data of a station were continuously missing, in the same time range, we filled the missing value with the observation data of the stations within 0.3^∘. This method was mainly based on the principle that the closer the distance between stations, the stronger the spatial consistency and correlation of temperature changes. (3) When the station data were continuously missing and the adjacent station data could not be filled, other relevant data were used for repair within the same time and space. In this study, we estimated the weather state from the T_s monitored by the same station. This method theoretically originates from the approximate consistency between the daily variation ranges of T_s and T_a and is suitable for situations where there are many missing values and incomplete time series at meteorological stations and adjacent meteorological stations. Many studies have analyzed the correlation between the daily trend of T_a and T_s and found strong consistency. The T_s retrieved by remote-sensing satellites is also widely used to estimate T_a, which proves the reliability of determining the pixel weather state through the T_s time series (He et al., 2020; Yoo et al., 2018; Johnson and Fitzpatrick, 1977; Caesar et al., 2006; Mostovoy et al., 2006). (4) When there is no meteorological station at the pixel location and the distance from the meteorological station is less than 0.3^∘, we use the inverse distance weighting method to perform spatial interpolation on adjacent pixels. Notably, before interpolation, we need to consider the impact of elevation differences. To improve the interpolation accuracy, we first correct the data of the observation station to a uniform sea level, and then perform further calculations according to the elevation of the interpolation point to obtain the corresponding temperature.

The second strategy was to target areas where the distribution of stations was sparse and the Euclidean distance between two adjacent stations was greater than 0.3^∘. To compensate for the insufficient coverage and uneven distribution of stations in these areas, we used hourly data from ERA5 to determine the approximate time of occurrence of T_max and T_min. Because of a certain difference between the spatial resolution of ERA5 and this dataset, it was difficult to fulfill our demand for higher spatial resolution. Consequently, we developed an effective downscaling process based on the spatial correlation between the ERA5 data and CMFD temperature data. ERA5 data (with a spatial resolution of 0.3^∘) were spatially downscaled with the aid of the CMFD data (with a spatial resolution of 0.1^∘). The downscaling process is illustrated in Fig. 4. First, quality control of the ERA5 data and CMFD was performed to eliminate temperature outliers. Second, the ERA5 data and CMFD were matched according to time series and central latitude and longitude to construct pixel pairs. Subsequently, we weighted the high-resolution data to the low-resolution ERA5 data pixel by pixel. Finally, the weight was used to downscale the ERA5 data to the same spatial resolution of the CMFD. The ERA5 downscaling was computed using Eqs. (1) and (2),

$$\begin{array}{}\text{(1)}& {\displaystyle }{T}_{\mathrm{E}}\left(x_o,y_o\right)={\displaystyle \frac{T_{\mathrm{C}}\left(x_o,y_o\right)}{\sum _{i=\mathrm{1}}^m\sum _{j=\mathrm{1}}^n{T}_{\mathrm{C}}\left(x_i,y_j\right)}}\cdot T_{\mathrm{E}}\left(x_m,y_n\right),\text{(2)}& {\displaystyle }{T}_{\mathrm{E}}\left(x_o,y_o\right)={\displaystyle \frac{T_{\mathrm{M}}\left(x_o,y_o\right)}{\sum _{i=\mathrm{1}}^m\sum _{j=\mathrm{1}}^n{T}_{\mathrm{M}}\left(x_i,y_j\right)}}\cdot T_{\mathrm{E}}\left(x_m,y_n\right),\end{array}$$

where T_E,T_C, and T_M represent the ERA5 data, CMFD, and MODIS data, respectively. T_E(x_o,y_o) is the temperature data after downscaling; T_E(x_m,y_n) is the temperature data before downscaling; i, j are pixel coordinates; and m, n are the pixel coordinates before downscaling.

Figure 4 Flowchart for spatial downscaling, where nv represents the number of valid values.

4.1.2 T_max and T_min estimation under clear-sky conditions

In addition to the severe temperature fluctuations caused by abnormal weather phenomena, the daily temperature changes under clear-sky conditions have a certain regularity, periodicity, and asymmetry (Leuning et al., 1995; Johnson and Fitzpatrick, 1977). According to the similarity between the surface temperature and the diurnal variation trend of air temperature, a method of estimating T_a is established by the daily air temperature variation model. Verified by in situ data, this method is feasible (Du et al., 2020; Zhu et al., 2013; Perkins et al., 2007; Cesaraccio et al., 2001; Serrano-Notivoli et al., 2019). However, using the surface temperature retrieved by remote-sensing methods to estimate the changing trend of air temperature is complicated, additional parameters need to be input, and the relationship between T_s and T_a is not fixed. Therefore, it is difficult to unify the types and quantities of parameters and ensure accuracy. Thus, we established a piecewise local sine function of temperature under clear-sky conditions for each pixel, which can simulate the change in T_a and calculate T_max and T_min (Mao et al., 2016; Jiang et al., 2010). First, according to the approximate periodicity of daily temperature changes and the asymmetry of H_max and H_min, we derived the T_a piecewise sine function of the adjacent regions of H_max and H_min (Eqs. 3 and 4). Second, using a method similar to that in Sect. 4.1.1, we obtained H_max and H_min for each pixel. These H_max and H_min values are entered as parameters into the piecewise sine function. The CMFD (3 h data) are used as T_a data, each pixel H_max and H_min are used as time, and the values of A_t and B_t are obtained by the least squares method. Finally, H_max and H_min values were substituted into the derivation formula to obtain T_max and T_min as preliminary results for subsequent correction and analysis. We constructed a temperature model, pixel by pixel, to fulfill the temporal and spatial heterogeneity of each region.

$$\begin{array}{}\text{(3)}& {\displaystyle }{T}_{max}& {\displaystyle }=A_t\cdot \sin \left[{\displaystyle \frac{(H_o-H_{\max })\mathit{\pi }}{H_{\max }-H_{\min }}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right]+B_t,\text{(4)}& {\displaystyle }{T}_{min}& {\displaystyle }=A_t\cdot \sin \left[{\displaystyle \frac{(H_o-H_{\max })\mathit{\pi }}{\mathrm{24}-H_{\max }+H_{\min }}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\right]+B_t,\end{array}$$

where H_max is the occurrence time of the daily maximum temperature, H_min is the occurrence time of the daily minimum temperature, H_o is the input time, and A_t and B_t are unknown parameters.

4.1.3 T_max and T_min estimation under non-clear-sky conditions

The daily temperature fluctuations in non-clear-sky conditions are relatively large, and there may be large-scale cooling or sudden temperature changes within a short period. Based on the spatial location information of each pixel, the most reliable and representative data source are the in situ data. Therefore, if there are in situ data for the pixel location, the temperature data at the same time will be directly obtained from the station to replace the pixel values T_max and T_min. For the pixels corresponding to non-meteorological stations, similar to the method of spatial downscaling for the pixel positions of non-meteorological stations in the weather condition judgment, we used ERA5 data to perform spatial downscaling with the assistance of the CMFD. By adding high-spatial-resolution MODIS data, the downscaling method was further expanded to improve the accuracy of each pixel. We mainly wanted to fully exploit the advantages of various data, especially with the help of high-resolution MODIS data. According to the QC field of MODIS data, we used MODIS data with high spatio-temporal resolution to improve local accuracy while ensuring high-quality MODIS data. The corresponding time of the effective pixel was matched with the ERA5 data according to the nearby time, to obtain the data weight for spatial downscaling. The downscaling process and the validity determination of MODIS data are shown in Fig. 4, and the downscaling formulas are shown in Eqs. (1) and (2).

4.1.4 T_avg estimation

Usually, the aim of calculating average temperature is to use the temperature value observed every day to obtain an arithmetic average. If each pixel has hourly temperature data, the calculated daily average temperature is the most representative. Because the observational conditions are limited, hourly temperature data is difficult to obtain; thus, often, the temperature values of four observation times (e.g., 02:00, 08:00, 14:00, and 20:00) are used to obtain the daily average temperature, or the daily maximum and minimum temperatures are directly averaged to obtain the daily average temperature. To improve the accuracy of the average temperature as much as possible, we used the 3 h temperature data provided by the CMFD and the maximum and minimum values we have calculated to conduct an arithmetic average to obtain the daily average temperature. Finally, to improve the accuracy, we performed multiple linear regression correction on the T_avg output value according to the in situ data (the linear correction method was the same as that described in Sect. 4.2) and obtained the daily T_avg dataset.

4.2 T_a data calibration scheme

Surface temperature is sensitive to changes in altitude and easily affected by the surrounding environment. For non-meteorological station pixels, we use interpolation to fill in the pixel values based on the principle of regional consistency. To improve the accuracy of the pixel temperature at non-meteorological stations, we fully considered the influence of altitude on temperature. First, the in situ T_a was unified to sea level according to the vertical rate of temperature drop. Next, the non-station pixels were interpolated according to the station data, and finally, the interpolated pixel values were restored to the corresponding elevation. This method can reduce the influence of altitude on temperature to a certain extent and improve the accuracy of the dataset. In this study, we used a uniform vertical temperature drop rate (γ), i.e., for every 100 m increase in altitude, the atmospheric temperature decreases vertically by 0.65 ^∘C, and vice versa. The height correction formula is provided by Eq. (5) (He and Wang, 2020; Schicker et al., 2015; Wang, 2013):

$$\begin{array}{}\text{(5)}& T_{\mathrm{SL}}=T_{\mathrm{a}}-\mathit{\gamma }\cdot \left(H_{\mathrm{SL}}-H_{\mathrm{a}}\right),\end{array}$$

where T_SL is the sea-level temperature, T_a is the temperature of the meteorological station, and H_SL is the sea-level height, where the value of γ is approximately 0.0065 ^∘C m⁻¹.

We used the jackknife method as follows: 699 in situ stations across China were divided into 140 verification points and 559 calibration points according to the ratio of 20 to 80 to establish a multiple linear regression equation (Benali et al., 2012; Xu et al., 2017). The preliminary accuracy results (Sect. 5.1) show that although the overall accuracy was high, there remains the problem of abnormal temperature values of the model output data caused by the violent fluctuations in daily temperature changes. Further correction is required to reduce the deviation and improve the accuracy of the dataset. The data correction process is illustrated in Fig. 5. For the abnormal temperature value, we replaced the T_a at the pixel location with the observation T_a from the meteorological station and performed the adjacent pixel temperature correction for the pixel without the meteorological station at the pixel location. The multiple linear regression method was used to process the original temperature, and the stepwise regression relationship between the measured value of the station and the fitted value of the corresponding pixel was established. Next, we calculated the predicted value of the regression temperature according to the regression equation and obtained the temperature residual value by calculating the observed value and the predicted value to obtain the final corrected temperature (Cristobal et al., 2006). The modified expression is shown in Eq. (6):

$$\begin{array}{}\text{(6)}& {\displaystyle }V\left(x,y\right)=\widehat{m}\left(x,y\right)+\widehat{\mathit{\epsilon }}\left(x,y\right),\end{array}$$

where x and y are the numbers of rows and columns of pixels, respectively; V(x, y) is the correction value of the regression equation; $\widehat{m}\left(x,y\right)$ is the regression prediction value of air temperature; and $\widehat{\mathit{\epsilon }}\left(x,y\right)$ is the residual value.

Figure 5 Flowchart for calibration of T_a model data.

4.3 Evaluation metrics

We mainly selected areas with a single surface type and flat terrain under clear skies as the comparative study area to verify the original dataset and reconstructed dataset. A scatter diagram can represent the overall distribution and aggregation of the data and intuitively convey accurate information from the data; thus, we used a scatter chart to display the accuracy range of this product. In addition, before establishing the model, we retained a part of the reanalyzed data excluded from the calculation and used it for cross-validation. We used three indicators as metrics to measure the accuracy of variables, i.e., R², MAE, and RMSE.

We compared T_max and T_min with the ERA5 data and CMA data. Notably, the ERA5 reanalysis dataset is an hourly temperature grid dataset; thus, we obtained the highest and lowest temperature values of ERA5 by constructing a local sine function similar to that in the prior section and further calculated the average daily temperature. The accuracy of T_avg products in this study was verified with the ERA5 data, CMA data, and CMFD daily temperature data. Because the spatial resolution of CMA is 0.5^∘, to facilitate comparison, we resampled the spatial resolution of all datasets to 0.5^∘.

4.4 Analysis of the T_a series trend

We not only compared the output T_a data with the in situ data, but also assessed the climate change trends of T_max, T_min, and T_avg in various regions of China, and further tested the effectiveness and regional applicability of the dataset through various climate variables. The World Meteorological Organization defined a series of extreme climate indexes, including 27 core indexes. We used four of them (TXx, TNn, TX90p, and TN10p) to analyze the trend of extreme temperature changes in T_max and T_min (Karl et al., 1999; Peterson et al., 2001). Specifically, the TXx (TNn) anomaly refers to the difference between the sum of monthly T_max (T_min) and the multi-year average of monthly T_max (T_min) in each year. The multi-year period of this study is 40 years. In addition, linear regression was performed on the TXx (TNn) anomaly to analyze the interannual variation trend. The TX90p (TN10p) means that the daily T_max (T_min) of each month during the study period is arranged in ascending order, and the 90 % (10 %) corresponding value in the time series is used as the threshold for judging warm days (cold nights; Zhang et al., 2005).

To study the spatio-temporal variation trend of T_avg, we used linear regression analysis (K), correlation coefficient analysis (R), and the T test (Du et al., 2020; Yan et al., 2020; Cao et al., 2021). The interannual change rate and correlation of T_avg were calculated by K and R, and the formulae are provided by Eqs. (7) and (8), respectively. We performed a two-tailed significance test on the T test to measure the significance of the temperature and time series changes (Eq. 9):

$$\begin{array}{}\text{(7)}& {\displaystyle }K& {\displaystyle }={\displaystyle \frac{n{\sum }_{i=\mathrm{1}}^n\left(i{T}_i\right)-{\sum }_{i=\mathrm{1}}^{n}i{\sum }_{i=\mathrm{1}}^n{T}_i}{n{\sum }_{i=\mathrm{1}}^n{i}^{\mathrm{2}}-{\left({\sum }_{i=\mathrm{1}}^{n}i\right)}^{\mathrm{2}}}},\text{(8)}& {\displaystyle }R& {\displaystyle }={\displaystyle \frac{n{\sum }_{i=\mathrm{1}}^n\left(i{T}_i\right)-{\sum }_{i=\mathrm{1}}^{n}i{\sum }_{i=\mathrm{1}}^n{T}_i}{\sqrt{n{\sum }_{i=\mathrm{1}}^n{i}^{\mathrm{2}}-{\left({\sum }_{i=\mathrm{1}}^{n}i\right)}^{\mathrm{2}}}\cdot \sqrt{n{\sum }_{i=\mathrm{1}}^n{T}_i^{\mathrm{2}}-{\left({\sum }_{i=\mathrm{1}}^n{T}_i\right)}^{\mathrm{2}}}}},\text{(9)}& {\displaystyle }& {\displaystyle }T\mathit{\_}\mathrm{test}\left(R\right)={\displaystyle \frac{R\sqrt{n-\mathrm{2}}}{\sqrt{\mathrm{1}-R^{\mathrm{2}}}}},\end{array}$$

where n represents the total number of years of the time series length, i represents the year, and T_i represents T_avg in the ith year. K>0 indicates that the temperature increases within the time series, and K<0 indicates that the temperature decreases within the time series.

5 Results

5.1 Evaluation of the original product

According to the six subregions in Fig. 1, comparative analyses of this product (T_max, T_min, and T_avg) based on in situ data were conducted. Figure 6 shows the accuracy scatter plot between the original data of T_max and the in situ data. The R² fluctuated from 0.91 to 0.99, the MAE ranged from 1.69 to 2.71 ^∘C, and the RMSE ranged from 2.15 to 3.20 ^∘C. Figure 7 shows the accuracy scatter plot of T_min. The R² fluctuated from 0.93 to 0.97, the MAE ranged from 1.34 to 2.17 ^∘C, and the RMSE fluctuated from 1.68 to 2.79 ^∘C. Figure 8 shows the accuracy scatter plot of T_avg. The R² fluctuated between 0.97 and 0.99, the MAE ranged from 0.58 to 0.96 ^∘C, and the RMSE fluctuated from 0.86 to 1.60 ^∘C. As shown in Figs. 6, 7, and 8, the R² of T_max, T_min, and T_avg, and the temperature measured at the meteorological station, were all greater than 0.90. In general, our method performed well in estimating the daily temperature values. However, due to the impact of complex changes in weather, the distribution of temperature values on certain days is discrete, especially in study areas V and VI. Further corrections are necessary to reduce errors and improve the accuracy of the dataset.

Figure 6 Scatter diagrams of the T_max output from the T_a model against ground station data; statistical accuracy measures (R², MAE, and RMSE) are also indicated.

Figure 7 Scatter diagrams of the T_min output from the T_a model against ground station data; statistical accuracy measures (R², MAE, and RMSE) are also indicated.

Figure 8 Scatter diagrams of the T_avg output from the T_a model against ground station data; statistical accuracy measures (R², MAE, and RMSE) are also indicated.

5.2 Evaluation of the new product

The temperature was further corrected using the linear correction method. The data verification results of T_a after correction are shown in Figs. 9, 10, and 11. The results show that the corrected data had a higher consistency with the in situ data. The fitted and observed temperatures were linearly distributed and gradually approached the regression line, and the outliers were significantly reduced. Figure 9 shows the corrected scatter plot of T_max for each study area. The R² fluctuated from 0.96 to 0.99, the MAE ranged from 0.63 to 1.40 ^∘C, and the RMSE fluctuated from 0.86 to 1.78 ^∘C. Figure 10 shows the corrected scatter plot of T_min for each study area. The R² fluctuated between 0.95 and 0.99, the MAE ranged from 0.58 to 1.61 ^∘C, and the RMSE fluctuated from 0.78 to 2.09 ^∘C. Figure 11 depicts the corrected scatter plot of T_avg in each study area, where R² fluctuated between 0.99 and 1.00, the MAE ranged from 0.27 to 0.68 ^∘C, and the RMSE fluctuated from 0.35 to 1.00 ^∘C. The results show that the distribution of numerical points in each area after the correction was denser, mostly concentrated near the 1:1 line, and the degree of clustering with the measured data was higher than before calibration. Our detailed analysis of the daily temperature in the six study areas demonstrated that the accuracy measurement values differed significantly between the east and west. For example, the accuracy error of study area IV is small, and the accuracy error of study areas VI and V is large, which may be affected by the regional topography and the distribution of meteorological stations. Study area IV is in the tropical monsoon climate zone, affected by latitude and topography, and the temperature is relatively high throughout the year. Moreover, the area is in Eastern China and has densely distributed meteorological stations and relatively flat terrain. Linear correction can significantly improve the agreement between the estimated value and the observed value. Study areas VI and V have the highest RMSE. They are in the Qinghai–Tibet Plateau in Southwest China and Xinjiang in the northwest. Such areas have similar characteristics, such as high altitude, large spatial heterogeneity, and few meteorological stations. This result shows that the temperature has strong spatial heterogeneity. In general, the corrected dataset has higher accuracy than the original dataset, satisfies the spatial heterogeneity of different regions, and better estimates the temperature under different weather conditions.

Figure 9 Scatter diagrams of the original T_max and reconstructed results versus their corresponding ground station data in six natural subregions (I, II, III, IV, V, and VI). Gray points indicate low-quality pixel values in the original T_max data, orange points represent the values in the after-calibration T_max dataset; the statistical accuracy measures (R², MAE, and RMSE) are also indicated.

Figure 10 Scatter diagrams of the original T_min and reconstructed results versus their corresponding ground station data in six natural subregions (I, II, III, IV, V, and VI). Gray points indicate low-quality pixel values in the original T_min data, blue points represent the values in the after-calibration T_min dataset; the statistical accuracy measures (R², MAE, and RMSE) are also indicated.

Figure 11 Scatter diagrams of the original T_avg and reconstructed results versus their corresponding ground station data in six natural subregions (I, II, III, IV, V, and VI). Gray points indicate low-quality pixel values in the original T_avg data, green points represent the values in the after-calibration T_avg dataset; the statistical accuracy measures (R², MAE, and RMSE) are also indicated.

To further verify the robustness and accuracy of this product, Table 1 shows the cross-validation results of this product and other datasets, the mean average precision (MAP) of each region, and that this product has a high regional consistency with other datasets. Study area IV in the tropical monsoon climate zone has the highest accuracy, and study area VI located in the Qinghai–Tibet Plateau region of China has the lowest data accuracy. This result may be because the reanalysis dataset is also affected by the number and distribution of meteorological stations and the spatial heterogeneity. The accuracy and robustness of the product were confirmed from another perspective. The accuracy comparison of each area shows that this product has higher accuracy and spatial representation than other datasets. R² is closer to 1, and MAE and RMSE remain low. Through the accuracy evaluation and data comparison between this product and the existing dataset, we found that our product has a better temperature estimation of each area, and the overall accuracy and accuracy of the dataset are higher.

Table 1 Cross-validation results of this product and other datasets. Values in bold indicate study areas with the highest precision, and values in italics indicate the lowest precision.

Figure 12 Multi-axis diagram of TXx anomaly, TX90p, and T_max linear trend graphs. The broken black line represents the TXx anomaly, the red line represents the linear regression of the TXx anomaly, and the orange histogram represents the TX90p change trend.

Figure 13 Multi-axis diagram of TNn anomaly, TN10p, and T_min linear trend graphs. The broken black line represents the TNn anomaly, the red line represents the linear regression of the TNn anomaly, and the blue histogram represents the TN10p change trend.

Figure 14 Multi-year climate change trends in T_avg. Panel (a) K, as calculated by Eq. (7); (b) R between temperature change and time series development, calculated by Eq. (8); (c) T test (R), calculated by Eq. (9). Panels (a.i), (b.i), and (c.i) represent the distribution of pixel values in the corresponding (a), (b), and (c) spatial images.

5.3 Application of the product for trend analysis

We analyzed temperature changes in various regions of China through extreme climate indexes and change trend values to further test the validity and regional applicability of the dataset. As shown in Figs. 12 and 13, the TXx anomalies and TNn anomalies are consistent in the regional change trend. Although the annual anomalies fluctuated during the study period, they gradually changed from negative to positive. This phenomenon confirmed that the temperature fluctuated and increased, and that the T_max and T_min gradually increased, which is consistent with the global warming trend. The average temperature rise of TXx anomalies in each study area was 0.42 ^∘C a⁻¹, and the average temperature rise of TXx anomalies was 0.47 ^∘C a⁻¹. The histograms in Figs. 12 and 13 show that the number of warm days and cold nights fluctuates in an increasing and decreasing trend, respectively. In addition, similarities are seen in the change trends between warm days and cold nights. For example, in 1980, under the continual influence of strong cold air in the north, low-temperature weather occurred continuously in most areas of China and many areas experienced low-temperature disasters, which led to a decrease in the number of warm days and an increase in the number of cold nights. In 2015, 2016, and 2017, the temperature continued to rise, with high temperatures that occurred once in decades. This finding is closely related to the severe El Niño events that occurred in 2015 and 2016, the impact of the subtropical high in 2017, and the overall global warming trend. From 1979 to 2018, there has also been an increase in the number of warm days and a decrease in the number of cold nights. Meteorological events can indirectly verify the accuracy of this product, indicating that the corrected data can be used to analyze long-term temporal and spatial changes in temperature.

To further analyze the change rate and regional differences in T_avg during the study period, we analyzed the temperature change rate (K), correlation coefficient (R), and significance test of the correlation coefficient (T test(R)). As shown in Fig. 14a and a.i, the T_avg in most regions of China shows a weak positive warming trend, accounting for 92.13 % of the total, and the average temperature of T_avg in each region increased by 0.03 ^∘C yr⁻¹. The analysis of R in Fig. 14b and b.i shows that they show a strong correlation of approximately 48.77 % and a correlation in the region of 84.06 %, which shows that there is a high correlation between temperature changes and time. Figure 14c and c.i show that after performing a significance test on the R between temperature and time, 83.17 % of the area passed the 95 % significance test and 75.23 % of the area passed the 99 % significance test, which shows that the correlation between temperature and time development is significant.

6 Data availability

The daily T_a products at 0.1^∘ resolution from 1979 to 2018 are freely available to the public in tif format at https://doi.org/10.5281/zenodo.5502275 (Fang et al., 2021a), and are distributed under a Creative Commons Attribution 4.0 License.

7 Code availability

The technical code of the T_a dataset based on the reconstruction model and verification can be downloaded at https://doi.org/10.5281/zenodo.5513811 (Fang et al., 2021b). We have been finishing and improving the code and plan to upload it as a supplementary version.

8 Conclusions

T_a is an indispensable variable for global climate change research. Therefore, how to obtain high-precision and high-temporal-resolution air temperature data products is an important issue. Many researchers have endeavored to produce datasets by using different data sources for the global or local region. However, because of the need for refinement of research, further improvements in accuracy and spatio-temporal resolution are necessary. Based on the full analysis of the advantages and disadvantages of various datasets and data sources, this study integrated various data sources, such as in situ data, remote-sensing data, and reanalysis data, and proposes a reconstruction model of T_a under clear-sky and non-clear-sky weather conditions. A multiple linear regression model was used to further improve the accuracy of the data, and we obtained a new set of gridded high-resolution daily temperature datasets in China from 1979 to 2018. For T_max, validation using in situ data shows that the RMSE ranges from 0.86 to 1.78^∘, the MAE varies from 0.63 to 1.40^∘, and the R² ranges from 0.96 to 0.99. For T_min, the RMSE ranges from 0.78 to 2.09^∘, the MAE varies from 0.58 to 1.61^∘, and the R² ranges from 0.95 to 0.99. For T_avg, the RMSE ranges from 0.35 to 1.00^∘, the MAE varies from 0.27 to 0.68^∘, and the R² ranges from 0.99 to 1.00. Furthermore, we verified the T_a dataset with the existing reanalysis dataset and found that the proposed dataset has credibility and accuracy. Moreover, based on the particularity of geographic climate change in different regions, we used four extreme climate indicators (TXx and TNn anomalies, TX90p, and TN10p) and three climate change indices (K, R, and T test) to analyze the trend changes of T_max, T_min, and T_avg. In summary, the temperature in most regions of China has been gradually increasing. The number of cold nights and warm days has gradually decreased and increased, respectively, and T_max and T_min have gradually increased, which is consistent with the general trend of global warming.

However, due to various factors, the weather may occasionally change drastically, such as to hail. Historical data cannot provide weather information to a greater specificity than was possible at that time; thus, particularly in areas without meteorological stations, refining past data is difficult. However, further research should consider more meteorological satellite data, especially geostationary meteorological satellite data, to improve the accuracy of surface temperature datasets used to monitor climate change.