Thin-plate Smoothing Spline Research Articles

Predicting lodgepole pine site index from climatic parameters in Alberta

We sought to evaluate the impact of climatic variables on site productivity of lodgepole pine (Pinus contorta var. latifolia Engelm.) for the province of Alberta. Climatic data were obtained from the Alberta Climate Model, which is based on 30- year normals from the provincial weather station network. Mapping methods were based on ANUSPLIN, Hutchinson's thin-plate smoothing spline in four dimensions (latitude, longitude, elevation, climatic variable). Site indices based on stem analysis (observed dominant height at an index age of 50 years at breast height) were used as a measure of forest site productivity. A total of 1145 site index plots were available for lodgepole pine, the major forest species in Alberta. Regression analyses were used to predict site index as a function of climatic variables for each plot. The strongest linear predictors of site index were growing degree days > 5° (GDD5), the Julian date when GDD5 reaches 100 (D100), and July mean temperature (MTWM).A nonlinear model with D100 as the predictor variable was chosen as the final model. Both the observed and the predicted site indices from the 1145 locations were interpolated using ANUSPLIN and mapped using ArcView. We concluded that climate is an important component of site productivity, accounting for about one quarter of the variation in lodgepole pine site index across the province. Key words: site index, climate variable, Alberta Climate Model, ANUSPLIN, site productivity map, lodgepole pine

Biomass and biomass change in lodgepole pine stands in Alberta

We describe methods and results for broad-scale estimation and mapping of forest biomass for the Canadian province of Alberta. Differences over successive decades provided an estimate of biomass change. Over 1500 permanent sample plots (PSP) were analyzed from across the range of lodgepole pine (Pinus contorta var. latifolia Engelm.), the major forest tree species of Alberta. The PSP network is densest in stands aged between 70 and 100 years and is well-represented by stands of all ages to 150 years of age. Stand biomass (Mg ha(-1)) was estimated for each PSP plot as the sum of the respective biomass components for each tree (live and standing dead). The biomass components for live trees were stem, bark, branches, foliage and roots. The components for standing dead trees excluded foliage. Equations from previous biomass studies were used for biomass component estimation. Biomass estimates of additional non-tree components were attempted, but without much success. Biomass of the soil organic layer was estimated once on 452 PSPs and a mean estimate of total dead fuels on the ground (28.4 Mg ha(-1)) was available only for the entire distribution of lodgepole pine. However, values of these two components were essentially constant over time and therefore did not alter the analysis or conclusions obtained by analyzing total tree biomass alone. We then used this spatial network of 1549 plots as the basis for mapping biomass across Alberta. Mapping methods were based on Australian National University SPLINe (ANUSPLIN) software, Hutchinson's thin-plate smoothing spline in four dimensions (latitude, longitude, elevation and biomass). Total tree biomass (mean = 172 Mg ha(-1)) was dominated by stem biomass (mean = 106 Mg ha(-1)), which was an order of magnitude greater than the mean estimates for the bark (11 Mg ha(-1)), branch (12 Mg ha(-1)) and foliage (12 Mg ha(-1)) components. A close relationship was found between total tree biomass and stand stem volume (R(2) = 0.992 with n = 3585; note that volume and biomass were calculated independently). We compared total tree biomass for two decades, the 1980s and the 1990s. After correcting for changes in harvest removals over time, the mean change in total biomass was positive (0.99 Mg ha(-1) year(-1)) and differed significantly from zero (n = 421; P < 0.001). Estimates ranged from -13.9 to 8.0 Mg ha(-1) year(-1). The heart of the lodgepole pine distribution (primarily the Foothills subregions) showed an increase in biomass, whereas isolated pockets of lodgepole pine in the boreal northern subregion indicated a decline in biomass.

Thin plate smoothing spline interpolation of daily rainfall for New Zealand using a climatological rainfall surface

AbstractThis study presents a method for estimating daily rainfall on a 0.05° latitude/longitude grid covering all of New Zealand for the period 1960–2004 using a second order derivative trivariate thin plate smoothing spline spatial interpolation model. Use of a hand‐drawn (and subsequently digitised) mean annual rainfall surface as an independent variable in the interpolation is shown to reduce the interpolation error compared with using an elevation surface. This result is confirmed when long‐term average annual rainfall data, derived from the daily interpolations, are validated using long‐term river flow data. Copyright © 2006 Royal Meteorological Society

Automatic computation of hierarchical biquadratic smoothing splines with minimum GCV

A computationally efficient numerical strategy for fitting approximate minimum GCV bivariate thin plate smoothing splines to large noisy data sets was developed. The procedure discretises the bivariate thin plate smoothing spline equations using biquadratic B-splines and uses a nested grid SOR iterative strategy to solve the discretised system. For efficient optimisation, the process incorporates a double iteration that simultaneously updates both the discretised solution and the estimate of the minimum GCV smoothing parameter. The GCV was estimated using a minimum variance stochastic estimator of the trace of the influence matrix associated with the fitted spline surface. A Taylor series expansion was used to estimate the smoothing parameter that minimises the GCV estimate. The computational cost of the procedure is optimal in the sense that it is proportional to the number of grid points supporting the fitted biquadratic spline. Convergence was improved by adding a first order correction to the solution estimate after each smoothing parameter update. The algorithm was tested on several simulated data sets with varying spatial complexity and noise level. An accurate approximation to the analytic minimum GCV thin plate smoothing spline was obtained in all cases.

On the Horizontal Scale of Elevation Dependence of Australian Monthly Precipitation

Abstract Determination of the scale of the interaction between precipitation and topography is important for the accurate interpolation of rainfall in mountainous areas and also provides insight into the physical processes involved. In this paper, trivariate thin-plate smoothing splines are used to investigate the scale of interaction between monthly precipitation and topography by interpolating monthly rainfall over three subregions of the Australian continent, incorporating different climatic conditions and rainfall types. The interpolations are based upon elevations derived from digital elevation models (DEMs) of various resolutions. All of the DEMs are local averages of version 2.0 of the 9-s-resolution DEM of Australia. The results suggest that the optimal scale of the interaction between precipitation and topography, as it pertains to the elevation-dependent interpolation of monthly precipitation in Australia, is between 5 and 10 km. This is in agreement with results of similar studies that addressed daily precipitation over Switzerland.

Spatial interpolation of large climate data sets using bivariate thin plate smoothing splines

Thin plate smoothing splines are widely used to spatially interpolate surface climate, however, their application to large data sets is limited by computational efficiency. Standard analytic calculation of thin plate smoothing splines requires O( n 3) operations, where n is the number of data points, making routine computation infeasible for data sets with more than around 2000 data points. An O( N) iterative procedure for calculating finite element approximations to bivariate minimum generalised cross validation (GCV) thin plate smoothing splines operations was developed, where N is the number of grid points. The key contribution of the method lies in the incorporation of an automatic procedure for optimising smoothness to minimise GCV. The minimum GCV criterion is commonly used to optimise thin plate smoothing spline fits to climate data. The method discretises the bivariate thin plate smoothing spline equations using hierarchical biquadratic B-splines, and uses a nested grid multigrid procedure to solve the system. To optimise smoothness, a double iteration is incorporated, whereby the estimate of the spline solution and the estimate of the optimal smoothing parameter are updated simultaneously. When the method was tested on temperature data from the African and Australian continents, accurate approximations to analytic solutions were obtained.

Spatial interpolation of monthly mean climate data for China

AbstractSpline interpolation techniques are used to develop a gridded climate database for China at a resolution of 0.01° in latitude and longitude. A digital elevation model (DEM) was developed at the same resolution to improve the accuracy of interpolation based upon the general spatial dependence of climate on topography. Climate data for the period 1971–2000 from meteorological stations in China were used to develop thin‐plate smoothing spline surfaces for monthly mean temperature and precipitation. A regularly gridded climate database was produced by coupling the spline surfaces with the underlying DEM. The summary statistics show interpolation errors for monthly temperatures varying within 0.42–0.83 °C and 8–13% for monthly precipitation. These estimates are superior to results produced by methods commonly used in China. The fine‐resolution spatial climate database has many potential applications in natural resource management. For example, it can be used as a baseline for climate change studies, in which potential distributions of flora and fauna can be predicted under the impact of climate change and priority areas for biodiversity conservation can be identified. Copyright © 2005 Royal Meteorological Society.

Very high resolution interpolated climate surfaces for global land areas

We developed interpolated climate surfaces for global land areas (excluding Antarctica) at a spatial resolution of 30 arc s (often referred to as 1-km spatial resolution). The climate elements considered were monthly precipitation and mean, minimum, and maximum temperature. Input data were gathered from a variety of sources and, where possible, were restricted to records from the 1950–2000 period. We used the thin-plate smoothing spline algorithm implemented in the ANUSPLIN package for interpolation, using latitude, longitude, and elevation as independent variables. We quantified uncertainty arising from the input data and the interpolation by mapping weather station density, elevation bias in the weather stations, and elevation variation within grid cells and through data partitioning and cross validation. Elevation bias tended to be negative (stations lower than expected) at high latitudes but positive in the tropics. Uncertainty is highest in mountainous and in poorly sampled areas. Data partitioning showed high uncertainty of the surfaces on isolated islands, e.g. in the Pacific. Aggregating the elevation and climate data to 10 arc min resolution showed an enormous variation within grid cells, illustrating the value of high-resolution surfaces. A comparison with an existing data set at 10 arc min resolution showed overall agreement, but with significant variation in some regions. A comparison with two high-resolution data sets for the United States also identified areas with large local differences, particularly in mountainous areas. Compared to previous global climatologies, ours has the following advantages: the data are at a higher spatial resolution (400 times greater or more); more weather station records were used; improved elevation data were used; and more information about spatial patterns of uncertainty in the data is available. Owing to the overall low density of available climate stations, our surfaces do not capture of all variation that may occur at a resolution of 1 km, particularly of precipitation in mountainous areas. In future work, such variation might be captured through knowledge-based methods and inclusion of additional co-variates, particularly layers obtained through remote sensing. Copyright © 2005 Royal Meteorological Society.

A COMPARISON OF SPATIAL PREDICTION METHODS USING INTENSE SPATIALLY-ACQUIRED WATER QUALITY DATA

Water quality information obtained through intensive spatial sampling using automated devices provides opportunities to monitor and forecast the spatial distribution of nutrients and phytoplankton concentrations, and help establish water circulation patterns in estuarine and coastal waters. To be cost effective, efficient sampling designs and estimation methodologies must first be developed. As a starting basis, we applied an original transect sampling design that was used to estimate the spatial distribution of chlorophyll a, salinity, and temperature in the Cienaga Grande de Santa Marta, a coastal lagoon in Colombia. We superimposed the transects over satellite images of the lagoon obtained in the period 1993-2001 to evaluate the efficiency and accuracy of using such transects to estimate the distribution of water quality variables. The satellite images were taken in 1993 (SPOT-3), 1995 (Landsat-6), and 1999 (Landsat-6), and water reflectance values were used as a “proxy” for the water quality variables. Spatial prediction using kriging and thin-plate smoothing splines were used to predict reflectance for a grid network of points taken from the images, and predictions were compared with observed values to compare methods and transect routes. Rapid changes in reflectance in short distances (for example , caused by phytoplankton blooms), complicated the analysis, and neither method proved superior over all transect routes and images, although the kriging predictor remained relatively consistent in performance over the various selected sampling routes.

An iterative procedure for calculating minimum generalised cross validation smoothing splines

This study analyses a simple iterative procedure for estimating minimum generalised cross-validation (GCV) univariate smoothing splines. The results provide guidelines for the development of a similar methodology to estimate minimum GCV bivariate thin plate smoothing splines. The methodology is based on the techniques described in Hutchinson [ ANZIAM J , 42(E):C774--C796, 2000], which uses nested grid SOR iterative methods in order to solve finite element thin plate smoothing spline systems efficiently for large data sets. The method also uses the stochastic approximation to the GCV developed by Hutchinson [ Commun. Stats --- Sim. and Comp. , 18:1059--1076, 1989]. A double iteration is used to produce increasingly accurate estimates of the minimum GCV smoothing parameter and the smoothing spline. First and second derivatives of the GCV with respect to the smoothing parameter are used to update the smoothing parameter. Convergence of the SOR iteration is improved significantly by correcting the solution estimate after each smoothing parameter update using the estimate of the derivative of the solution with respect to the smoothing parameter.

Using spatial interpolation to construct a comprehensive archive of Australian climate data

A comprehensive archive of Australian rainfall and climate data has been constructed from ground-based observational data. Continuous, daily time step records have been constructed using spatial interpolation algorithms to estimate missing data. Datasets have been constructed for daily rainfall, maximum and minimum temperatures, evaporation, solar radiation and vapour pressure. Datasets are available for approximately 4600 locations across Australia, commencing in 1890 for rainfall and 1957 for climate variables. The datasets can be accessed on the Internet at http://www.dnr.qld.gov.au/silo. Interpolated surfaces have been computed on a regular 0.05° grid extending from latitude 10°S to 44°S and longitude 112°E to 154°E. A thin plate smoothing spline was used to interpolate daily climate variables, and ordinary kriging was used to interpolate daily and monthly rainfall. Independent cross validation has been used to analyse the temporal and spatial error of the interpolated data. An Internet based facility has been developed which allows database clients to interrogate the gridded surfaces at any desired location.

Computational Aspects of Approximation to Scattered Data by Using ‘Shifted’ Thin-Plate Splines

A new multivariate approximation scheme to scattered data on arbitrary bounded domains in R d is developed. The approximant is selected from a space spanned (essentially) by corresponding translates of the ‘shifted’ thin-plate spline (‘essentially,’ since the space is augmented by certain functions in order to eliminate boundary effects). This scheme applies to noisy data as well as to noiseless data, but its main advantage seems to be in the former case. We suggest an algorithm for the new approximation scheme with a detailed description (in a MATLAB-like program). Some numerical examples are presented along with comparisons with thin-plate spline interpolation and Wahba's thin-plate smoothing spline approximation.

Geographic patterns and dynamics of Alaskan climate interpolated from a sparse station record.

Data from a sparse network of climate stations in Alaska were interpolated to provide 1-km resolution maps of mean monthly temperature and precipitation--variables that are required at high spatial resolution for input into regional models of ecological processes and resource management. The interpolation model is based on thin-plate smoothing splines, which uses the spatial data along with a digital elevation model to incorporate local topography. The model provides maps that are consistent with regional climatology and with patterns recognized by experienced weather forecasters. The broad patterns of Alaskan climate are well represented and include latitudinal and altitudinal trends in temperature and precipitation and gradients in continentality. Variations within these broad patterns reflect both the weakening and reduction in frequency of low-pressure centres in their eastward movement across southern Alaska during the summer, and the shift of the storm tracks into central and northern Alaska in late summer. Not surprisingly, apparent artifacts of the interpolated climate occur primarily in regions with few or no stations. The interpolation model did not accurately represent low-level winter temperature inversions that occur within large valleys and basins. Along with well-recognized climate patterns, the model captures local topographic effects that would not be depicted using standard interpolation techniques. This suggests that similar procedures could be used to generate high- resolution maps for other high-latitude regions with a sparse density of data.

Interpolation of Daily Global Solar Radiation with Thin Plate Smoothing Splines

Thin plate smoothing splines incorporating topographic dependence were used to interpolate daily global solar radiation in the Bavarian forest ecosystem monitoring network, with the degree of data smoothing determined by minimizing the generalized cross validation. A simple cross validation method was used to discuss the spatial distribution of mean relative errors at 18 forest climate stations. The results show that, from this network 14%–30% mean relative errors can be expected for most of these forest climate stations in summer, and 20%–30% mean relative errors can be found at a few of forest climate stations in winter. Time-averaging can reduce these interpolation errors. In this network a mean relative error of 10% can be expected for weekly and biweekly mean solar radiation at most of forest climate stations in summer. Large errors are related to low radiation amount under heavy cloud cover. Mean relative errors increase as daily global solar radiation decreases.

A comparison of two statistical methods for spatial interpolation of Canadian monthly mean climate data

Two methods for elevation-dependent spatial interpolation of climatic data from sparse weather station networks were compared. Thirty-year monthly mean minimum and maximum temperature and precipitation data from regions in western and eastern Canada were interpolated using thin-plate smoothing splines (ANUSPLIN) and a statistical method termed ‘Gradient plus Inverse-Distance-Squared’ (GIDS). Data were withheld from approximately 50 stations in each region and both methods were then used to predict the monthly mean values for each climatic variable at those locations. The comparison revealed lower root mean square error (RMSE) for ANUSPLIN in 70 out of 72 months (three variables for 12 months for both regions). Higher RMSE for GIDS was caused by more frequent occurrence of extreme errors. This result had important implications for surfaces generated using the two methods. Both interpolators performed best in the eastern (Ontario/Québec) region where topographic and climatic gradients are smoother, whereas predicting precipitation in the west (British Columbia/Alberta) was most difficult. In the latter case, ANUSPLIN clearly produced better results for most months. GIDS has certain advantages in being easy to implement and understand, hence providing a useful baseline to compare with more sophisticated methods. The significance of the errors for any method should be considered in light of the planned applications (e.g., in extensive, uniform terrain with low relief, differences may not be important).

Mapping rainfall fields and their ENSO variation in data-sparse tropical south-west Pacific Ocean region

Rainfall fields for the data-sparse tropical south-west Pacific Ocean region have been mapped by partial thin-plate smoothing spline surface modelling applied to island rainfall measurements, enhanced by the use of satellite observations of outgoing longwave radiation (OLR) as a regression covariate. The aim is to obtain spatially realistic rainfall maps, especially in the data-sparse areas between island groups, through a fully objective and statistically valid method that includes error estimates. The method has been applied to the region 4°N–24°S, 168°E–154°W. The rainfall data set initially comprised 57 stations, most with 40 year records. As a first step, a regression of annual OLR and rainfall for atolls only was formed and used to eliminate ‘outlier’ rainfall stations, all of which are on mountainous islands and thus are probably influenced orographically. The maps clearly show the spatial patterns and seasonal behaviour of the region's key meteorological features, namely, the South Pacific Convergence Zone (SPCZ), the southern edge of the Inter-tropical Convergence Zone (ITCZ), and the wedge shaped region of divergent easterlies lying between them. To identify ENSO variations, maps of 3-month seasonal rainfall were constructed from composites of eight El Nino (negative SOI) episodes and nine La Nina (positive SOI) episodes. These maps are relatively rough in appearance, but nevertheless they show the evolution of the spatial patterns through each composite episode and the strong and symmetrically opposite differences between them. Marked variations in the strength and position of the SPCZ are evident and the isohyets in the equatorial dry zone exhibit east–west shifts of nearly 3000 km relative to the average field. The rainfall variation at a particular location may be understood in terms of competition of influence among the changing features of the pattern, rather than as a simple linear function of the SOI. © 1998 Royal Meteorological Society

Site regions revisited: a climatic analysis of Hills' site regions for the province of Ontario using a parametric method

not available

SPATIAL MODELLING OF NEW ZEALAND TEMPERATURE NORMALS

Spatial modelling and prediction of normals of monthly temperatures in relation to geographical factors is examined for a data set of 100 New Zealand climate stations, most of which have 30-year records. Results are given for maximum and minimum temperature and for mean temperature and diurnal range. Latitude and altitude are well established as predictors of temperature. The application of a partial thin-plate smoothing spline model shows that there is no advantage in using this model to represent the temperature response as a surface in the latitude and longitude field. using additive models to examine non-linear responses, we find that the normals are linear with altitude and mostly exponential with distance-to-sea. Based on these findings, multiple regression analysis and the Akaike Information Criterion are used to find optimal prediction models. New predictors are identified as: distance-to-sea in non-linear exponential form to capture the oceanic influences in coastal regions; a regional index that divides the country into two parts based on exposure to, or sheltering from, the predominant west to south-westerly moist air flow; and a forest index that models the lower temperatures experienced at forest stations. The resulting regression models for the monthly mean temperature normals explain 86–95 per cent of the variance and give root mean square residuals within the range of 0ċ 43°C to 0ċ 62°C. We find that the use of the annual mean temperature normal as a predictor results in substantially improved models, which implies that there are non-random factors affecting site temperature regimes that remain unmodelled by the existing predictors. The magnitude, sign and seasonal variation of the coefficients of the regression models provide insights into spatial temperature variability and its associated climate processes.

Site regions revisited: a climatic analysis of Hills' site regions for the province of Ontario using a parametric method

A climatic analysis was conducted` of Hills' site regions for Ontario using spatially distributed models of long-term, mean monthly climate data. Gridded estimates of eight climatic variables were generated by coupling thin-plate smoothing spline surfaces (fitted as a trivariate function of longitude, latitude, and elevation) to a 1-km resolution digital elevation model of Ontario. The gridded climate estimates were used to characterize the climate of each site region. New climatic classifications for the province were calculated from these data using a clustering algorithm. They were also compared with the existing site region boundaries. The analyses confirm many of the general climatic gradients established by Hills, but also suggest where major revisions should be made. The climatic models present new analytical opportunities for mapping plant–climate response in the context of forest planning and management.

The implementation of a Generalized Cross Validation algorithm using deflation techniques for linear systems

The fitting of a thin plate smoothing spline to noisy data using the method of minimizing the Generalized Cross Validation (GCV) function is computationally intensive involving the repeated solution of sets of linear systems of equations as part of a minimization routine. In the case of a data set of more than a few hundred points, implementation on a workstation can become unrealistic and it is then desirable to exploit high performance computing. The usual implementation of the GCV algorithm performs Householder reductions to tridiagonalize the influence matrix and then solves a sequence of tridiagonal linear systems which are updated only by a scalar value (the minimization parameter) on the diagonal. However, this approach is not readily parallelizable. In this paper the deflation techniques described by Burrage et al. (1994), which are used to accelerate the convergence of iterative schemes applied to linear systems, will be adapted to the problem of minimizing the GCV function. This approach will allow vector and parallel architectures to be exploited in an efficient manner.