Spatial interpolation techniques: their applications in regionalizing climate-change series and associated accuracy evaluation in Northeast China

ABSTRACT An accurate gridded climatological temperature data-set can be a reliable basis for studying the issues that concern climate change, natural disasters, among others. In this study, the climate standard value data and annual observed data collected by 104 meteorological sites in Northeast China from 1971 to 2000 are used to interpolate the annual mean air temperature of the region with different spatial interpolation techniques. Efforts have also been made to verify the accuracy of the results generated by each interpolation method with error indicator, temperature characteristic value and temperature variation curves. The results show at the spatial scale, the partial thin-plate smoothing spline (PTPSS) scheme generates the most desirable temperature interpolations with a root-mean-square error at 0.34 °C and an average standard error at 0.52 °C. At the temporal scale, the PTPSS-based 1971–2000 temperature curves agree the best with the gridded temperature curves, with the correlation coefficients of means, minimums and maximums between the two being larger than 0.9. Meanwhile, the PTPSS produce the smallest error approaching zero among all the interpolation schemes tested in the study. In this context, the PTPSS scheme stands out as the most desirable interpolation method for the annual mean temperature in Northeast China.


Introduction
In the past 50 years or so, China has witnessed significant changes of climate under the influence of global climate change. In a National Assessment Report on Climate Change released by China in 2007, it is stated that climate change will cast a major impact on China's natural ecosystems and socioeconomic systems, particularly on agricultural, animal raising and water supply activities. A range of climate predictions made so far have indicated that China will see a further rise of surface temperatures in the coming 50-100 years (Wang 2001;China's National Climate Change Programme 2007), and a noticeably enhanced frequency and intensity of extreme weather/climate events. The increasingly deepened climate change study calls for the extraction of the spatial and temporal distribution of climatic elements of a geographical region from the observed meteorological data collected by the sparsely and unevenly distributed meteorological observed stations. As a major atmospheric variable, air temperature may directly affect the physiological and ecological processes on the Earth, including evapotranspiration and plant growth (Stahl et al. 2006;Li et al. 2013). In this context, building an accurate and comprehensive regional temperature data-set would create a reliable basis for modelling the regional ecosystems, and is hence meaningful for studies in a variety of areas, for example, ecology, hydrology, climate change, among others.
A lot of temperature data-sets have been produced, such as the 2.0 £ 2.0 global daily temperature data-set (GISS (Goddard Institute for Space Studies) Surface Temperature Analysis, GISTEMP) released by the U.S. National Oceanic and Atmospheric Administration (NOAA), and the 1.0 £ 1.0 gridded daily, monthly and annual temperature data-set produced by China Meteorological Administration (CMA). In addition, some research teams developed the precipitation data-sets of different resolutions to accommodate their own research needs (Shen et al. 2005;Sinha et al. 2006). Though difference in spatial and temporal scales, those data-sets are unanimously spatially interpolated. Spatial interpolation is a mathematical process designed to predict the value at a point or in a region where the value is unknown based on the known discrete data or partition data from a set of sample points that are calculated on their mathematical relationships (Qian et al. 2010). Spatial interpolation can either be a deterministic or a stochastic (also known as geostatistical) process (Piazza et al. 2011), depending on the mathematical principles used. The kriging technique is mostly used in geostatistical interpolation (Figure 1). The study of interpolation techniques mainly works on the following two aspects: (1) interpolate meteorological elements in a region using different interpolation techniques, in an attempt to compare the accuracy of the techniques used and analyse the spatial and temporal variations of meteorological elements in the region (Zhuang & Wang 2003;Cai et al. 2006;Piazza et al. 2011;Bostan et al. 2012;Wagner et al. 2012;Teegavarapu et al. 2012); and (2) raised interpolation accuracy by improving the exiting interpolated methods but the improved interpolation methods are usually only fit for the specific study area not for other regions (Liu et al. 2004;Pan et al. 2004;Shi et al. 2011;Tang et al. 2012;Kilibarda et al. 2014). As a matter of fact, the variation of meteorological elements, especially at a large scale up to a region, could be very sophisticated. For example, China has a vast territory that allows large geographic and climatic differences. To improve the spatial accuracy of climate data, one has to make climate zone based on geographical divisions and climatic characteristics of each division before comparing the accuracy of different interpolation techniques (Figure 1) for the most desirable results. The results, as such, can be used to justify the application of a given interpolation technique in developing a large-scale meteorological data-set. On the other hand, one has to recognize that different meteorological elements may have a different climate division standard. In this context, the interpolation shall be made in line with the zoning of specific meteorological elements. For example, the north-eastern part of China is a typical climate zone that embraces three provinces (Heilongjiang, Jilin and Liaoning) and four cities of Inner Mongolia. The region is noticeably differed in air temperature, precipitation and other meteorological elements compared with other regions. Meanwhile, the region is a major commodity grain producer in the country. Meteorological elements are the major factors that affect the crop growth. Temperature variations would directly affect the length of crop growth period, and hence the crop yield. Therefore, an accurate temperature data-set of long time series is important for studying climate change, crop growth/development and natural disasters in the region. In recent years, a range of studies have been launched to regionalize the meteorological elements of north-eastern part of China (Zhuang & Wang 2003;Cai et al. 2006;Fu et al. 2009;He et al. 2013), though most of them employed only one interpolation technique and took the interpolation results as the input data for basic research works without the verification of accuracy, especially at the temporal scale. A few others worked to evaluate interpolation accuracy based on the comparisons of different interpolation techniques (Piazza et al. 2011;Bostan et al. 2012;Teegavarapu et al. 2012). Unfortunately, those studies have not sorted out a desired interpolation technique that can be applied in an accurate and universal manner.
The current study works to interpolate the annual air temperatures of north-eastern part of China, using different interpolation techniques (Figure 1) in an attempt to sort out the most desirable one that can mirror the spatial and temporal variations of the annual air temperatures, and to build an accurate and reliable annual air temperature data-set of long time series that can be used by climate-change study, ecosystem modelling, among others.

Data
The meteorological data used in the study are obtained from CMA (http://cdc.cma.gov.cn/home. do), including (1) Chinese surface climate standard value data-set between 1971 and 2000. The climate standard value is defined as the average value of 30-year climatic elements that can be used to describe the climate as well as reference the object current weather conditions. The World Meteorological Organization recommends that each country take 30 years as the standard and sets the initial stage of the time period for 1901-1930. The next period of time is the interval of 30 years (i.e. 1931-1960, 1961-1990) and the climatic standard value is updated every 10 years. Climatic standard value gets more and more attention which can reflect the climate anomalies (http://www.nmic.gov.cn/ news/news-2012-07-11-313419677295479.htm). The data-set contains many climatic elements, such as temperature, precipitation, pressure, wind speed and others. (2) Chinese surface climate data-set year by year contains the mean climatic value of 12 months in every year (from 1971 to 2000). (3) The 1 £ 1 gridded Chinese surface temperature data-set year by year which produced by the improved Kriging method and 1 £ 1 digital elevation model data based on GTOPO 30 data (http://data.cma.cn/). Both the data-sets (2) and (3) contain the similar climatic elements in the Chinese surface climate standard value data-set, and we only used mean temperature at 104 meteorological stations ( Figure 2) in the north-eastern part of China. The stations are fairly uniformly distributed across the region, with a somewhat higher spatial density in the southern part than the northern part. The 1-km elevation data of north-eastern China are obtained from National Aeronautics and Space Administration (NASA, http://srtm.csi.cgiar.org).

Interpolation techniques
The interpolation methods used in this study ( Figure 1) are not new and have been frequently described in the literatures and monographs (Yao et al. 2002;Zhuang & Wang 2003;Denby et al. 2005;Cai et al. 2006;Li & Heap 2008;Knnotters et al. 2010). Therefore, we restrict ourselves to a brief description.

Deterministic interpolation.
(1) Global polynomial (GP) GP interpolation fits a smooth surface that is defined by a mathematical function (a polynomial) to the input sample points. The GP surface changes gradually and captures coarse-scale pattern in the data. Conceptually, GP interpolation is like taking a piece of paper and fitting it between the raised points (raised to the height of value). The result from GP interpolation is a smooth surface that represents gradual trends in the surface over the area of interest. For detailed information, see ArcGIS Resources (http://resources.arcgis.com/en/home/).
(2) Local polynomial (LP) While GP interpolation fits a polynomial to the entire surface, LP interpolation fits many polynomials, each within specified overlapping neighbourhoods. The search neighbourhood can be defined by using the size and shape, number of neighbours, and sector configuration. LP interpolation, on the other hand, fits the specified-order (zero, first, second, third and so on) polynomial using points only within the defined neighbourhood. The neighbourhoods overlap, and the value used for each prediction is the value of the fitted polynomial at the centre of the neighbourhood. For detailed information, see ArcGIS Resources (http://resources.arcgis.com/en/home/).
(3) Inverse distance weighting (IDW) IDW method estimates the values of an attribute at unsampled points using a linear combination of values at sampled points weighted by an inverse function of the distance from the point of interest to the sampled points. IDW is based on the assumption that the similarity of values of a spatially distributed numerical variable z decreases with distance d or a power p of the distance. This decline of similarity is incorporated in the interpolation aŝ (1) where s 0 indicates the location to be interpolated to, s i (i D 1, …, n) indicates the n locations, where z has been observed, and d i is the distance between x 0 and x i . The main factor affecting the accuracy of IDW is the value of power parameter. Weights diminish as the distance increases, especially when the value of the power parameter increases, so nearby samples have a heavier weight and have more influence on the estimation, and the resultant spatial interpolation is local. For detailed information, see e.g. Isaaks and Srivastava (1989) and Knnotters et al. (2010).
(4) Spline functions (SF) SF are polynomials which are fitted in a flexible way through observations of z which are commonly applied to perform the smooth multivariate interpolation of the data that are irregularly distributed (Liu et al. 2008). Its smoothing parameters can be borrowed to achieve the optimal balance between the fidelity and the surface smoothness, ensuring the smoothness and continuity of the interpolated surface and associated accuracy and reliability (Yan 2004). Partial thin-plate smoothing splines (PTPSS) can in fact be viewed as a generalization of standard multivariate linear regression, in which the parametric model is replaced by a suitably smooth non-parametric function. The degree of smoothness, or inversely the degree of complexity, of the fitted function is usually determined automatically from the data by minimizing a measure of predictive error of the fitted surface given by the generalized cross-validation (GCV) (Li & Heap 2008). For detailed information, see Denby et al. (2005) and Hutchinson (2004).
The output statistics are best interpreted in relation to the partial spline model for N observed data values z i given by where each x i is a d-dimensional vector of spline-independent variables, f is an unknown smooth function of the x i , each y i is a p-dimensional vector of independent covariates, b is an unknown pdimensional vector of coefficients of the y i and each e i is an independent, zero-mean-error term with variance w i s 2 , wherew i is termed the relative error variance (known) and s 2 is the error variance which is constant across all data points, but normally unknown (Hutchinson 2004). The model reduces, on the one hand, to an ordinary thin-plate spline model when there are no covariates (p D 0) and to a simple multivariate linear regression model, on the other hand, when f ðx i Þ is absent. The latter possibility is not currently permitted by ANUSPLIN.

Geostatistical interpolation.
Kriging refers to a group of geostatistical interpolation methods in which the value at an unobserved location is predicted by a linear combination of the values at surrounding locations, using weights according to a model that describes the spatial correlation (Knotters et al. 2010). The statistical basis of kriging enables to quantify the accuracy of the predicted values by means of the kriging variance. The kriging variance is a measure of the accuracy of the interpolated values, or, in other words, a measure of the uncertainty about the true values. Statistically speaking, it starts from variable correlation and variability, allowing an unbiased optimal regionalization of the variables from a limited area. It asks for the spatial correlation between the regionalized variables and makes projections for certain points in the study area, based on the semivariogram and the spatial distribution of neighbouring areas. Kriging, like most interpolation techniques, is built on the assumption that things that are close to one another are more alike than those farther away (quantified here as spatial autocorrelation). The semivariogram is a means to explore this relationship. Pairs that are close in distance should have a smaller measurement difference than those farther away from one another. The extent that this assumption is true can be examined in the empirical semivariogram (Johnston et al. 2001).
Semivariance (g) of Z between two data points is an important concept in geostatistics and is defined as where h is the distance between points x i and x 0 and gðhÞ is the semivariogram (commonly referred to as variogram) (Webster & Oliver 2001). Variogram modelling and estimation is extremely important for structural analysis and spatial interpolation (Burrough & McDonnell 1998). The formula involves calculating half the difference squared between the values of the paired locations.
Once you have created the empirical semivariogram, you can fit a line to the points forming the empirical semivariogram model. The modelling of a semivariogram is similar to fitting a leastsquares line in regression analysis. The empirical semivariogram is fitted by an analytical function (model), the geostatistical analyst provides the following functions to choose from the model empirical semivariogram: circular, spherical, tetraspherical, pentaspherical, exponential, Gaussian, rational quadratic, hole effect, K-Bessel, J-Bessel and stable. The selected model influences the prediction of the unknown values, particularly when the shape of the curve near the origin differs significantly, a spherical type that rises at first and then levels off for larger distances beyond a certain range. The steeper the curve near the origin, the more influence the closest neighbours will have on the prediction. As a result, the output surface will be less smooth. Each model is designed to fit different types of phenomena more accurately (Burrough & McDonnell 1998;Johnston et al. 2001;Webster & Oliver 2001;Pebesma 2004;Denby et al. 2005). The formula for the spherical model is given as follows: where u s is the sill value, H is the lag vector and h is the length of H(distance between two locations), and u r is the range of the model. The basic goal is to calculate the parameters (i.e. range, sill and nugget) of the curve to minimize the deviations from the points according to some criterion. When the fitted curve does not pass through the origin, it denotes the existence of the so-called nugget effect (Johnston et al. 2001). Variogram modelling and estimation is extremely important for structural analysis and spatial interpolation (Burrough & McDonnell 1998).
(1) Ordinary kriging (OK) OK uses a geostatistical model of reality and makes the following assumption: where m represents the constant mean structure of the concentration field, eðsÞ is a smooth variation plus measurement error (both zero-mean) and D is the examining area (Denby et al. 2005;Hengl 2009). The OK estimate is a linear weighted average of the available n observations defined in Equation (3) aŝ whereẑðs 0 Þ is the OK estimate at location s 0 , zðs i Þ is the measured value in the ith observation point, i D 1, …, n, n is the number of surrounding stations from which the interpolation is computed, λ 1 , …, λ n are the weights assumed at the basis of variogram.
Kriging is an often used standard method; to apply OK, it is necessary that some essential assumptions be fulfilled, especially stationarity (i.e. the assumptions given in Equations (3) and (5) are fulfilled). For detailed information, see e.g. Lloyd (2005) and Denby et al. (2005).
(2) Simple kriging (SK) SK is mathematically the simplest, but the least general. It assumes the expectation of the random field to be known, and relies on a covariance function. However, in most applications, neither the expectation nor the covariance is known beforehand.
The practical assumptions for the application of simple Kriging are as follows: (a) Wide-sense stationarity of the field (b) The expectation is zero everywhere: mðxÞ D 0 (c) Known covariance function: cðx; yÞ D CovðZðxÞ; ZðyÞÞ SK is based on a stationary random function model ZðsÞ D m C RðsÞ (7) with s being a vector of spatial coordinates, m the stationary mean or systematic component, and R (s) the stochastic component with zero mean, constant variance and with spatial covariance C(h), the vector h representing the lag distance between two locations. The random function Z is assumed to have generated the unknown reality z, which is observed at a limited set of point locations and which needs to be interpolated. For detailed information, see e.g. Li and Heap (2008) and Knnotters et al. (2010).
(3) Universal kriging (UK) The kriging with a trend (KT) is normally called UK that was proposed by Matheron (1969). It is an extension of OK by incorporating the local trend within the neighbourhood search widow as a smoothly varying function of the coordinates. UK estimates the trend components within each search neighbourhood window and then performs SK on the corresponding residuals.
In UK, the following model assumption is considered, which is similar to Equation (7): ZðsÞ D mðsÞ C RðsÞ where m(s) is a trend component. The trend component can be defined as m(s) D E{Z(s)} and is modelled as a (local) function from the coordinates vector, whose unknown parameters are fitted from the data: Here, the functions f k are so-called base functions and must be known. The residual component RðsÞ is modelled as a stationary random function with zero mean and covariance C R ðhÞ. The value of z at an unsampled location s 0 is estimated by (4) Disjunctive kriging (DK) DK is used for the primary variable that the conventional transformations (e.g. logarithm or square-rooting) cannot yield a near-normal distribution. In DK, the primary variable is transformed into Hermite polynomials, which are a series of normally distributed subvariables that are kriged separately. DK also provides an estimate of the conditional probability that a random variable located at a point, or averaged over a block in two-dimensional space, exceeds certain thresholds (Li & Heap 2008). DK produces a nonlinear unbiased, distribution-dependent estimator with the characteristics minimum variance of errors (Burrough & McDonnell 1998).
DK assumes the model where m 1 is a known constant and f ½ZðsÞ is some arbitrary function of ZðsÞ. Notice that you can write f ½ZðsÞ D I½ZðsÞ > c t , so indicator kriging is a special case of DK. In the geostatistical analyst, you can predict either the value itself or an indicator with DK.
(5) Co-kriging (CK) CK uses information on several variable types. CK improves over kriging only when the secondary variables are better sampled than the primary variable, or more accurately reflect the real world. CK is most effective when the covariance is highly correlated with the primary variable (Hartkamp et al. 1999). The main variable of interest is Z 1 , and both autocorrelation for Z 1 and crosscorrelations betweenZ 1 and all other variable types are used to make better predictions. It is appealing to use information from other variables to help make predictions, but it comes at a price. CK requires much more estimation, including estimating the autocorrelation for each variable as well as all cross-correlations. Theoretically, you can do no worse than kriging because if there is no crosscorrelation, you can fall back on autocorrelation for Z 1 . However, each time you estimate unknown autocorrelation parameters, you introduce more variability, so the gains in precision of the predictions may not be worth the extra effort.
CK is proposed to use non-exhaustive secondary information and to explicitly account for the spatial cross-correlation between the primary and secondary variables (Goovaerts 1997): where m 1 is a known stationary mean of the primary variable, Z 1 ðx i 1 Þ is the data of the primary variable at pointi 1 , m 1 ðx i 1 Þ is the mean of samples within the search window, n 1 is the number of sampled points within the search window for point x 0 used to make the estimation, λ i 1 is the weight selected to minimize the estimation variance of the primary variable, n v is the number of secondary variables, n j is the number of jth secondary variable within the search window, λ i j is the weight assigned to i j th point of jth secondary variable, Z j ðx i j Þ is the data at i j th point of jth secondary variable, and m j ðx i j Þ is the mean of samples of jth secondary variable within the search widow . The cross-semivariance (or cross-variogram) can be estimated from data using the following equation:ĝ 12 ðhÞ D 1 2n where n is the number of pairs of sample points of variables Z 1 and Z 2 at point x i , x i C h separated by distance h (Burrough & McDonnell 1998). Cross-semivariances can increase or decrease with h depending on the correlation between the two variables, and the Cauchy-Schwartz relation must be checked to ensure a positive CK estimation variance in all circumstances (Burrough & McDonnell 1998).

Cross-validation methodology
Leave-one-out cross-validation method was used to validate the interpolation results. The method works as follows: the model predicts a complete annual mean temperature at all stations using only observations from neighbouring stations. The predicted values were compared with the actual observations of the target stations to derive cross-validation statistics (Kilibarda et al. 2014). The error indicators applied include (1)root-mean-square error (RMSE); (2) index of agreement (IOA); (3) Nash-Sutcliffe efficiency (NSE); (4) average standard error (ASE); and (5) root-mean-square standardized error (RMSSE). The criteria for selecting an optimal model are set as follows: the RMSE is smaller, the interpolated result is better; IOA (d) varies between 0 and 1, and a value of 1 indicates a perfect match, and 0 indicates no agreement at all (Willmott 1981); NSE (E) ranges from ¡1 to 1. An efficiency of 1 (E D 1) corresponds to a perfect match of modelled discharge to the observed data, an efficiency of 0 (E D 0) indicates that the model predictions are as accurate as the mean of the observed data, whereas an efficiency less than zero (E < 0) occurs when the observed mean is a better predictor than the mode (Nash & Sutcliffe 1970); ASE near RMSE; and RMSSE near 1. The above-mentioned error indicators and optimal model choosing criteria are applicable to the interpolation techniques listed in Figure 1. Additionally, all the interpolation processes other than the spline process that uses the ANUS-PLIN program are realized through the Arcgis platform: whereẐðs i Þ is the predicted value, zðs i Þ is the observed value,ŝðs i Þ is the prediction standard error for location s i , m is the mean value of the observed value andm is the mean value of the predicted value.
The aim of the ANUSPLIN package is to provide a facility for transparent analysis and interpolation of noisy multivariable data using thin-plate smoothing splines. The package supports this process by providing comprehensive statistical analyses, data diagnostics and spatially distributed standard errors. For the spline process, the ANUSPLIN program provides a number of criteria that can be used to work on an optimal PTPSS model: (1) no asterisk ( Ã ) logo; (2) signal smaller than half of the station number (signal is the number of degrees of freedom of the fitted spline, or the effective number of parameters); (3) the smallest signal-to-noise ratio (SNR); and (4) the smallest GCV in the cross-verification results (Hutchinson 2004).

Comparison of semivariance model for geostatistic methods
As above mentioned, semivariogram modelling and estimation is extremely important for structural analysis and spatial interpolation. Compare the semivariance model for different geostatistic methods to determine the best temperature interpolation scheme.

Comparison and validation of interpolation methods
First, make a lot of tests of each interpolation method and compare the error indicators of the interpolation results to find the best fitted model to each method (Table 1). Then, efforts are made to sort out the optimal interpolation method by comparing the error indicators of all of the best-fitted model in Step (1) (Table 2). Finally, use the best interpolation method above to give the distribution map of the spatial interpolation results and associated prediction error.

Comparison of time-series interpolation results between 1971 and 2000
Each year temperature-gridded surface between 1971 and 2000 was obtained by interpolating annual mean temperature data at 104 stations based on the optimal interpolation methods determined in Table 2. Temperature characteristic value (i.e. means, minimums and maximums) curves (name 'interpolated temperature curves') were extracted from above-30 temperature-gridded surfaces from 1971 to 2000. In the study, the following three approaches are also used to validate the interpolation methods by temperature variation curves over 30 years: (1) extract the means, minimums and maximums from the 1 £ 1 gridded Chinese surface temperature data ( name 'reference-gridded data' ) limited in Northeast China between 1971 and 2000, which was made of temperature characteristic value curves (name 'reference temperature curves') , and compare the correlation coefficients between the reference temperature curves and the interpolated temperature curves in the context of annual means, minimums and maximums; (2) compare the predicted error variation curves of mean temperature based on different methods (the predicted mean temperature (interpolated with different interpolation methods) minus the reference-gridded mean temperature value in each year); and (3) The best interpolated result based on each method is in bold font. ( Figure 3). It is apparent that PTPSS saw a difference between CK and OK with the former having a more allowable prediction error range compared with the latter two (Table 2). All the three scenarios have their large-error areas along the borders of the study area except the southern section (Figure 3). The large error could be caused by (1) the very limited number of meteorological observing sites along the borders; and (2) some large-error areas sitting on high elevations. For example, the western part of Chifeng claims an elevation of 1000 m or above. Both CK and PTPSS scenarios take elevation as a covariate. However, PTPSS produces a prediction error that is noticeably smaller compared with CK, especially at higher elevations, demonstrating the suitability of the PTPSS scenario. At the spatial scale, the PTPSS with two-order interpolation proves most desirable for annual mean temperature in Northeast China.

Analysis of interpolation results based on different methods between 1971 and 2000
Compare the correlation coefficients between the interpolated results and the reference-gridded data in the context of annual means, minimums and maximums. All the interpolated scenarios, except the one using GP technique, show a fine correlation between the interpolated results and the referencegridded data at 0.9 or above, with PTPSS claiming the best at 0.996. For the annual maximums, only IDW, UK and PTPSS stand at 0.9 or above, with OK, SK, DK, CK and TPS at 0.9 or above when the minimums are involved. Apparently, the PTPSS scenario renders fine annual means, minimums and maximums. Of all the interpolation scenarios, PTPSS produces the smallest error approaching zero for the annual mean temperature prediction (Figure 4) by CK and OK. A comparison with the results (figure omitted) of the studies made by He et al. (2013) and Fu et al. (2009) shows that the results unveiled by the current study share the similarity to the previous results for the period of 1971-2000, confirmed the same climate-change trend featured with a fluctuating rise of annual mean temperature, though the annual mean temperatures from the previous two studies are higher than the one from the current study, with the reference-gridded 1 £ 1 temperatures being universally higher than the current study but lower than the previous two studies. In Figure 5, it is easy to see that the PTPSS scheme has the mean temperature prediction errors that are universally smaller than 0, suggesting a systematic error against the reference-gridded data. One can correct the errors by averaging the difference between the mean temperatures based on the gridded and PTPSS interpolated temperatures. The average values, when added to the PTPSS temperatures, can generate a corrected curve that almost overlaps the reference-gridded one. At the temporal scale, the PTPSS scenario also proves able to generate better results.

Discussion
Spatial interpolation is a technique that can be used to regionalize meteorological elements. One can easily see from the current study that an interpolation technique that involves multiple variables, such as PTPSS and CK, is able to generate more reasonable results, compared with the one that has a single variable. This study made a comparison of the results derived from the above two interpolation scenarios under different combinations of covariates. Of them, the one using latitude and longitude as independent variables and elevation as a covariate stands out with desirable results. Of the interpolation scenarios using the single variable, the OK scenario brings up the highest accuracy. The results bear a similarity with the one from the study made by Atkinson and Lloyd (1998) who considers that of all the interpolation scenarios tested, the geostatistical kriging technique and the thin-plate spline technique are the most desirable tools for the interpolation of meteorological elements (Atkinson & Lloyd 1998). Comparing with other interpolation techniques that use the single variable, the kriging technique covers the spatial distribution of meteorological elements and their spatial autocorrelation, and produces the bias-free results with raised accuracy. Theoretically, the more covariates get involved in the interpolation, the higher the accuracy of the results would be. Unfortunately, the empirical interpolation would not tell the same situation. We may take three topographical factors: elevation, slope and aspect as an example. There is a strong correlation between the three. An interpolation process would ask to have the factors bearing the strongest correlation to be removed, while only retaining the factors that would have a significant impact on the interpolation results or would significantly improve the accuracy of the interpolation. In the study, the author also tried the schemes using only slope or aspect as a covariate, or both slope and aspect as covariates, and found that all the three interpolation schemes would produce very close results which were still more accurate than the interpolated results without any covariates, which suggested that adding terrain factors as covariates could improve the accuracy of interpolated results. However, when elevation is added, one would see the reduced error indicators of RMSE and ASE, suggesting that elevation is a more important factor that would affect the interpolated annual mean temperature, which is consistent with the results of the study made by Liu et al. (2008).
Additionally, most previous studies have worked on the verification of mean air temperatures, where error indicators would only be used to compare the mean air temperature so as to sort out an optimal interpolation technique. The selected optimal technique would then be used to study the dynamic interdecadal variation of temperatures. As a matter of fact, the approach to verify the accuracy of interpolation results, as such, would fail to see a full picture. Meanwhile, one cannot verify the accuracy of annualized interpolations only with a number of error indicators, as some extremes would be smoothed out in annualizing that may compromise the accuracy of the interpolated surface. A more desirable approach would be the combined verification using both error indicators and temperature characteristic values for an overall verification that involves both time and space.

Conclusion
In the study, different interpolation techniques are tested to assess the accuracy of the annual mean temperature in Northeast China. The results show that the PTPSS with two-order interpolation using latitude and longitude as independent variables and elevation as a covariate is the most desirable interpolation scheme for annual mean temperature in Northeast China. The tests come out with an RMSE prediction error at 0.34 C, an IOA at 0.998 and an NSE at 0.991. The PTPSS scenario generates a correlation coefficient at 0.9 or above for the time curves of the means, minimums and maximums of the temperature-gridded surfaces. Of all the interpolation methods tested, PTPSS scheme renders the smallest temperature prediction error approaching zero.