# Assessing error of a spatial interpolation algorithm

I have a set of 3D points. The points have three components $x$, $y$, $z$. You can think of these points as the surveyor points that one collects from measuring a terrain for GIS purposes.

I have a few algorithms that take in this set of 3D points as the reference and also a list of 2D points $p_{out}(x,y)$: the algorithms will then compute the $z_{out}$ for each of the $p_{out}(x,y)$.

Now, given that I have a few algorithms for this, is there anyway I can evaluate and compare the correctness of these algorithms? Is there any kind of "error index" formula that allows me to know which algorithms are better at interpolation?

Or, Is there any kind of formula that allows me to know which algorithms are better at generating the 3D surface? The two most intuitive criteria for the evaluation are:

1. How near the interpolated points are to the given input points

2. How smooth is the surface

I have no idea how to quantitatively define the above two, what do you think?

-
I changed the "3D" in the title to "spatial" because this is really 2D interpolation (for most terrain, anyway--excluding tunnels, arches, and overhanging cliffs). –  whuber Sep 30 '11 at 16:02
Could Procrustes analysis help your point 1 (comparing how close the model shape resembles the input shape)? –  ttnphns Oct 1 '11 at 7:04

One option may be to split the original data into two subsets: one that will be used in interpolating values and one that will be used to validate the interpolation results. The error is then estimated by comparing interpolated values at the validation point locations with the actual validation point values. Note that the appropriateness of this approach is largely driven by the sample point density and distribution vis-à-vis the type and scale of the underlying process you are attempting to model.

Edit: This is an expansion of the original answer following @whuber’s comments.

As noted by @whuber, one disadvantage with the aforementioned technique is the degradation in the interpolation’s quality when a subset of the sampled sites are removed. A solution to this problem is described in Maciej Timczak’s 1998 paper. The author applies a cross-validation technique to estimate the optimal interpolation parameters. He then uses a jackknife method to estimate the predicted value and uncertainty at an unsampled site $Z_J$. A brief summary of the techniques described in the paper with a simple example in R follows.

In the cross-validation (leave-one-out) method, one data point $s_i$ is removed from the point data set $S$ and its interpolated value is computed using all other $(n-1)$ points of $S$. The interpolated value is then compared with the actual value $s_i$. This process is repeated with all other data points from $S$. The performance of the interpolator is evaluated via the root-mean of squared residuals (RMSE). The RMSE can be computed for different interpolation parameters (or different interpolators) then compared. The interpolator with the lowest RMSE is usually desired. $$RMSE=\sqrt{\frac{\sum_{i=1}^n (Z_{i(int)} -Z_{i})^2}{n}}$$

Once the interpolation technique is chosen, the jackknife technique is used to estimate the unknown value $Z_j$ at an unsampled location $s_j$ along with its confidence interval. The method, as implemented by the author, involves first interpolating the unknown value $Z_j$ at unsampled location $s_j$ using all sample points from $S$, then interpolating a pseudo-value $Z_i^*$ using $(n-1)$ points from the sample dataset $S$ (i.e. $Z_i^*$ is computed once for each omitted point sample $s_i$): $$Z_{i}^* = n Z_{all} - (n-1) Z_{-i}$$ where i = 1,2,..,n The jackknifed estimator of $Z_j$ at $s_j$ is computed by averaging all of the pseudo-values following: $$Z_{J} =\frac{ \sum_{i=1}^{n} Z_{i}^*}{n}$$ A confidence interval is then computed as follows: $$\sigma_J=\sqrt{\frac{1}{n(n-1)}\sum_{i=1}^{n}(Z_i^*-Z_J)^2}$$ The estimated value is therefore $Z_j \pm- t_{(\alpha/2,n-1)}\sigma_J$

A simple R example follows:

X <- c(100.0,19.4,9.0,64.4,39.4,50.7,99.0,44.4,82.5,55.9,
56.2,54.0,14.9,54.8,35.5,34.6,15.2,32.0,23.8,87.4,
49.4,77.9,63.9,14.8,5.9,45.3,95.6,10.3,59.5,47.2,
26.7,46.5,41.3,62.9,34.2,3.7,57.7,78.5,73.1,28.3,
13.1,49.4,24.2,99.2,76.3,93.2,71.6,28.8,49.4,94.0,
84.4,0.0,90.3,48.4,44.8,5.1,29.8,27.7,93.8,25.6)
Y <- c(0.0,1.9,3.2,6.0,12.4,13.3,13.7,15.3,15.6,
16.0,18.0,22.3,22.9,23.0,24.3,26.3,26.6,27.4,
31.6,33.1,33.8,35.0,35.2,42.0,44.9,45.3,45.8,
48.8,50.5,58.3,60.2,60.8,60.8,61.5,64.4,64.6,
65.5,69.2,69.3,69.4,71.2,73.3,78.5,80.1,83.6,
84.7,84.7,91.1,92.1,92.7,93.0,93.7,93.8,95.1,
96.0,96.4,96.7,97.6,99.8,100.0)
Z <- c(209,478,424,817,866,720,327,833,731,1488,562,868,318,496,488,
1146,369,735,593,778,771,304,538,669,368,474,391,346,872,556,
348,765,779,809,357,720,416,544,338,560,455,555,340,307,589,
280,745,452,1116,442,659,343,385,655,828,490,425,665,276,333)
library(akima)
n <- 60 # Number of points in S
#
# Cross validation (leave-one-out) method to compare interpolators
#
P.spl <- vector()
sum.dif2 <- numeric()
for (i in 1:n){
P.spl[i] <- interpp(X[-i],Y[-i],Z[-i],X[i],Y[i],linear=F,extrap=T)$z sum.dif2 <- (P.spl[i] - Z[i])^2 } rmse = sqrt(sum.dif2 / n) rmse # # Jackknife to estimate the uncertainty in the interpolated value P.spl <- vector() Zi <- numeric() jx <- 40 # X coordinate of parameter Zj to be estimated at sj jy <- 40 # Y coordinate of parameter Zj to be estimated at sj Zall <- interpp(X,Y,Z,jx,jy,linear=T)$z  # Interp value from all si
for (i in 1:n){
Z1 <- interpp(X[-i],Y[-i],Z[-i],jx,jy,linear=T)$z # Interp value from si-1 # Calculated pseudo-value Z at j Zi[i] <- n * Zall - (n-1) * Z1 } # # Jackknifed estimator of parameter Z at location j # Zj <- sum(Zi) / n # Estimated value Zj # Estimated standard error sig.j <- sqrt(1/(n*(n-1)) * sum((Zi-Zj)^2)) # The confidence interval on Zj alpha <- 1 - .05 / 2 t.value <- qt(alpha,n-1) ci <- t.value * sig.j # Confidence interval for Zj  - Sometimes people do this, but it's problematic. The issue is that the quality of any spatial interpolation relies on spatial correlation, which tends to decrease with distance ("Tobler's First Law of Geography"). If you split randomly, you tend to increase average distances (by about$\sqrt{2}\$ on average), degrading the quality. If you split spatially (e.g., east versus west), you introduce a confounder (because the correlation of terrain in one half could differ from that in the other half). That's why jackknifing (leave-one-out cross-validation) is usually recommended here. –  whuber Sep 30 '11 at 17:40
+1 Thank you for the detailed account of this method. –  whuber Oct 13 '11 at 18:34