A global dataset of daily maximum and minimum near-surface air temperature at 1 km resolution over land (2003–2020)

Near-surface air temperature (T_a) is a key variable in global climate studies. A global gridded dataset of daily maximum and minimum T_a (T_max and T_min) is particularly valuable and critically needed in the scientific and policy communities but is still not available. In this paper, we developed a global dataset of daily T_max and T_min at 1 km resolution over land across 50^∘ S–79^∘ N from 2003 to 2020 through the combined use of ground-station-based T_a measurements and satellite observations (i.e., digital elevation model and land surface temperature) via a state-of-the-art statistical method named Spatially Varying Coefficient Models with Sign Preservation (SVCM-SP). The root mean square errors in our estimates ranged from 1.20 to 2.44 ^∘C for T_max and 1.69 to 2.39 ^∘C for T_min. We found that the accuracies were affected primarily by land cover types, elevation ranges, and climate backgrounds. Our dataset correctly represents a negative relationship between T_a and elevation and a positive relationship between T_a and land surface temperature; it captured spatial and temporal patterns of T_a realistically. This global 1 km gridded daily T_max and T_min dataset is the first of its kind, and we expect it to be of great value to global studies such as the urban heat island phenomenon, hydrological modeling, and epidemic forecasting. The data have been published by Iowa State University at https://doi.org/10.25380/iastate.c.6005185 (Zhang and Zhou, 2022).

1 Introduction

Near-surface air temperature (T_a) refers to the atmospheric temperature 1.5–2 m above surfaces, and it is an important variable for numerous applications, especially those pertinent to climate and environment change (Huang et al., 2019; Zhang et al., 2018), terrestrial hydrology and phenology (Lin et al., 2012; Ren et al., 2019), public health (Lan et al., 2010, 2022; Zhang et al., 2019), disease vector propagation (Lowen et al., 2007; Petrova and Russell, 2018; Wu et al., 2020), and epidemic forecasting (Verma et al., 2017; Connor et al., 1998). T_a generally varies across space and time dramatically due to the spatial heterogeneity and temporal dynamics of environmental factors such as solar radiation, wind speed, land cover, cloud cover, and vegetation phenology (Benali et al., 2012; Chen et al., 2015, 2021; Prihodko and Goward, 1997). At the global scale, a T_a dataset will be of limited or no use if it does not characterize and capture such fine spatial details and continuous temporal coverage. A high-/medium-resolution global T_a dataset at the daily interval is highly desirable.

Many global or regional T_a datasets have been previously published (Chen et al., 2021; Crespi et al., 2021; Fang et al., 2021; Hersbach et al., 2018; Hooker et al., 2018; Kalnay et al., 1996; MacDonald et al., 2020; Meyer et al., 2019; Nashwan et al., 2019; Oyler et al., 2015; Thornton et al., 2021; Werner et al., 2019); however, these have either coarse spatiotemporal resolutions or only cover specific regions (Table S3 in the Supplement). For example, some global T_a datasets have daily frequencies but at coarse spatial resolutions (e.g., 0.05^∘ or even coarser) (Hersbach et al., 2018; Hooker et al., 2018; Kalnay et al., 1996); other T_a datasets with medium spatial resolutions (∼ 1 km) are only available for specific regions such as North America and mainland China (MacDonald et al., 2020; Oyler et al., 2015; Thornton et al., 2021; Chen et al., 2021). There are also several T_a datasets at even finer spatial resolutions but generated only for much smaller spatial regions (Crespi et al., 2021; Meyer et al., 2019; Nashwan et al., 2019; Werner et al., 2019). Despite the increasing need for a global daily T_a at finer resolutions (e.g., 1 km), such global products do not exist yet – a gap still not filled yet.

Methodologically speaking, a range of techniques have been proposed and applied to generate T_a products; the majority of them rely on combining weather station data and gridded auxiliary datasets to simply make spatial interpolation or build empirical predictive models (Chen et al., 2015; Goward et al., 1994; Hengl et al., 2012; Hou et al., 2013; Hrisko et al., 2020; Li and Zha, 2019; Li et al., 2018; Nemani and Running, 1989; Rao et al., 2019; Shen et al., 2020; Shi et al., 2017; Sun et al., 2005; Yoo et al., 2018; Zhu et al., 2013). Common spatial interpolation algorithms, such as inverse distance weighting (IDW), spline, and kriging, are unlikely to be applicable at the global scale, for example, due to the relative sparsity of weather stations and the high spatial heterogeneity of T_a (Chai et al., 2011; Dodson and Marks, 1997; Li and Heap, 2011; Stahl et al., 2006). Model-based approaches are often a better choice to capture the true spatial variability in T_a; these methods are roughly divided into three groups. The first is the temperature–vegetation index (TVX) method, which estimates T_a from the maximum normalized difference vegetation index (NDVI) based on the assumption that the canopy temperature over fully covered vegetation approximates near-surface T_a (Goward et al., 1994; Nemani and Running, 1989; Zhu et al., 2013). An apparent weakness of the TVX method is its large uncertainty or unsuitability for regions with low vegetation cover. The second group is the energy balance method, leveraging the explicit modeling of surface energy balance and the quantification of net radiation (e.g., the sum of sensible, soil, and latent heat fluxes) (Sun et al., 2005; Zhang et al., 2015). Energy-based methods are physically based, requiring detailed characterization of surface biophysical conditions and thereby making it difficult to implement for large areas due to the lack of such detailed biophysical parameters.

Of the three groups, the third category is statistical methods that estimate T_a via statistical relationships empirically calibrated between T_a and other covariates. Common algorithms used include geographically weighted regression (GWR), cubist, random forests, and deep learning, among others (Chen et al., 2015, 2021; Hooker et al., 2018; Li et al., 2018; Rao et al., 2019; Shen et al., 2020; Yoo et al., 2018). Compared to physical-based methods, statistical methods have fewer restrictions on data requirements and better applicability for large spatial extents (Noi et al., 2017). However, direct applications of common statistical methods often fail to capture and preserve relationships between T_a and auxiliary covariates (e.g., an unrealistically positive relationship between T_a and elevation), thereby leading to large uncertainties or even incorrect results. To overcome such drawbacks, we recently proposed a class of Spatially Varying Coefficient Models with Sign Preservation (SVCM-SP) (Kim et al., 2021; Zhang et al., 2022b), which can capture and preserve relationships between T_a and explanatory variables. The SVCM-SP algorithm was originally implemented for estimating T_a over mainland China, with significant improvement in terms of both accuracies and computational efficiency compared to alternative methods such as GWR (Zhang et al., 2022b). The potential of SVCM-SP as a routine algorithm to generate global T_a products is still untapped and addressed here.

Here we aim to generate a global dataset of daily maximum and minimum T_a (T_max and T_min) at 1 km resolution across 50^∘ S–79^∘ N from 2003 to 2020 by integrating ground-station-based T_a measurements and gridded satellite-observed covariates (i.e., digital elevation model, DEM; land surface temperature, LST). We employed our newly developed SVCM-SP algorithm that, for example, can preserve negative and positive relationships with elevation and LST, respectively. Zhang et al. (2022b) successfully estimated and validated gridded T_a using the SVCM-SP algorithm and demonstrated its novelty through the comparison with the GWR model, while in this study, we developed the global product of gridded T_a, performed extensive model calibration and accuracy assessment at the global scale, and provided details on accuracy and spatial and temporal patterns of the global gridded T_a. Our dataset aims to provide the first ever 1 km resolution daily maximum and minimum T_a dataset with a global coverage.

Figure 1 Regions and locations of weather stations in this study. Red points are the locations of weather stations, and polygons are the boundaries of regions used in the SVCM-SP algorithm. Specifically, polygons of red, purple, orange, blue, and black represent the boundaries of North America, South America, Africa, Australia, and Europe and Asia, respectively.

2 Study area and data

Land areas covered by our global dataset span from 50^∘ S to 79^∘ N. We divided the global lands roughly into five regions: North America, South America, Europe and Asia, Africa, and Australia and New Zealand. To encompass all the land areas resolved at 1 km resolution, as well as to cover all the possible weather stations, the boundaries of the five regions were irregular. There also exist some overlaps between the regions. Our analysis considered three major data sources: ground-station-based T_a observations, satellite-derived LSTs, and elevation. In our algorithm, ground-station-based T_a is assumed to be statistically related to satellite LST and elevation. Details about each data source are further described below.

Ground-station-based daily T_max and T_min were compiled from a total of 103 156 weather stations from 2003 to 2020. These are obtained from two climatology networks: the Global Historical Climatology Network daily (GHCNd) across the world and the China Meteorological Data Service Centre (CMDC) across mainland China. The LST dataset is a global seamless 1 km resolution LST dataset at a daily (mid-daytime and mid-nighttime) frequency from 2003 to 2020, which was gap-filled from the MODIS LST products (Zhang et al., 2022a). Both the mid-daytime and mid-nighttime LSTs were considered in our analysis. The DEM layer we used is the SRTM30_PLUS product at 1 km resolution (Becker et al., 2009), which has been generated from the combination of the Shuttle Radar Topography Mission (SRTM30) topography (collected in 2000) (Hennig et al., 2001; Rosen, 2000) within a latitude of ± 55^∘, ICESat-derived topography (collected from 1 February 2003 to 30 June 2005) (Dimarzio et al., 2007) in Antarctica, and the GTOPO30 topography (completed in late 1996) (Danielson and Gesch, 2011) in the Arctic. Besides, the Köppen–Geiger climate zones and MODIS land cover data (MCD12Q1) were acquired to divide the world into zones for accuracy assessment (Sulla-Menashe and Friedl, 2018). Specifically, our dataset covers a small portion of Greenland which is constrained by the extent of the global seamless 1 km daily LST dataset.

Figure 2 Framework for implementing the SVCM-SP algorithm in a region (e.g., Africa). β₀, β_elev, and β_lst are the intercept, coefficients of elevation, and LST, respectively.

Figure 3 Scatter plots between estimated and observed T_a in five regions in the year 2010. Each point represents the estimated and observed T_a (T_max or T_min) on a specific day at a weather station. The color of points represents the density, in which red and blue points represent the high and low densities, respectively. The red line is the regression line, and the black line is the 1:1 line.

3 Methods

The core of our methodological framework is the SVCM-SP algorithm that correlates ground-station-based T_a with satellite LST and elevation. We applied the SVCM-SP algorithm to estimate T_max and T_min separately. To capture potential non-stationarity, the algorithm was trained for each day of the period 2003–2020, as well as for each of the five regions. Before applying the SVCM-SP algorithm, weather station T_a data were first pre-processed and filtered for quality control to ensure the high fidelity of reference T_a observations.

Table 1 Multi-year average accuracies for T_max and T_min in different regions from 2003 to 2020.

Note: the selected testing stations were within 50 km surrounding the training stations. SD represents the corresponding standard deviation.

Table 2 Model performance for different land cover types in 2003–2020.

Note: SD represents the corresponding standard deviation.

Table 3 Model performance for different elevation ranges in 2003–2020.

Note: SD represents the corresponding standard deviation.

More specially, we first processed the weather station data in three ways. First, the locations of many weather stations in China, especially those located in complex terrains, are not accurately documented, geo-referenced only at the level of arc degrees and minutes in the metadata. Such location errors have to be corrected, and we manually corrected the locations of those weather stations located over complex terrains by searching the meteorological observation fields near the reported locations of weather stations with the help of high-spatial-resolution images from ArcGIS base map or Google Maps (Zhang et al., 2022b). Second, there are missing values, especially in stations in Africa and South America (Fig. S1 in the Supplement). We filled these data gaps using a 5 d local moving window (Fig. S2 in the Supplement). Accordingly, the number of records largely increased (Figs. S3 and S4 in the Supplement) with reasonable error ranges (Fig. S5 in the Supplement). Third, the processed ground-station-based T_a data from the two steps were overlaid and matched with satellite LST and elevation to extract pairs of ground-station-based T_a and satellite covariates as inputs to the SVCM-SP algorithm. Specifically, mid-daytime and mid-nighttime LSTs were used to develop their relationship with air temperature to interpolate station T_max and T_min, respectively. The actual time of T_max and T_min may be slightly different from mid-daytime and mid-nighttime LSTs. Within the small difference in time between LST and $T_{max}/T_{min}$, there will not be significant change in the spatial variations in LST. Therefore, the impact of time difference between LST and $T_{max}/T_{min}$ on the accuracy of the estimated T_a is minor as shown by shifting LST for time difference (Fig. S8 in the Supplement).

The SVCM-SP algorithm seeks to build a spatially varying relationship between ground-station-based T_a with LST and elevation with sign preservation (Kim et al., 2021; Zhang et al., 2022b). A salient feature distinguishing it from conventional regression approaches is the spatially varying nature with constraints of estimated coefficients in the predictive relationship:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{T}_{\mathrm{a}}(u_i,v_i)=& {\mathit{\beta }}_{\mathrm{0}}(u_i,v_i)+{\mathit{\beta }}_{\text{elev}}(u_i,v_i)\\ & \times \text{ Elevation }(u_i,v_i)+{\mathit{\beta }}_{\text{lst}}(u_i,v_i)\\ & \times \text{ LST }(u_i,v_i)+{\mathit{\epsilon }}_i,\end{array}\end{array}$$

where both the variables (e.g., T_a, elevation, and LST) and the model parameters are functions of locations/coordinates (u_iv_i). More importantly, the two slope parameters (e.g., β_elev and β_lst) are constrained to be negative and positive, respectively. ε_i is the normal random error with mean zero and finite variance. These unknown parameters were estimated with a penalized bivariate spline method based on the triangulation technique under constraints. Details about the SVCM-SP algorithm are reported in Kim et al. (2021) and Zhang et al. (2022b). To estimate T_a across the globe, we applied the SVCM-SP algorithm to develop region-specific relationships for the five regions (Fig. 2). Also, two separate sets of equations were developed, one for T_max using mid-daytime LST as the explanatory variable, and another for T_min using midnight LST as the explanatory variable. The model performance for estimating gridded T_a was assessed based on root mean square error (RMSE) and mean square error (MAE) using the 10-fold cross-validation in these regions in each day. Taking the RMSE as an example, a RMSE was generated in each test of the 10-fold cross-validation, and all RMSEs from the 10 tests were averaged as the final RMSE on a specific day in a specific region. This accuracy assessment using the 10-fold cross-validation was implemented based on independent validation data and can provide a reliable evaluation of the accuracy. For each station, we can also calculate RMSE based on the time series of estimated and validation T_a from the 10-fold cross-validation. Accordingly, we can calculate mean RMSE and corresponding standard deviation in each land cover type, climate type, and elevation range. Specifically, this accuracy assessment represents conservative estimates of the uncertainties of our data because when producing the final results, we used all the available data, more than those in the 10-fold cross-validation.

Table 4 Model performance for different climate types in 2003–2020.

Note: SD represents the corresponding standard deviation.

Figure 4 Accuracy of estimated T_a in climate zones in 2003–2020. Climate zones with black boundaries are areas with low densities of weather stations (i.e., distances between training and validation sites are larger than 50 km). The white regions on land are areas without reliable evaluations due to the lack of weather stations.

4 Results and discussion

4.1 Accuracy of the estimated T_a

The results of the 10-fold cross-validation indicate the accuracy of estimated T_a varies across regions within a reasonable range (Fig. 3 and Table 1). The estimated and observed T_a in different regions scattered along the 1:1 line with the RMSE ranging from 1.17 to 2.38 and 1.59 to 2.34 ^∘C, respectively, for T_max and T_min in 2010 (Fig. 3). As shown in Table 1, the estimated average RMSE and MAE from 2003 to 2020 ranged from 1.20 to 2.44 and 0.89 to 1.82 ^∘C, respectively. The highest accuracy was obtained in Australia for T_max, with a RMSE and MAE of 1.20 and 0.89 ^∘C, respectively. The lowest accuracy was obtained in North America for T_max, with a RMSE and MAE of 2.44 and 1.82 ^∘C, respectively. Meanwhile, the variation in accuracy across years in each region is smaller compared to spatial variations in the accuracy across regions (Tables S1 and S2 in the Supplement). The variations in accuracy may be caused by the differences in climate and topography in these regions (Hooker et al., 2018). For example, Australia is a continent with the gentlest undulations of terrains with about 87 % of the land below 500 m a.s.l. and is dominated by hot arid desert and steppe climates, leading to the smallest spatial variations in T_a. However, other regions contain a variety of dominant climate types and geomorphic types, contributing to the large spatial variability observed in T_a.

The accuracy of estimated T_a varied in different land cover types, elevation ranges, and climate types (Tables 2–4). RMSE and MAE for T_max ranged from 2.06 to 2.56 and 1.54 to 1.97 ^∘C, respectively, and these indicators for T_min range from 1.84 to 2.83 and 1.40 to 1.96 ^∘C, respectively. The model performs well for an impervious surface (with the lowest RMSE), cropland, water, and wetland, whereas RMSE values were higher for the tundra and bare land, which was generally consistent with the findings of existing studies in mainland China (Chen et al., 2021; Shen et al., 2020; Zhang et al., 2022b). As shown in Table 3, RMSE and MAE values vary with elevation ranges but did not increase with the increase in elevation ranges, which is different from existing findings (Chen et al., 2015; Rao et al., 2019). This is because we only used weather stations within the distance of 50 km from the training sites to evaluate the accuracy of estimated T_a in this study, which can mitigate the effects of sparse weather stations at high elevations on accuracy assessment, as reported in existing studies (Chen et al., 2015; Rao et al., 2019). RMSE and MAE values in equatorial climate zones are distinctly lower than those of other climate zones (Table 4), indicating the highest accuracies for both T_max and T_min possibly due to T_a near the Equator being generally warmer and less intra-annual variations compared to other climate zones (Legates and Willmott, 1990).

Spatial distributions of RMSE illustrate that most of the climate zones show reasonable accuracies (RMSE < 3.0 ^∘C) for T_max and T_min (Fig. 4). The lower-accuracy climate zones (RMSE > 3.0 ^∘C) mainly occur where there are low station densities (Fig. S6 in the Supplement), which is consistent with the finding of decreasing accuracy with the increase in station density (Shen et al., 2020). Meanwhile, these lower-accuracy climate zones are generally located at the boundary of regions where some directions have no weather stations.

Figure 5 Temporal patterns of accuracies in estimated T_a in different regions in the year 2010.

The RMSE values generally show distinctly seasonal patterns in the five regions within reasonable ranges (Fig. 5). Taking the year 2010 as an example, RMSEs in summer (June, July, and August) are generally lower than those in winter (December, January, and February) in North American, European, and Asian regions (Fig. 5) possibly due to plant phenology, which leads to a closer relationship between T_a and LST in the summer than in the winter (Benali et al., 2012; Cai et al., 2017; Lin et al., 2012). This seasonal variation is less obvious in Africa and South America possibly due to weaker correlations between plant phenology and air temperature in the two regions, which are located across the Equator (Adole et al., 2019; Sakai and Kitajima, 2019). Specifically, the Australian region, which is located in the Southern Hemisphere, shows higher RMSEs in summer (December, January, and February) than in winter (June, July, and August) for T_max. This may be caused by more homogeneous spatial variations in T_max in winter than in summer in the Australian region.

Figure 6 Spatial pattern of estimated T_a at the global scale on an example day, day 200, in 2010.

4.2 Spatial and temporal patterns of T_a

The estimated T_a shows significant spatial variations at the global scale (Fig. 6). Taking the estimated T_a on one July day as an example, both T_max and T_min decrease from about 30^∘ N to the North and South poles (Fig. 6). Meanwhile, lower T_a values also occur at higher-elevation regions such as the Tibetan Plateau in the center of Asia and the Andes Mountains in the west of South America. Therefore, the characteristics of T_a change with latitude and elevation (i.e., the trend of lower T_a in higher-latitude/elevation areas), which is consistent with the existing studies (Chen et al., 2015; Zhang et al., 2022b). The highest T_a values occur in northern Africa and the Arabian Peninsula, as these regions are mainly covered by the Gobi Desert.

Figure 7 Spatial pattern of estimated T_a in five representative cities on day 200 of the year 2010. A city shape includes the urban T_a extracted by using nighttime light data (Zhou et al., 2018) and a surrounding buffer of equal size. Black polygons are the urban extents extracted by using global 30 m resolution artificial impervious area data (Li et al., 2020).

The spatial patterns of estimated T_a in selected cities with clear weather around the world illustrate that the urban heat island (UHI) phenomenon (i.e., the higher temperature in urban than in the surrounding rural areas) has been well captured at the city scale (Fig. 7). On an example day of July in 2010, the estimated T_a in these cities shows an obvious UHI phenomenon, which is reasonable with the transition from urban centers to surrounding rural areas. The estimated T_a in Changsha, China, shows several hotspots because some nearby cities (such as Xiangtan and Zhuzhou) have also been included in the buffer of Changsha, indicating the effectiveness of the estimated T_a for presenting UHI in small urban areas. Specifically, as a coastal city, estimated T_a in Melbourne, Australia, shows decreasing trends from the coast, and the UHI phenomenon is not obvious in surrounding small cities. This is because there is also an increasing trend of elevation from the coast in Melbourne, leading to the mixed spatial patterns of T_a due to the UHI phenomenon and elevation changes.

Figure 8 Temporal pattern of estimated and observed T_a at the weather station of Beijing (black point) in the year 2010. The black rectangles are example days showing maps of estimated data in Beijing.

The comparison of the temporal pattern between estimated T_a and ground-station-based measurements from an example of weather stations in a mega-city (Fig. 8) illustrates that the SVCM-SP algorithm can effectively (RMSE of 1.25 ^∘C and 1.53 ^∘C, respectively, for T_max and T_min) estimate T_a for the entire period. As shown in Fig. 8, the estimated T_a based on 10-fold cross-validation and T_a observations from the weather station in Beijing, China, show similar temporal patterns and very close values for both T_max and T_min in 2010. For both clear weather (days 28 and 130 in Fig. 8) and overcast weather (days 219 and 293 in Fig. 8) (Zhang et al., 2022a), the gridded T_a can illustrate the UHI phenomenon. An existing study has found that the estimated T_a in urban areas was more accurate than those of other regions (Zhang et al., 2022b), specifically suggesting its great value for urban applications.

4.3 Comparison with existing T_a datasets

The gridded T_a data in this study have advantages regarding spatiotemporal resolutions (i.e., 1 km and daily maximum and minimum) and its global coverage (Table S3). The spatial resolution of existing global T_a datasets with daily frequencies and long-term coverage is generally low (e.g., 0.25^∘) (Hersbach et al., 2018; Kalnay et al., 1996). T_a datasets with improved spatial resolutions (e.g., 1 km) are usually only available at the continental or national scales (Chen et al., 2021; Fang et al., 2021; MacDonald et al., 2020; Oyler et al., 2015; Thornton et al., 2021).

The gridded T_a in this study can effectively capture the spatial variation in T_a by preserving physical relationships between T_a and response variables (Fig. S7 in the Supplement). In other T_a datasets, such physical relationships (e.g., positive relationship between T_a and LST) cannot always be preserved in some situations because these datasets were created using methods without explicit constraints on the relationships between T_a and response variables. Efforts have been made to build vertical lapse models to estimate gridded T_a according to adiabatic lapse rate (ALR) (Dodson and Marks, 1997; Rhee and Im, 2014; Thornton et al., 2021; Zhu et al., 2017), but the generalization of these models is limited because it is difficult to accurately capture ALR due to its spatial change.

The accuracy of the resulting gridded T_a from this study is comparable to several other reported gridded T_a datasets (e.g., Chen et al., 2021; Oyler et al., 2015; Thornton et al., 2021). Among them, the 1 km daily T_a from Daymet (Thornton et al., 2021) reaches an MAE of 1.52 and 1.78 ^∘C for T_max and T_min, respectively, and the 30 arcsec (∼ 800 m) daily T_a from TopoWx (Oyler et al., 2015) reaches 1.03 and 1.06 ^∘C, while in this study, the average MAE is 1.82 and 1.78 ^∘C in North America. However, Daymet failed to capture the UHI phenomenon due to the spatial interpolation of T_a being implemented based only on elevation (Menne et al., 2012) and did not consider the impact of biophysical and socioeconomic factors on spatial variations in T_a (Li et al., 2018). Therefore, Daymet has difficulties in capturing the spatial variation in T_a in urban areas, although its accuracy is comparable to our dataset. The estimated T_a from TopoWx can display the UHI phenomenon but tends to overestimate the impact of topographical features and show fewer temporal variations in the spatial pattern of T_a within a month than that in this study, as a 10-year average of monthly LSTs was used as a covariate in TopoWx (Li et al., 2018; Oyler et al., 2015) instead of daily LST data in this study. The 1 km daily average T_a data by Chen et al. (2021) reach a RMSE of 1.615 to 1.957 K using leave-location-out cross-validation in mainland China, while the average RMSE of estimated T_max and T_min is 1.80 and 1.75 ^∘C, respectively, in Europe and Asia. While the accuracy of T_a obtained in this study is comparable to the other large-scale T_a datasets, our dataset is produced at the global scale using consistent modeling and assessment approaches.

There are some limitations in the SVCM-SP algorithm used in this study for creating the gridded T_a dataset, and future work can focus on improving the accuracy of the estimated T_a with an improved SVCM-SP algorithm. First, we only considered the linear relationship between T_a and covariates. However, nonlinear relationships may exist between T_a with elevation and LST when other factors, such as winds, clouds, snow, and land cover types, have non-negligible impacts on T_a (Cai et al., 2017; Good, 2016). Second, we only used two covariates in the SVCM-SP algorithm, although the potential of generalization of our framework is large. Additional covariates (e.g., other surface characters such as Geoscience Laser Altimeter System (GLAS)-derived canopy height and vegetation parameters) can be explored in the SVCM-SP algorithm to further improve the model performance. Third, the limited number of valid station observations on specific days might introduce larger uncertainties in interpolating T_a using the SVCM-SP algorithm. In future studies, station observations from neighboring days can be explored to improve the interpolation of T_a.

5 Data availability

Data described in this paper can be accessed at Iowa State University's DataShare at https://doi.org/10.25380/iastate.c.6005185 (Zhang and Zhou, 2022). The dataset contains 36 sub-datasets. Each sub-dataset contains T_max or T_min of each year from 2003 to 2020 in five regions (i.e., North America, South America, Europe and Asia, Africa, and Australia and New Zealand). The data are in GeoTIFF with the georeferenced information embedded. The MODIS ellipse sinusoidal projection with a spatial resolution of 1 km is used in the data. The unit of LST in GeoTIFF is 0.1^∘ Celsius (^∘C), and the naming rule can be found in the file “README.pdf”.

6 Conclusions

We generated a global land (50^∘ S–79^∘ N) 1 km daily maximum and minimum T_a (i.e., T_max and T_min) dataset from 2003 to 2020 based on ground-station-based T_a measurements and gap-filled LST dataset using the Spatially Varying Coefficient Models with Sign Preservation (SVCM-SP) algorithm. The dataset showed acceptable accuracies based on the 10-fold cross-validation for five regions of the globe compared to existing T_a datasets. The RMSEs of estimated T_max and T_min ranged from 1.20 to 2.44 and 1.69 to 2.39 ^∘C, respectively. The estimated T_a was affected by land cover types, elevation ranges, and climate types, with varying accuracies but within reasonable ranges. Our gridded T_a dataset effectively captured the spatial variation in T_a under clear physical meanings (i.e., negative and positive relationships with elevation and LST, respectively), which is not always true in other gridded T_a datasets. The new dataset is unique in terms of spatiotemporal resolutions (i.e., 1 km daily maximum and minimum), global coverage, and temporal span and should be useful for a wide range of applications such as the urban heat island phenomenon, hydrological modeling, and epidemic forecasting. Future work can focus on improving the accuracy of the gridded T_a dataset using the SVCM-SP algorithm by exploring more explanatory variables which are available over large areas.