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 of Snowdonia, tested against 12,000 held-out cells (Spatial interpolation explained):

kriging (fitted spherical model)   41.3 m RMSE
TIN, linear                        49.1 m
IDW, power 2                       56.1 m
nearest neighbour                  65.4 m

Kriging beats IDW by 26%. Meanwhile, doubling the sample from 500 to 1,000 points improved IDW from 55.4 m to 40.6 m — a 27% gain, from the same method. And in the bake-off in Example 1, a thin-plate spline beat all four at 36.8 m: elevation is the smooth kind of field splines suit.

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

No structure, categories, IDW, TIN or kriging

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  75.3   67.2   56.1   52.5   52.4   54.7

k      1      3      6     12     24     48    500
RMSE  65.4  53.7   52.5   56.1   62.2   69.9   97.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 85% 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        94.7%
    250        94.5%
    500        97.4%
  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, in How to fit a variogram and krige a surface:

mean kriging variance     4,450 m²  ->  predicted sd 66.7 m
actual RMSE                                          40.9 m

The kriging variance overstated the real error by 63%. 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.

Choosing by measurement, and a hybrid

The data used on this page

The surface is real: the Copernicus DEM GLO-30 over a 10 × 10 km window centred on the summit of Snowdon (Yr Wyddfa), warped to a 30 m grid in EPSG:32630 (x 422,880–432,870, y 5,875,440–5,885,430; 333 × 333 cells, 99.8 km²). The 500 samples are cell centres drawn at random with a fixed seed; 12,000 further cells, none of them samples, are held out as test_points/test_values. grid_points/truth are every cell centre and its elevation, and truth_fn reads the DEM bilinearly at any location, which the error tables score against. idw is the function from Spatial interpolation explained, and surface is its power-2, 12-neighbour surface on that 30 m grid. The comparison tables quoted in the steps were measured on these same samples in Spatial interpolation explained and How to fit a variogram and krige a surface.

import os

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy.interpolate import RegularGridInterpolator

bounds = (422_880.0, 5_875_440.0, 432_870.0, 5_885_430.0)
transform = from_origin(bounds[0], bounds[3], 30, 30)
shape = (333, 333)
if not os.path.exists("snowdon_dem_utm30.tif"):
    COP = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{0}_00_DEM/"
           "Copernicus_DSM_COG_10_N53_00_{0}_00_DEM.tif")
    dem = np.full(shape, np.nan, "float32")
    for tile in ("W005", "W004"):
        part = np.full(shape, np.nan, "float32")
        with rasterio.open(COP.format(tile)) as src:
            reproject(rasterio.band(src, 1), part, dst_transform=transform, dst_crs="EPSG:32630",
                      dst_nodata=np.nan, resampling=Resampling.bilinear)
        dem = np.where(np.isnan(dem), part, dem)
    with rasterio.open("snowdon_dem_utm30.tif", "w", driver="GTiff", width=shape[1], height=shape[0], count=1,
                       dtype="float32", crs="EPSG:32630", transform=transform, nodata=np.nan) as dst:
        dst.write(dem, 1)
dem = rasterio.open("snowdon_dem_utm30.tif").read(1)

xs = bounds[0] + (np.arange(shape[1]) + 0.5) * 30
ys = bounds[3] - (np.arange(shape[0]) + 0.5) * 30
gx, gy = np.meshgrid(xs, ys)
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
truth = dem.ravel().astype("float64")

rng = np.random.default_rng(0)
order = rng.permutation(truth.size)
sample_idx, test_idx = order[:500], order[500:12_500]
points, values = grid_points[sample_idx], truth[sample_idx]
test_points, test_values = grid_points[test_idx], truth[test_idx]
_bilinear = RegularGridInterpolator((ys[::-1], xs), dem[::-1], bounds_error=False, fill_value=None)
truth_fn = lambda xy: _bilinear(np.column_stack([xy[:, 1], xy[:, 0]]))
from scipy.spatial import cKDTree


def idw(points, values, targets, power=2, k=12):
    tree = cKDTree(points)
    distance, index = tree.query(targets, k=min(k, len(points)))
    distance, index = distance.reshape(len(targets), -1), index.reshape(len(targets), -1)
    weights = 1.0 / np.maximum(distance, 1e-9) ** power
    predicted = (weights * values[index]).sum(1) / weights.sum(1)
    exact = distance[:, 0] < 1e-9
    predicted[exact] = values[index[exact, 0]]
    return predicted, distance[:, 0]


surface = idw(points, values, grid_points)[0].reshape(shape)
print(f"DEM {np.nanmin(dem):.0f} to {np.nanmax(dem):.0f} m; {len(points)} samples, {len(test_points):,} held out")

DEM 52 to 1072 m; 500 samples, 12,000 held out

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


_ = bake_off(points, values, test_points, test_values)

  rbf_thin_plate   RMSE   36.79   coverage 100.0%
  tin_cubic        RMSE   40.18   coverage  97.8%
  tin_linear       RMSE   49.06   coverage  97.8%
  idw_p3           RMSE   52.49   coverage 100.0%
  idw_p2           RMSE   56.08   coverage 100.0%
  nearest          RMSE   65.41   coverage 100.0%

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]


_ = tune_neighbours(points, values)

  k=  3  LOO RMSE   52.33
  k=  6  LOO RMSE   51.09
  k=  9  LOO RMSE   54.02
  k= 12  LOO RMSE   55.84
  k= 18  LOO RMSE   58.71
  k= 24  LOO RMSE   61.96
  k= 36  LOO RMSE   65.89
  best k = 6

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


_ = local_mean_fallback(points, values, grid_points, radius=500)      # 500 m: about the sample spacing

  97.4% interpolated, 2.6% local mean, 0.0% masked

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^2$ produce a confident number four kilometres from anything.

Why the methods mostly agree

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 there all of them are wrong: error grows with distance from the nearest sample whatever the method, which is the subject of spatial interpolation explained. The choice of method matters most in between — and the next section is where kriging earns its keep.

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 26% 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.4 m, only 17% behind IDW at 56.1 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.

Neighbour count before power, and spline overshoot

• Categorical data has no average. Nearest neighbour only.
• Tune k before tuning the power. Neighbour count cost 85% here; power cost 6%.
• A large k drags 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 63% 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, a thin-plate spline at 36.8 m RMSE, then kriging at about 41 m, linear TIN at 49.1 m, IDW at 56.1 m and nearest neighbour at 65.4 m. The ranking depends on the variable — splines win on smooth fields and overshoot on stepped ones — so run the bake-off on your own data.

Is kriging worth the extra effort?

It bought 26% 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 85%.

How many neighbours should IDW use?

Six to twelve. Using all 500 samples raised RMSE from 52.5 m to 97.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 63% 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 — 94.7% of the held-out cells at 100 points, 97.4% at 500. That is the method refusing to extrapolate.