The application of trivariate thin-plate smoothing splines to the interpolation of daily weather data is investigated. The method was used to develop spatial models of daily minimum and maximum temperature and daily precipitation for all of Canada, at a spatial resolution of 300 arc s of latitude and longitude, for the period 1961–2003. Each daily model was optimized automatically by minimizing the generalized cross validation. The fitted trivariate splines incorporated a spatially varying dependence on ground elevation and were able to adapt automatically to the large variation in station density over Canada. Extensive quality control measures were performed on the source data. Error estimates for the fitted surfaces based on withheld data across southern Canada were comparable to, or smaller than, errors obtained by daily interpolation studies elsewhere with denser data networks. Mean absolute errors in daily maximum and minimum temperature averaged over all years were 1.1° and 1.6°C, respectively. Daily temperature extremes were also well matched. Daily precipitation is challenging because of short correlation length scales, the preponderance of zeros, and significant error associated with measurement of snow. A two-stage approach was adopted in which precipitation occurrence was estimated and then used in conjunction with a surface of positive precipitation values. Daily precipitation occurrence was correctly predicted 83% of the time. Withheld errors in daily precipitation were small, with mean absolute errors of 2.9 mm, although these were relatively large in percentage terms. However, mean percent absolute errors in seasonal and annual precipitation totals were 14% and 9%, respectively, and seasonal precipitation upper 95th percentiles were attenuated on average by 8%. Precipitation and daily maximum temperatures were most accurately interpolated in the autumn, consistent with the large well-organized synoptic systems that prevail in this season. Daily minimum temperatures were most accurately interpolated in summer. The withheld data tests indicate that the models can be used with confidence across southern Canada in applications that depend on daily temperature and accumulated seasonal and annual precipitation. They should be used with care in applications that depend critically on daily precipitation extremes.
Spatially continuous models of daily weather variables are increasingly in demand to support many applications. Daily inputs are needed for water budget modeling (Akinremi et al. 1996), extreme event models (Wagner 1996), analysis of crop and tree growth (Easterling and Apps 2005), and modeling the occurrence of insect pests and disease (Boland et al. 2004). Crop yield can be significantly influenced by daily temperature extremes, such as frosts, which are masked by monthly climate averages (Lobell et al. 2007). Similarly, the water balance of crops is sensitive, not only to the total amount of precipitation received over a period of time, but also to the finer-scale temporal patterns of precipitation events (Shen et al. 2001).
Spatial interpolation of daily weather offers a number of challenges beyond those encountered with long-term monthly mean climatologies. Daily models are required to track complex spatial patterns due to elevation gradients, weather fronts, vegetation cover, and local exposure to water bodies. Daily models also require a large amount of computing time and storage space. Further requirements are operational simplicity and reliability of the model optimization procedures that need to be applied systematically over many days. Thus, despite the demand for daily data, there are relatively few high-resolution models that provide extensive finescale spatiotemporal coverage of both daily temperature and daily precipitation.
These considerations have led Thornton et al. (1997) to devise a local, elevation detrended, distance-weighted method to interpolate daily meteorological variables over complex terrain. Thornton (2007) used this procedure to generate surfaces of daily meteorological variables for the period 1980–97 over the conterminous United States and Hasenauer et al. (2003) applied the same method to a dense network over Austria. Such analyses are a significant advance over earlier coarse-resolution studies using relatively simple nearest-neighbor and Cressman methods (Piper and Stewart 1996; Higgins et al. 2000). Other recent studies with more restricted spatial coverage include Hunter and Meentemeyer (2005) who combined ordinary kriging with long-term monthly mean climate maps to interpolate 2-km grids of daily temperature and precipitation over California for the period 1980–2003 and Shen et al. (2001) who used a hybrid inverse distance-weighted and nearest-neighbor method to interpolate daily temperature and precipitation values onto ecodistrict polygons in Alberta for the period 1961–97.
Here we apply trivariate thin-plate smoothing splines, as implemented in the “ANUSPLIN” software (Hutchinson 1995a, 2004), to model the complex spatial patterns associated with daily data across Canada as spatially continuous functions of longitude, latitude, and elevation. These models are required particularly for agricultural and forestry applications in the southern half of Canada. Thin-plate splines have been found to be superior to the gradient plus inverse distance squared (GIDS) method of Nalder and Wein (1998) in interpolating monthly mean temperature and precipitation over Canada (Price et al. 2000). They have been less commonly applied to daily data, although separate comparative daily studies by Xia et al. (2001) and Jarvis and Stuart (2001) found trivariate thin-plate spline models to yield the most accurate results.
Thin-plate spline smoothing splines can be viewed as a generalization of multivariate linear regression in which a parametric model is replaced by a smooth nonparametric function (Wahba 1990). The method optimizes the amount of data smoothing to minimize the predictive error, as measured by the generalized cross validation (GCV). This is performed automatically and simultaneously as each surface is fitted, making it well suited to the application here to many days of data. The method is robust to varying underlying spatial statistical models of the data (Hutchinson 1993) and adapts automatically to the spatial density of the data. The latter feature is particularly relevant for the large variation in density displayed by the Canadian meteorological network. The method also supports the calculation of spatially distributed standard errors (Wahba 1990; Hutchinson 1993, 1995a).
The usual approach is to fit smoothing spline functions of three independent variables representing the location of each station: longitude (in decimal degrees), latitude (in decimal degrees), and elevation above sea level (in kilometers). This relative scaling of elevation has been optimized in precipitation analyses at various time scales and locations to obtain surfaces with minimum error (Hutchinson 1995a, 1998b). It represents the impact of elevation on surface climate as ∼100 times greater than the impact of horizontal position. This is consistent with the ratio between the commonly accepted synoptic horizontal and vertical distance scales of 1000 km and 10 km, respectively (Daley 1991). Thus, temperature and precipitation patterns usually reflect a topographic landscape that is exaggerated in the vertical, leading to significant influence on precipitation patterns by relatively modest topographic features (Barros and Kuligowski 1998). The method is extended here for daily precipitation to a two-step procedure in which a binary surface of precipitation occurrence is first generated and then combined with a surface of positive precipitation values. This accommodates the differing spatial coherence that has been observed in precipitation occurrence and intensity data (Hutchinson 1995b). It is similar to the strategy employed by Thornton et al. (1997) who also separately interpolated precipitation occurrence and positive precipitation values.
The accuracy of the fitted daily models is assessed by examining statistics provided by the thin-plate smoothing spline procedure and by withholding from the analyses 50 stations broadly representing the southern half of Canada. This is where most agricultural and forestry activity occurs and where the overwhelming majority of Canadian meteorological stations are located. Comprehensive error statistics for the withheld stations are presented, including assessments of bias, mean absolute error, root-mean-square error, and daily and seasonal extremes. These error statistics are compared with corresponding error statistics presented by existing daily interpolation studies, bearing in mind the differing data densities for the different studies. A comparative assessment of accuracy in the northern half of Canada is also provided. The paper concludes by discussing the limitations of the current methodology, particularly in the representation of daily precipitation extremes, and briefly discusses possible improvements.
2. Data preparation
Considerable effort went into preparing and checking the data (Hopkinson 2005). Daily records of precipitation, and maximum and minimum temperature from Environment Canada were examined for accuracy in geographic position (i.e., latitude, longitude, and elevation), and for the presence of flags indicating uncertainty in the associated data values. Modifications were made to data values to address missing values, estimated values, trace precipitation amounts, and precipitation accumulated over multiple days. These efforts substantially increased the number of station days available and the final product is the most complete daily dataset currently available for Canada for the period.
The period 1961–2003 was selected for the current work. The number of precipitation stations that were active in any given year over this period ranged from 2000 to 3000 while the number of temperature stations varied from about 1500 to 2200, as shown in Fig. 1. For both temperature and precipitation the period of maximum station coverage extends from the early 1970s to the early 1990s. Figure 2 shows the stations used for modeling on yearday 250 in 1975. Coverage is relatively dense in the southern portion of the country but becomes very sparse in the north. In fact, around 95% of Canadian weather stations lie south of the line shown in Fig. 2 defined by
where latitude and longitude are measured in degrees. Approximately 50% of Canada lies south of this line, corresponding to a land area of about 5.0 × 106 km2. The density of station coverage in this southern half of the country is thus around 1 station per 2500 km2. This is sparser than the data network densities for all non-Canadian interpolation studies cited in this paper.
An examination of geographic coordinates of the data revealed a number of issues. There were stations with missing coordinates and some with default coordinates such as elevations of 0 m. Gross errors in latitude, longitude, and elevation were identified and corrected. These errors were readily identified as large outliers from preliminary smoothing spline analyses. For all stations with missing elevation, and some stations with zero elevation, the digital elevation model developed by Great Lakes Forestry Centre of Natural Resources Canada (more information is available online at http://www.glfc.cfs.nrcan.gc.ca/landscape/topographic_models_e.html; see also Lawrence et al. 2008) was used to estimate station elevations.
The daily temperature and precipitation data had various flags indicating missing and estimated values. The estimated values were accepted as valid data via a standard quality control process. This included comparisons with observations from surrounding stations and interpolation between preceding and subsequent observations at the same station. The precipitation data had flags indicating trace values and accumulated values over multiple days. In the Canadian data archive, trace precipitation values have been assigned a value of zero. Mekis and Hogg (1999) suggested a value of 0.1 mm for rainfall traces, 0.07 mm for snowfall traces, and 0.03 mm for ice crystal traces (very cold temperatures). Coincident precipitation, weather, and temperature observations were not available at many climate stations so a simple method, based on Mekis and Hogg (1999), was used to estimate precipitation associated with trace observations. The method used latitude and maximum temperature at each station to assign a water equivalent ranging from 0.07 mm at lower latitudes to 0.03 mm at higher latitudes and/or very cold temperatures.
The temporal distribution of accumulated precipitation was estimated from the temporal distribution of daily precipitation at neighboring stations within a radius of 100 km. The total precipitation at each of the neighboring stations was determined for the accumulation period and then the daily amounts expressed as a fraction of the total. An inverse distance-squared-weighting scheme (Cressman 1959) was used to estimate the daily fraction at the accumulation station. If there were no stations with valid data within 100 km, then values for all days in the accumulated event were set to missing. Accumulation events longer than 4 days were set to missing because of the possible impact of evaporation from the gauge. This approach enabled the use of thousands of precipitation values that would otherwise have been set to missing.
One issue that was not addressed was the difference in the climate day between ordinary climate stations and principal meteorological stations. Since July 1961 all principal stations have a climate day ending at 0600 UTC (around 0000 LST). On the other hand, ordinary climate stations, which make up the overwhelming majority of stations, usually end their day for maximum temperature and precipitation around 0700 or 0800 LST and for minimum temperature around 1700 LST (Hopkinson 2005). These time differences can result in occasional large differences in observed minimum and maximum temperature and precipitation, leading to negative biases in winter daily minimum temperatures recorded at principal stations. Work is under way to adjust daily data at principal stations to be consistent with the recording times of ordinary stations. These data will be included in subsequent studies.
a. Thin-plate smoothing spline models of surface climate
A comprehensive introduction to the technique of thin-plate smoothing splines is given in Wahba (1990) and early applications to the spatial interpolation of monthly mean climate are given in Hutchinson (1991, 1995a). Associated statistical analyses and comparisons with kriging are given in Hutchinson (1993) and Hutchinson and Gessler (1994). The basic partial spline model for N observed data values zi is given by
where each xi is a d-dimensional vector of independent variables, f is an unknown smooth function of the xi, each yi is a p-dimensional vector of independent covariates, b is an unknown p-dimensional vector of coefficients of the yi, and each ɛi is an independent, zero mean error term. The model reduces to an ordinary thin-plate spline model when there are no covariates (p = 0). In the applications here, ordinary splines are considered, with the xi representing the three coordinates longitude, longitude, and appropriately scaled elevation. The function f and the coefficient vector b are determined by minimizing
where Jm( f ) is a nonnegative measure of the complexity of f, the “roughness penalty” defined in terms of an integral of mth-order partial derivatives of f, and ρ is a positive number called the smoothing parameter. The value of the smoothing parameter is normally determined by minimizing the GCV. Operational details of the ANUSPLIN package can be found in Hutchinson (2004).
The method includes a computationally efficient implementation for larger datasets. This is based on restricting the complexity of the fitted surfaces by using a set of “knots” that are chosen to equisample the independent variable space (Wahba 1990; Luo and Wahba 1997; Hutchinson 2004). Since the number of computations is proportional to the cube of the number of knots, this can lead to a 10-fold or more reduction in computation time. An alternative computational simplification using regression splines has been described by Wood (2003).
A regular grid of estimates can be generated from the fitted thin-plate spline surface coefficients, provided a regular grid of values for the third independent variable of elevation is supplied. Thus, the resolution of the final climate grid is determined by the resolution of the supplied elevation grid. This involves a trade-off between managing the resultant data files and the variation–error in the map output. In the present study, the independent variable was provided from a digital elevation model (DEM) with 300 arc s resolution (∼10 km), although finer resolutions are possible. This DEM was derived from a Canadian DEM generated from Canada’s National Topographic Series 1:250 000 scale data using the “ANUDEM” program (Hutchinson 1989, 2008; see also Lawrence et al. 2008).
b. Modeling daily maximum and daily minimum temperature
The daily maximum and daily minimum temperature surfaces were calculated by fitting standard trivariate smoothing splines with every data point a knot. These surfaces effectively fitted spatially varying, but regionally defined, elevation lapse rates. This approach is indicated by the analysis of Bolstad et al. (1998) who found regional regression models superior to bivariate kriging models in interpolating daily temperature. Rolland (2003) similarly found significant spatial and seasonal variation in lapse rates of monthly temperature in alpine regions. The chosen model is also indicated by the analysis of monthly mean temperature by McKenney et al. (2001) who found trivariate splines to marginally outperform partial thin-plate spline analyses where the elevation dependence was restricted to be a linear dependence specified by a fitted constant value.
c. Modeling daily precipitation
On average over 50% of days in Canada receive no precipitation. A two-stage approach was adapted to model precipitation occurrence and positive precipitation separately. As noted above, this approach accommodates the differing spatial coherence that has been observed in precipitation occurrence and intensity data (Hutchinson 1995b). It also overcomes instabilities encountered in applying thin-plate splines directly to all zero and nonzero data simultaneously. We found greater spatial coherence, and hence longer spatial correlation scales, in the occurrence data than in the positive precipitation data. Thus, the occurrence analyses could be effectively performed by a thin-plate spline with simpler structure defined by 1000 knots. This was applied to all 2000–3000 precipitation observations after first converting all positive values to 1 and leaving zero values unchanged. The interpolated surface then indicated that precipitation occurred where the interpolated value exceeded a threshold of 0.5 and that precipitation did not occur where the interpolated value was less than 0.5. The final interpolated precipitation values were set to zero where the interpolated occurrence surface was less than or equal to 0.5 and set to the values from the positive precipitation surface, described below, where the occurrence surface was greater than 0.5.
The positive precipitation surface was generated from the relatively smaller number of positive precipitation values. These surfaces were fitted by a standard thin-plate spline with every data point a knot. Importantly, the surfaces were fitted to the square roots of the positive data values. This is consistent with the observation that the square root of daily rainfall is distributed in time as approximately the upper tail of a normal distribution (Hutchinson 1995b). The transport processes associated with daily precipitation imply that precipitation data are similarly distributed over the data network. The effect of using the square root transformation is to apply more smoothing to large precipitation values and less smoothing to small precipitation data values. The interpolated values are corrected for the small negative bias that the square root transformation introduces (Hutchinson 2004). Though commonly considered, a logarithmic transformation applied to precipitation data leads to spline surfaces that oversmooth larger values and, unlike the square root transformation, is undefined when applied to zero. Hutchinson (1998a) found that applying the square root transformation can reduce interpolation error by about 10%. Tait et al. (2006) have confirmed that the square root transformation can yield a significant reduction in daily precipitation interpolation error. By equilibrating the variability of daily precipitation, the square root transformation also leads to efficient detection of precipitation data errors, particularly false zeros, which give rise to significant outliers from the fitted surface. False zeros are quite common in precipitation data but tend not to give rise to large residuals if the data are not transformed.
d. Assessing the daily climate surfaces
The accuracy of the fitted daily climate surfaces was assessed in two ways: by examining a range of standard diagnostic statistics provided by the surface fitting procedure (Hutchinson 2004) and by examining residuals from stations withheld from the analyses.
The primary surface diagnostic is the GCV. This is a measure of the predictive error of the fitted surface that is calculated by implicitly removing each data point and summing with appropriate weighting the square of the difference of each omitted data value from the spline fitted to all other data points. The GCV is a robust measure of the predictive error of the spline surface, provided there is no short-range correlation in the data. The square root of the GCV can be directly compared with the root-mean-square residual from withheld data.
The signal, given by the trace of the influence matrix (Wahba 1990), is a second informative surface diagnostic. It is an explicit measure of the complexity of the fitted surface that varies between a small positive integer and the number of stations used to generate the surface. Hutchinson and Gessler (1994) suggest that the signal should be no greater than about half the number of data points. Models satisfying this condition tend to be more robust and reliable in data-sparse regions. Higher signals can indicate that the climate field being analyzed is too complex to be adequately represented by the data or that there is significant short-range correlation in the data. The latter shortcoming can be addressed by simply removing closely spaced data points from the analysis (Hutchinson and Gessler 1994; Hutchinson 1998a) or by explicitly including short-range correlation structure in the model (Wang 1998).
e. Withheld data
A set of 50 high-quality stations were withheld from the data in the southern half of Canada defined by Eq. (1). As noted above, this is where data are reasonably dense and where agricultural and forestry applications are needed. Of 368 reference climate stations available for Canada, there were 150 stations in the southern half of Canada with near-complete data coverage for the years 1961–90. The “SELNOT” program in the ANUSPLIN package was used to choose from these stations 50 stations that equisampled three-dimensional longitude, latitude, and elevation space, with longitude and latitude in degrees and elevation in kilometers. These withheld stations had virtually complete daily records and evenly sampled the full range of latitude, longitude, and elevation in the station data south of the line represented by Eq. (1). This method was preferred to a random selection process that would bias the withheld dataset to areas with relatively more stations at lower elevations. The locations of the withheld stations are shown in Fig. 3. Of the 50 withheld stations, 36 were principal climate stations, leading to slightly inflated estimates of interpolation error. The withheld stations were returned to the full dataset to construct the final models.
An important consideration when testing the accuracy of precipitation occurrence patterns is delineating wet days from dry. Inaccuracies in observing small precipitation amounts are common so standard meteorological practice was adopted in defining a wet day to be one where at least 0.2 mm of precipitation was recorded. Wet days for the interpolated precipitation surfaces were defined similarly.
The final models yielded grids of minimum and maximum temperature and precipitation for each day from 1 January 1961 to 31 December 2003—a total of 46 989 grids. The grids include the area covering 40°–84°N, 50°–153°W with a cell size of 300 arc s. These grids are available through the Canadian Forest Service (contact authors) and also through Agriculture and Agri-Food Canada (2007). Point estimates and finer-resolution grids for specific locations or applications may be generated upon special request.
a. Model assessment—Surface diagnostics
Diagnostic measures for the daily temperature surfaces are summarized by season in Table 1 and over all years in Figs. 4a,b. The ratios of the signal to the number of data points varied between about 0.3 and 0.8 and were consistent across seasons and over years. This indicates that the analyses were generally stable and robust, since ratios in this range are in reasonable agreement with the maximum value of 0.5 recommended by Hutchinson and Gessler (1994). Square root GCV values for maximum temperature averaged around 1.6°C and were generally less than 2.0°C. Square root GCV values for minimum temperature averaged around 2.0°C and were generally less than 2.5°C. These values showed a slight reduction over the first 15 yr of the analyses in line with the increase in temperature station numbers from less than 1600 to over 2000. The GCV values did not increase over the final 10 yr of analyses as station numbers declined, suggesting that data quality may have improved slightly over time. This is supported by the slight increase in temperature signal ratios after around 1980. The smaller GCV values and larger signal ratios for the maximum temperature surfaces reflect the simpler processes controlling maximum temperature. Root mean GCV values for maximum temperature were smallest in autumn and summer, and for minimum temperature were smallest in summer. Root mean GCV values for both maximum and minimum temperatures were largest in winter.
Diagnostic measures for the daily precipitation surfaces are summarized by season in Table 2 and over all years in Figs. 4c,d. As for the temperature surfaces, signal ratios for precipitation were generally stable across seasons and over years. The signal ratio for positive precipitation was generally less than 0.5 and showed no response to changes in station numbers over time. On the other hand, the signal ratios for precipitation occurrence were generally less than 0.2 and showed a consistent reduction in the middle decade when station numbers were at their largest. The actual signals (not plotted) were reasonably constant after the 1960s suggesting that there are intrinsic limits to the overall complexity in occurrence structure that can be obtained from the station network with this technique. The small occurrence signals, which were generally less than 400, confirmed the appropriateness of fitting the occurrence surfaces with an approximate thin-plate smoothing spline using 1000 knots. Square root GCV values for the positive precipitation surfaces ranged between about 2 and 5 mm over all years and were slightly larger in summer. These values are small in absolute terms, although it should be noted that, because of the preponderance of small precipitation values, surface means generally ranged between 1 and 6 mm. Further examination of precipitation errors is undertaken below to examine errors in seasonal and annual totals and in precipitation extremes.
b. Model assessment—Withheld data for southern Canada
Withheld data errors for the daily temperature surfaces are summarized by season in Table 3. The values are generally consistent with the surface diagnostics in Table 1 and Fig. 4. Predicted values of maximum temperature tracked actual values closely. Mean absolute errors (MAE) varied between about 0.5° and 2.0°C, with an overall average of 1.1°C. Root-mean-square errors (RMSE) averaged around 1.8°C, in close agreement with the average square root GCV values of 1.6°C. The maximum temperature withheld errors showed the same seasonal variation exhibited by the root GCV values, with errors generally smallest in summer and autumn and largest in winter. The biases, though well below measurement error, are statistically significant at the 1% level except for summer. The largest bias occurred in winter.
Mean absolute errors were slightly higher for minimum temperature, varying between about 0.5° and 2.5°C, with an overall average of 1.6°C. RMSEs averaged around 2.2°C, in good agreement with the average square root GCV values of 2.0°C. Again, there was evidence of seasonal variation, with the smallest errors occurring in the summer and the largest in the winter. These biases were again well below measurement error but statistically significant at the 1% level.
Withheld errors in temperature extremes are summarized by midseason months in Table 4. The temperature extremes were reasonably well matched, although the mean errors differ significantly from zero at the 1% level. The upper 95th percentiles of daily maximum temperature were underestimated on average by 0.5°C in October (midautumn) up to 1.0°C in January (midwinter). The lower 95th percentiles of daily minimum temperature were overestimated on average by 0.6°C in July (midsummer) to 1.0°C in October (midautumn).
Withheld data errors for interpolated daily precipitation occurrence and amount are summarized by season in Table 5. Daily precipitation occurrence errors are also summarized over all years in Fig. 5. On average, the models correctly predicted the occurrence of daily precipitation at withheld stations with 83% accuracy. There was some seasonal variation, with predictive accuracy generally exceeding 83% in spring, summer, and autumn whereas predictive accuracy in winter was around 80%, reflecting observational difficulties in winter. Of the cases in which precipitation occurrence was predicted incorrectly, 7% were “false negatives” and 10% were “false positives.” This suggests that the models have little bias with respect to predicting precipitation occurrence. However, Fig. 5 also shows that the false positives increased and the false negatives decreased after 1977 when precipitation observations changed from imperial to metric. MAE in daily positive precipitation averaged around 2.9 mm and RMSE averaged around 5.2 mm. The average RMSE were somewhat larger than the average root GCV errors of 3.7 mm in Table 2, but, as in Table 2, the withheld errors in daily positive precipitation were largest in summer. There was a small negative bias in predicted daily positive precipitation, statistically significant at the 1% level across all months, with an overall average of around −0.6 mm.
We also compared actual and predicted seasonal and annual totals at the withheld stations by summing daily values over each season. These are summarized in Table 6. The seasonal MAE ranged between 5% and 30% with an overall average of around 14%. Percentage errors were distinctly larger in winter and smallest in summer and autumn. The models tended to slightly underestimate seasonal precipitation in spring and summer, and overestimate winter precipitation. The biases are statistically significant at the 1% level except in fall. The larger errors in winter are consistent with observational difficulties. The average percentage error in the predicted annual totals was only 9% with a small negative bias overall of −1.6%.
Withheld errors in daily, seasonal, and annual precipitation extremes are shown in Table 7. Both daily and seasonal extremes were underestimated, with all mean differences statistically different from zero at the 1% level except for winter. Extremes were the least underestimated in winter and most underestimated in summer, reflecting the spatial complexity of summer convective weather systems. The errors in predicting daily 95th upper percentiles were small in absolute terms, averaging −2.9 mm, but these are relatively large in percentage terms, averaging around −26%. However the average percentage error in predicting seasonal upper 95th percentiles was only −8%. The percentage error in estimating annual upper 95th percentiles was −7.3%.
c. Model assessment—Withheld data for northern Canada
While the main focus of the error analysis for this study is southern Canada, it is of interest to assess the accuracy of the interpolated surfaces over the data sparse region of northern Canada. Withheld validation errors can be calculated for each data point by extracting the diagonal values of the influence matrix for the fitted spline surfaces and applying the “leaving out one lemma” (Wahba 1990). This was done for selected days in midseason months for the surfaces for four selected years (1965, 1975, 1985, and 1995). Withheld errors in maximum temperature were around 100% larger for data points north of the line given by Eq. (1) than for data points south of the line. Withheld errors in minimum temperature were around 80% larger in the north. Withheld errors in precipitation occurrence were around 40% larger in the north, while withheld errors in daily positive precipitation were around 50% larger in the north. Thus, withheld errors for northern Canada were significantly larger, but not as large as might be expected, given that the data in the north are around a factor of 20 sparser than the data in the south. The surfaces can thus be used with moderate confidence in northern Canada, although daily extremes are likely to be poorly represented in the north.
The errors in the temperature and precipitation surfaces show consistent patterns reflecting contributing atmospheric processes, inadequacies in data networks, and observational difficulties in winter months. Precipitation when falling as snow is notoriously difficult to measure, with reported measurement errors of up to 50% (Sevruk 1982; Goodison 1978). The errors also reflect the minor impact of differing definitions of the climate. Overall, errors in predicted temperature values were about 0.5°C larger than those obtained when modeling monthly values with the same technique (McKenney et al. 2006). This is reasonable given the greater spatial complexity of daily weather patterns. As for most studies, errors in daily minimum temperature were consistently larger, by around 0.5°C, than errors in daily maximum temperature. Spatial patterns in maximum temperature are strongly influenced by ground elevation intersecting with linear atmospheric lapse rates but patterns in minimum temperature are influenced by additional processes including cold air drainage, which inverts local lapse rates, particularly in winter, and the latent heat of water condensation when the minimum temperature is close to the dewpoint temperature. Thus, Hutchinson (1991) found maximum temperature lapse rates for Tasmania of around 8°C (1000 m)−1, reasonably close to the dry adiabatic lapse rate of 9.8°C (1000 m)−1, whereas minimum temperature lapse rates varied between 3°C in winter and 5°C in summer. Bolstad et al. (1998) found similar behavior in maximum and minimum temperature lapse rates over the southern Appalachian Mountains. They also noted the impact of topographic aspect on maximum temperature.
The spatially varying dependence on elevation of the trivariate spline analyses fitted here permitted spatial variations in lapse rates to be accommodated. Plots of withheld daily maximum and minimum temperature errors versus elevation are shown in Figs. 6a,b. They indicate very little residual bias in the modeled temperature dependences on elevation. There is a slight positive bias in the maximum temperature dependence and a slight negative bias in the minimum temperature dependence. These biases are associated with data from the nonsummer months, consistent with the more complex controlling processes and observational difficulties in these months. The largest outliers from the fitted trend lines are associated with remote stations in the Yukon and British Columbia with few close neighbors and atypical elevation gradients due to pooling of cold air east of the Rocky Mountains.
The high spatial complexity of daily precipitation leads to significant errors in the daily precipitation surfaces, although these are ameliorated when data are summed to seasonal and annual totals. There are also moderate differences between daily root mean GCV values and daily RMSE values. These reflect the inadequacy of the data network to fully sample the spatial complexity of daily precipitation. They also reflect the more even sampling of elevation gradients by the withheld data than the main data network that is dominated by low-elevation stations. However, the overall elevation dependence of the interpolated precipitation appears to be sound, with Fig. 6c showing virtually no residual trend in bias on elevation. The largest residuals from the fitted trend line are associated with a remote west coastal station and a remote station on the northern shore of the Gulf of St. Lawrence. In fact, there is no residual trend on elevation in summer and autumn consistent with the smaller seasonal percentage precipitation errors in summer and autumn shown in Table 6. The smaller precipitation errors in autumn are consistent with the larger correlation length scales associated with the broadscale synoptic systems that occur in this season (Milewska and Hogg 2001). These broadscale systems also appear to have an impact on the generally smaller errors in maximum temperature in autumn, as shown in Tables 1, 2, and 4.
a. Comparison with other studies
Relatively few studies have generated historical daily climate grids of the spatial coverage and resolution presented here. We restrict attention here to studies that have provided representative error estimates based on a validation procedure that withheld data from the analysis. The key summary error statistics are presented in Table 8. The statistics include MAE of daily maximum and minimum temperature, MAE of daily precipitation on wet days, mean percentage error in predicting daily precipitation occurrence, and mean percentage error in annual precipitation. None of the comparison studies provided estimates of errors in predicting temperature or precipitation extremes.
The “DAYMET” procedure was developed and applied to the northwestern United States by Thornton et al. (1997). It shares some broad similarities with the present study in fitting local lapse rates on elevation and in adapting to the spatially varying density of the data. The method also fitted precipitation occurrence and wet day precipitation separately. As for the present study the interpolated values had minimal bias. The method had somewhat larger errors in temperature and in annual precipitation, despite the station density being 3 times higher than that of the present study. A likely contributor to these errors is the larger elevation coverage, about 3 times that of the present study.
This is confirmed by the application of DAYMET to Austria by Hasenauer et al. (2003). This study had 6 times the station density of the present study and less elevation coverage than the application to the northwestern United States. It obtained errors in temperature and daily precipitation occurrence that were slightly smaller than those of the present study. It was found to overestimate precipitation at elevations above 1500 m but, as for the present study, found no significant trend in precipitation errors on elevation for stations below 1500 m. The study by Xia et al. (2001) had limited spatial coverage over Bavaria but their preferred method, trivariate thin-plate splines, had similarly good results for temperature over a network that was half as dense as the network of Hasenauer et al. (2003).
The study by Hunter and Meentemeyer (2005) developed a climatologically aided method that applied ordinary kriging to daily differences from monthly values obtained by the Precipitation-Elevation Regressions on Independent Slopes Model (PRISM) method (Daly et al. 1994). The monthly PRISM values provided a spatially varying dependence on elevation for the interpolated daily values. The climatologically aided approach achieved a moderate reduction in error over using ordinary kriging alone, but the errors in maximum temperature were larger than those of the present study and the predicted maximum and minimum temperatures had a significant negative bias of around −0.4°C.
The study by Shen et al. (2001) suggests that inverse distance interpolation can achieve errors slightly larger than those obtained here when applied to denser data networks with limited elevation coverage. The GIDS method (Nalder and Wein 1998) was the method of choice by Stahl et al. (2006) in interpolating daily maximum and minimum temperature over British Columbia. It performed marginally better than the authors’ implementation of the DAYMET method. GIDS incorporates a spatially varying dependence on elevation by combining local multiple regression on position and elevation with inverse distance squared interpolation. It appears to achieve similar performance to that of the present study over a data network with similar density and elevation coverage. Its main deficiency, as noted by Stahl et al. (2006) and Price et al. (2000), is a tendency to produce occasional large outliers in data-sparse high-elevation areas where lapse rates can be poorly determined by local methods.
Two temperature interpolation studies with fixed lapse rates over extensive regions have yielded remarkably good results. Jarvis and Stuart (2001) applied a partial thin-plate smoothing spline with a fixed dependence on elevation to daily minimum and maximum temperature over the United Kingdom and Wales in 1976. In this analysis the lapse rate was determined automatically from data for each day. The RMS cross-validation errors were somewhat smaller than the other studies cited here, reflecting the limited elevation coverage of the Jarvis and Stuart study and perhaps broader-scale controlling processes and higher-quality data. Dodson and Marks (1997) applied inverse distance interpolation to elevation detrended temperatures with a fixed lapse rate of 3.9°C km−1 for both maximum and minimum temperature over the northwestern United States. Pooled MAE were around 1.3°C, quite similar to the result obtained by the present study. While it is plain in general that lapse rates should be spatially adaptive to spatially varying processes, these last two studies show that the elevation dependence should be quite stable to give good results.
Through a concerted modeling effort, daily grids of minimum and maximum temperature and precipitation have been produced for Canada for the period 1961–2003. Additional years will be added as data become available. The grids have already been used in preliminary studies to evaluate high-resolution precipitation products from satellites in high latitudes and in agricultural applications. The method used in the present study is well suited to a highly automated approach that is necessitated by the magnitude of daily data. The method appears to be more accurate than the other methods cited when data network density and elevation coverage are taken into account. However, as noted by Thornton et al. (1997), differences between methods are not particularly large. The key strength of the trivariate spline method appears to be its global dependence on all data. This permits robust and stable determination of spatially varying dependences on elevation. This is particularly important in the data-sparse areas that tend to prevail in high-elevation regions. Overall error estimates for southern Canada were quite low, with mean absolute errors of 1.1° and 1.6°C for daily maximum and minimum temperatures, respectively, and about 9% for annual precipitation. Daily temperature extremes were also well modeled because of the strong broadscale topographic control on daily temperature.
The results for daily precipitation compare well with other studies but the percentage errors in daily precipitation and daily precipitation extremes are relatively large. This is an inevitable consequence of the high spatial complexity associated with daily precipitation that is below the spatial resolution of most data networks. Further improvements in the current methodology are aimed at combining the interpolation of precipitation occurrence and positive amount so that each can inform the other and remove occasional sharp disjunctions between interpolated positive precipitation and interpolated dry areas. This may offer a modest improvement but the inadequate coverage of most data networks will remain. Methods that take account of long-term mean precipitation patterns, as suggested by Hutchinson (1995b), may also offer some improvement. The study by Tait et al. (2006) offers encouraging evidence in this regard. Nevertheless, daily precipitation surfaces interpolated from standard data networks should be used with care, particularly when daily precipitation extremes are of critical importance. Such applications may be better served probabilistically by interpolating the parameters of the statistical distributions describing the overall temporal and spatial patterns of extremes (Hutchinson 1995b). Fortunately, when the daily precipitation values of this study are summed to seasonal and annual values the errors are modest and these values can be used with some confidence across southern Canada.
The success of any effort involving historical data is also limited by deficiencies and inaccuracies in the data itself. Deficiencies in data quality and network coverage and inconsistencies in the definition of the climate day have contributed to the errors tabulated here. We put substantial effort into checking and correcting the Canadian historical daily data record, which provided a larger, higher-quality dataset for the current work. Efforts to improve the daily historical record such as addressing the climate day issue and correcting climate station locations are ongoing.
Daily climate data and station geographic coordinates were obtained from Environment Canada in Downsview. Thanks are given to Climate Services for supplying the data and others at Environment Canada for supporting our efforts in this area. Funding and support for the work was provided by Agriculture and Agri-Food Canada, Natural Resources Canada, Canadian Forest Service, Environment Canada, and the Australian National University. The authors thank the anonymous reviewers for comments that improved the manuscript but accept responsibility for any remaining errors.
Corresponding author address: Michael Hutchinson, Fenner School of Environment and Society, Australian National University, Canberra ACT 0200, Australia. Email: firstname.lastname@example.org