IDW, Kriging, Splines or TIN? Choosing an Interpolator
Problem statement
Every GIS offers four or five interpolation methods and no guidance about which one to use. The literature suggests the choice is consequential; the measurements suggest it is less consequential than almost everything else you could do instead.
On 500 sample points taken from a real 30 m elevation model, tested against 12,000 held-out cells:
kriging (fitted spherical model) 40.2 m RMSE
TIN, linear 46.7 m
IDW, power 2 54.0 m
nearest neighbour 65.1 m
Kriging wins by 25% over IDW. Meanwhile, doubling the sample from 500 to 1,000 points improved IDW from 58.4 m to 43.3 m β a 26% gain, from the same method.
The method choice is worth roughly one doubling of the sample. Both are worth having; only one of them is under your control after the survey is done.
Quick answer
| method | use it when | cost |
|---|---|---|
| Nearest neighbour | the variable is categorical, or you need exact blockiness | worst accuracy for continuous data |
| IDW | you need something defensible in ten lines, samples are dense and even | bullseyes; no uncertainty estimate |
| TIN / linear | you want a surface that refuses to extrapolate | no values outside the convex hull |
| Kriging | the variogram shows real structure and you need uncertainty | needs a variogram; slow; variance is not a confidence interval |
| Splines / RBF | the surface is genuinely smooth (a physical field) | can overshoot wildly near sharp changes |
Step-by-step solution
1. Rule out interpolation entirely if there is no spatial structure
Fit a variogram first. If semivariance is flat from the shortest lag, nearby samples are no more alike than distant ones and the best estimate everywhere is the mean. No method fixes that.
2. Categorical variable? Use nearest neighbour and stop
Soil class, land cover, geological unit β these have no meaningful average. Nearest neighbour (a Voronoi tessellation) is the only defensible choice, and the blockiness is honest: it says "this is the nearest observation" rather than pretending to a smooth transition that was never measured.
3. Need a number today? Use IDW, with k set properly
IDW is a weighted average with weights 1/d^p. Its parameters have a broad optimum in one and a sharp cliff in the other. Measured at 500 samples:
power 0.5 1 2 3 4 6
RMSE 72.2 64.4 54.0 50.7 50.9 53.2
k 1 3 6 12 24 48 500
RMSE 64.5 52.5 51.5 54.0 59.6 66.9 90.1
Power barely matters between 2 and 4. Neighbour count matters a great deal, and the failure is on the high side: using all 500 points cost 75% accuracy against using the nearest 6, because distant samples drag every prediction towards the global mean.
Default to power 2 and k between 6 and 12, and tune k by cross-validation before touching the power.
4. Want a surface that will not extrapolate? Use a TIN
Linear interpolation over a Delaunay triangulation produces values only inside the convex hull of the samples. Measured coverage:
samples TIN coverage of the study area
100 85.2%
250 93.5%
500 97.3%
1,000 99.0%
Everything outside is NaN. That is often exactly what you want β see My interpolated surface extends far beyond the data. It also means TIN's RMSE is not comparable with methods that fill the whole area, because it declined to answer where the answer was hardest.
5. Need uncertainty, or the best accuracy? Use kriging
Kriging derives its weights from the fitted variogram rather than from a typed-in exponent, which is why it won here. It also returns a variance per prediction.
Read that variance carefully. Measured on the same 12,000 held-out cells:
mean kriging variance 4,536 mΒ² -> predicted sd 67.3 m
actual RMSE 40.2 m
The kriging variance overstated the real error by 68%. It depends only on the geometry of the samples and the fitted variogram, never on their values β so it is a good relative map of where the surface is weak, and a poor absolute confidence interval.
Code examples
Example 1 β pick the method by measuring, not by reading
import numpy as np
from scipy.interpolate import griddata, RBFInterpolator
from scipy.spatial import cKDTree
def bake_off(points, values, test_points, test_values, k=12):
"""Run every candidate on your own data and rank them."""
def rmse(pred):
ok = np.isfinite(pred)
return float(np.sqrt(np.mean((pred[ok] - test_values[ok]) ** 2))), \
float(ok.mean())
tree = cKDTree(points)
def idw(power):
d, i = tree.query(test_points, k=min(k, len(points)))
d = d.reshape(len(test_points), -1); i = i.reshape(len(test_points), -1)
w = 1.0 / np.maximum(d, 1e-9) ** power
return (w * values[i]).sum(1) / w.sum(1)
candidates = {
"nearest": griddata(points, values, test_points, method="nearest"),
"tin_linear": griddata(points, values, test_points, method="linear"),
"tin_cubic": griddata(points, values, test_points, method="cubic"),
"idw_p2": idw(2),
"idw_p3": idw(3),
"rbf_thin_plate": RBFInterpolator(
points, values, neighbors=min(50, len(points)),
kernel="thin_plate_spline")(test_points),
}
rows = [(name, *rmse(pred)) for name, pred in candidates.items()]
for name, err, cov in sorted(rows, key=lambda r: r[1]):
print(f" {name:16} RMSE {err:7.2f} coverage {cov:6.1%}")
return candidates
Always print the coverage column beside the error. A method that answers 85% of the map with a low error has not beaten one that answered all of it.
Example 2 β tuning k properly, which is the parameter that matters
import numpy as np
from scipy.spatial import cKDTree
def tune_neighbours(points, values, candidates=(3, 6, 9, 12, 18, 24, 36)):
"""Leave-one-out error against neighbour count."""
tree = cKDTree(points)
best = None
for k in candidates:
# k+1 because the nearest neighbour of a sample is itself
d, i = tree.query(points, k=min(k + 1, len(points)))
d, i = d[:, 1:], i[:, 1:]
w = 1.0 / np.maximum(d, 1e-9) ** 2
pred = (w * values[i]).sum(1) / w.sum(1)
err = float(np.sqrt(np.mean((pred - values) ** 2)))
print(f" k={k:3d} LOO RMSE {err:7.2f}")
if best is None or err < best[1]:
best = (k, err)
print(f" best k = {best[0]}")
return best[0]
Excluding the point itself is the whole trick. Forget the [:, 1:] slice and every prediction is the sample's own value, LOO RMSE is zero, and the tuning silently picks the largest k.
Example 3 β a hybrid that is often the right answer
import numpy as np
from scipy.spatial import cKDTree
def local_mean_fallback(points, values, targets, radius, k=12, power=2):
"""IDW inside the variogram range, the local mean beyond it, NaN far out."""
tree = cKDTree(points)
d, i = tree.query(targets, k=min(k, len(points)))
d = d.reshape(len(targets), -1); i = i.reshape(len(targets), -1)
w = 1.0 / np.maximum(d, 1e-9) ** power
predicted = (w * values[i]).sum(1) / w.sum(1)
nearest = d[:, 0]
inside = nearest <= radius
far = nearest > radius * 3
# beyond the range, distance weighting is meaningless β average instead
predicted[~inside] = values[i[~inside]].mean(axis=1)
predicted[far] = np.nan
print(f" {inside.mean():.1%} interpolated, "
f"{(~inside & ~far).mean():.1%} local mean, {far.mean():.1%} masked")
return predicted
Beyond the variogram range, distance carries no information, so weighting by it is theatre. Switching to a plain local mean and then to NaN is more honest than letting 1/dΒ² produce a confident number four kilometres from anything.
Explanation
Why they all agree more than the literature suggests
All five methods produce sum(w_i * z_i) for weights that depend only on geometry. Where samples are dense, every reasonable weighting scheme puts almost all the weight on the same few nearby points, so they converge.
Where samples are sparse, they diverge β and where samples are sparse, all of them are wrong. Measured across distance bands from the nearest sample, IDW's error went from 18.9 m within 100 m to 334.3 m beyond a kilometre. A better interpolator moves the 334 m to perhaps 300 m. Another survey day moves it to 50 m.
Why kriging actually wins
Kriging solves for weights that minimise expected squared error under the fitted variogram. Three consequences:
- Declustering. Two samples close together carry partly redundant information, and kriging downweights them jointly. IDW counts both at full weight.
- A principled search radius. The variogram range says where correlation ends; kriging's weights fall to near zero there naturally.
- Anisotropy. If the structure is directional, the variogram model can say so and the weights follow.
The 25% gain measured here is a fair expectation for a well-behaved continuous variable with a clean variogram. It is not free: you must fit and inspect the variogram, and a bad fit produces a worse surface than IDW.
Why splines are the risky choice
Radial basis functions and thin-plate splines fit a smooth function through the samples exactly. Where data are dense and the underlying field is genuinely smooth β a pressure field, a geoid β they are excellent.
Where the field has a step, they overshoot: to pass through both sides of a cliff smoothly, the fitted surface must swing beyond both. The result is elevations above the summit and below the valley floor, with no warning. Always check the output range against the input range; if the surface exceeds the sample extremes, a spline has overshot.
Why nearest neighbour deserves more respect than it gets
It came last here at 65.1 m, only 20% behind IDW at 54.0 m, for a fraction of the thinking. And it has two properties nothing else has: it never invents a value that was not measured, and it works for categorical data.
For a first look at an unfamiliar dataset, a Voronoi map is more honest than a smooth surface. Its blockiness shows you the sampling design, which a smoothed surface hides β and the sampling design is the thing that actually determines your accuracy.
Edge cases or notes
- Categorical data has no average. Nearest neighbour only.
- Tune
kbefore tuning the power. Neighbour count cost 75% here; power cost 6%. - A large
kdrags predictions to the global mean and flattens the map. - TIN's error is not comparable with methods that fill the whole area.
- Kriging variance is geometry, not accuracy. It overstated real error by 68% here.
- Check splines for overshoot against the sample minimum and maximum.
- Duplicate points break kriging (singular matrix). Aggregate first.
- Use a projected CRS β every one of these methods weights by distance.
- Anisotropy is invisible to IDW. If structure is directional, only kriging can use that.
Internal links
- Spatial interpolation explained β why the method matters less than the sampling
- The variogram explained β the check that comes before the choice
- How to interpolate points to a grid with IDW in Python β the ten-line option
- How to fit a variogram and krige a surface in Python β the accurate option
- How to build a TIN and interpolate elevations β the option that refuses to extrapolate
- How to cross-validate an interpolated surface β how to run the bake-off honestly
- My IDW surface is full of bullseyes β IDW's signature artefact
- Sample design explained: where to measure β the lever that beats the method
FAQ
Which interpolation method is most accurate?
On the data measured here, kriging at 40.2 m RMSE, then linear TIN at 46.7 m, IDW at 54.0 m and nearest neighbour at 65.1 m. The ranking depends on the variable; run the bake-off on your own data.
Is kriging worth the extra effort?
It bought 25% over IDW here β about the same as doubling the sample size. Worth it when you also need the uncertainty map, or when the field is anisotropic.
What IDW power should I use?
Between 2 and 4; the difference was 6% here. Spend your tuning effort on the neighbour count instead, where the difference was 75%.
How many neighbours should IDW use?
Six to twelve. Using all 500 samples raised RMSE from 51.5 m to 90.1 m, because distant points pull every prediction towards the global mean.
Can I use kriging variance as a confidence interval?
Not directly. It depends only on sample geometry and the variogram, not on the data values, and it overstated the true error by 68% here. Use it as a relative map of where the surface is weak.
When should I use a spline?
When the underlying field is genuinely smooth and the samples are dense. Check for overshoot: if the surface exceeds the sample minimum or maximum, it has invented extremes.
Why does my TIN surface have NaN around the edges?
Because linear TIN interpolation only produces values inside the convex hull of the samples β 85.2% of the study area at 100 points, 97.3% at 500. That is the method refusing to extrapolate.