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 scientiﬁc and policy communities but is still not available. In this paper


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(Lan et al., , 2022Zhang 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., 2015Chen et al., , 2021Prihodko 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-/mediumresolution 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(Chen et al., , 2021Hooker 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 satelliteobserved 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.

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, satellitederived LSTs, and elevation. In our algorithm, groundstation-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 (middaytime 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. Framework for implementing the SVCM-SP algorithm in a region (e.g., Africa). β 0 , β elev , and β lst are the intercept, coefficients of elevation, and LST, respectively.

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. 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 highspatial-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-stationbased 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: where both the variables (e.g., T a , elevation, and LST) and the model parameters are functions of locations/coordinates (u i v 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 middaytime 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   Land cover type   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 crossvalidation.

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, re-  spectively. 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 loweraccuracy climate zones are generally located at the boundary of regions where some directions have no weather stations.
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.

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.
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.
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.

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 (Thorn- 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  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). ton 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 . 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 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 .

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".

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.
Author contributions. YZ designed the research, TZ implemented the research and wrote the original manuscript, and YZ and ZZ supervised the research. All co-authors revised the manuscript and contributed to the writing.
Competing interests. At least one of the (co-)authors is a member of the editorial board of Earth System Science Data. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Financial support. This research was supported by the College of Liberal Arts and Science (LAS) Dean's Emerging Faculty Leaders award at the Iowa State University and the National Science Foundation (2041859, 1803920, and 2203207).
Review statement. This paper was edited by Jing Wei and reviewed by three anonymous referees.