Spatial Interpolation Explained: Guessing Between the Samples #

Problem statement #

You have measurements at points — rainfall gauges, soil pits, boreholes, air-quality sensors — and you need a value everywhere. Interpolation fills the gap, and it always produces an answer, at every location, with no warning about how much of that answer is data and how much is assumption.

The honest way to think about it is: an interpolated surface is mostly a model, and the model's error grows with distance from the nearest sample. That growth can be measured. Taking a real 30 m digital elevation model as ground truth, sampling 400 points from it and interpolating back:

distance to nearest sample      cells   RMSE
        0 -   100 m             6,959   18.9 m
      100 -   250 m             7,286   54.0 m
      250 -   500 m             9,479  106.1 m
      500 - 1,000 m            19,802  169.8 m
    over   1,000 m             72,611  334.3 m

Same surface, same method, same colour ramp. Error varies by a factor of eighteen across it, and nothing in the output shows where.

Quick answer #

Interpolate, and return the distance to the nearest sample alongside the prediction:

import numpy as np
from scipy.spatial import cKDTree


def idw(points, values, targets, power=2, k=12):
    """Inverse distance weighting, with the support distance returned."""
    tree = cKDTree(points)
    distance, index = tree.query(targets, k=min(k, len(points)))
    distance = distance.reshape(len(targets), -1)
    index = index.reshape(len(targets), -1)

    exact = distance[:, 0] < 1e-9                    # a target on a sample
    weights = 1.0 / np.maximum(distance, 1e-9) ** power
    predicted = (weights * values[index]).sum(1) / weights.sum(1)
    predicted[exact] = values[index[exact, 0]]

    return predicted, distance[:, 0]                 # value, and its support

The second return value is not a diagnostic extra. It is the only thing that separates a prediction you can defend from one you cannot.

Step-by-step solution #

1. Decide whether the variable is interpolatable at all #

Interpolation assumes spatial continuity: nearby places are more alike than distant ones. That is true for elevation, temperature, soil pH and groundwater level. It is false for land ownership, species presence at fine scale, and anything dominated by a boundary.

The test is empirical — fit a variogram and look at whether variance rises with distance and levels off. If it is flat from the shortest lag, there is no spatial structure to exploit and the best prediction everywhere is the mean. See The variogram explained.

2. Choose a method, knowing they differ less than the literature suggests #

Measured on the same 500 sample points against 12,000 held-out cells of a real DEM:

kriging (fitted spherical model)   40.2 m
TIN, linear                        46.7 m
IDW, power 2                       54.0 m
nearest neighbour                  65.1 m

Kriging is genuinely the best here, by about 25% over IDW. That is a real gain and it is smaller than the gain from adding samples: going from 500 to 1,000 random points cut IDW's error from 58.4 m to 43.3 m — a 26% improvement from data rather than from method.

More samples beat a better interpolator. If you can influence the survey, spend there.

3. Tune the parameters, but expect a shallow optimum #

IDW's power, at 500 samples:

power  0.5     1      2      3      4      6
RMSE  72.2   64.4   54.0   50.7   50.9   53.2

And its neighbour count:

k      1      3      6     12     24     48    500
RMSE  64.5  52.5   51.5   54.0   59.6   66.9   90.1

Both have a broad minimum. Power 3 beats power 2 by 6%; using all 500 points instead of the nearest 6 costs 75%. The parameter that matters is k, and the failure mode is using too many neighbours, which drags every prediction towards the global mean.

4. Match the output cell size to the sample spacing #

With 500 samples over 105 km², the mean spacing is sqrt(area/n) = 457 m. Interpolating to finer and finer grids:

cell size    output cells    RMSE
     10 m      1,045,233    54.01 m
     30 m        116,137    53.96 m
    120 m          7,216    53.89 m
    250 m          1,638    54.17 m

A 10 m grid has ninety times the cells of a 120 m grid and is no more accurate — very slightly worse. The extra cells are interpolation, not information. See How to choose a cell size for an interpolated surface.

5. Validate against held-out data, and know the validation is optimistic #

Leave-one-out cross-validation on the same 500 points reported an RMSE of 60.8 m. The true error across every cell of the surface was 54.0 m.

Here LOO was 13% pessimistic, because removing a point from a random sample leaves a slightly larger gap than the average gap. Under clustered sampling it goes the other way and becomes dramatically optimistic, because each held-out point still has near neighbours from its own cluster.

Code examples #

Example 1 — interpolate and report support in one pass #

import numpy as np
from scipy.spatial import cKDTree


def interpolate_with_support(points, values, targets, power=2, k=12,
                             max_distance=None):
    """IDW plus the distance to the nearest observation for every prediction."""
    tree = cKDTree(points)
    distance, index = tree.query(targets, k=min(k, len(points)))
    distance = distance.reshape(len(targets), -1)
    index = index.reshape(len(targets), -1)

    weights = 1.0 / np.maximum(distance, 1e-9) ** power
    predicted = (weights * values[index]).sum(1) / weights.sum(1)
    predicted[distance[:, 0] < 1e-9] = values[index[distance[:, 0] < 1e-9, 0]]

    nearest = distance[:, 0]
    if max_distance is not None:
        predicted = np.where(nearest <= max_distance, predicted, np.nan)

    for lo, hi in ((0, 100), (100, 250), (250, 500), (500, 1000), (1000, np.inf)):
        band = (nearest >= lo) & (nearest < hi)
        if band.sum():
            label = f"{lo:>5.0f}-{'inf' if np.isinf(hi) else int(hi):>5}"
            print(f"  {label} m: {int(band.sum()):7,} cells")
    return predicted, nearest

Setting max_distance to something defensible — roughly the variogram range, or twice the mean sample spacing — turns "the model extrapolated confidently into a data desert" into a visible hole in the map.

Example 2 — comparing interpolators on your own data #

import numpy as np
from scipy.interpolate import griddata
from scipy.spatial import cKDTree


def compare_interpolators(points, values, test_points, test_values):
    """Which method actually wins on this dataset? Measure, do not assume."""
    def rmse(pred, truth):
        ok = np.isfinite(pred)
        return float(np.sqrt(np.mean((pred[ok] - truth[ok]) ** 2)))

    results = {
        "nearest": griddata(points, values, test_points, method="nearest"),
        "tin_linear": griddata(points, values, test_points, method="linear"),
        "idw_p2": idw(points, values, test_points, power=2, k=12)[0],
        "idw_p3": idw(points, values, test_points, power=3, k=12)[0],
    }

    for name, pred in sorted(results.items(),
                             key=lambda kv: rmse(kv[1], test_values)):
        coverage = float(np.isfinite(pred).mean())
        print(f"  {name:12} RMSE {rmse(pred, test_values):6.2f}  "
              f"coverage {coverage:5.1%}")
    return results

  idw_p3       RMSE  50.74  coverage 100.0%
  idw_p2       RMSE  53.96  coverage 100.0%
  nearest      RMSE  65.05  coverage 100.0%
  tin_linear   RMSE  46.72  coverage  97.3%

Note the coverage column. Linear TIN interpolation only produces values inside the convex hull of the samples — 97.3% of cells at 500 points, and 85.2% at 100 points. That is a feature: it refuses to extrapolate. Comparing its RMSE against methods that do extrapolate is not comparing like with like.

Example 3 — a sample-spacing report before you interpolate anything #

import numpy as np
from scipy.spatial import cKDTree


def sampling_report(points, bounds):
    """Is this sample dense enough, and is it evenly spread?"""
    left, bottom, right, top = bounds
    area = (right - left) * (top - bottom)
    n = len(points)

    ideal = np.sqrt(area / n)                     # spacing of a perfect grid
    tree = cKDTree(points)
    nn, _ = tree.query(points, k=2)
    observed = float(np.median(nn[:, 1]))         # median sample-to-sample gap

    print(f"  {n} samples over {area / 1e6:.1f} km²")
    print(f"  ideal grid spacing        {ideal:7.0f} m")
    print(f"  median nearest-neighbour  {observed:7.0f} m")
    print(f"  clustering ratio          {ideal / observed:7.2f}"
          f"  ({'clustered' if ideal / observed > 1.5 else 'well spread'})")
    return {"ideal_spacing": ideal, "observed_spacing": observed}

A clustering ratio well above 1 means the samples are closer to each other than a uniform design would put them — so they cover less of the area than their count suggests, and cross-validation on them will be optimistic.

Explanation #

Why every interpolator is a weighted average #

Nearest neighbour, IDW, TIN, splines and kriging all produce sum(w_i * z_i) for some weights that depend on geometry. They differ only in how the weights are chosen:

• Nearest neighbour gives weight 1 to the closest sample and 0 to the rest.
• IDW uses 1/d^p, decided by you.
• TIN uses barycentric coordinates within a triangle — three non-zero weights.
• Kriging solves for the weights that minimise expected squared error under a fitted model of spatial correlation.

That last one is why kriging usually wins: its weights come from measured spatial structure rather than from an assumption you typed in. It is also why it is slower, needs a variogram, and can fail when the variogram fit is poor.

Why error grows with distance, and why the hull is the wrong boundary #

The measured error bands at the top of this page are the honest summary of interpolation: 18.9 m within 100 m of a sample, 334.3 m beyond a kilometre.

A common heuristic is to clip the surface to the convex hull of the samples. Measured on a clustered sample of 400 points, that hull covered 59.0% of the study area, and the error inside it was 268 m against 286 m outside — barely different.

The hull is a poor proxy because a point can sit inside the hull and still be a kilometre from any sample, in the middle of a gap between clusters. Distance to the nearest sample is what predicts error, and it is what you should mask on.

Why more data beats a better method #

Doubling the sample from 500 to 1,000 improved IDW from 58.4 m to 43.3 m — 26%. Switching from IDW to kriging on 500 points improved it from 54.0 m to 40.2 m — 25%.

They are comparable in size, and only one of them is free. But the sampling gain compounds: 2,000 points brought IDW to 28.8 m, better than kriging on 500. No interpolator recovers information that was never measured.

Why the surface looks better than it is #

An interpolated surface is smooth, continuous and visually convincing everywhere. It carries no visual signal of its own uncertainty — the part built from thirty nearby samples looks identical to the part extrapolated across a five-kilometre gap.

That is the argument for always producing two rasters: the prediction, and the distance to the nearest sample (or the kriging variance, with the caveats in How to fit a variogram and krige a surface). Publishing the first without the second invites everyone downstream to trust the whole map equally.

Edge cases or notes #

• Duplicate coordinates break most solvers. Two samples at the same location give a singular kriging matrix; aggregate them first.
• IDW is an exact interpolator: every sample is reproduced exactly, which is what produces bullseyes. See My IDW surface is full of bullseyes.
• Use a projected CRS. Distance weighting in degrees weights latitude and longitude differently everywhere except the equator.
• Barnes and Cressman are IDW variants with a distance cutoff; the cutoff matters more than the kernel.
• Anisotropy is common — rainfall along a valley, contamination down a gradient. Isotropic IDW cannot represent it; kriging can.
• Interpolating a rate or a ratio is usually wrong. Interpolate the numerator and denominator separately, then divide.
• A finer output grid never adds information. Ninety times more cells gave a 0.1 m accuracy change here.
• Do not interpolate across a barrier. A cliff, a fault or a coastline breaks the continuity assumption the whole method rests on.

Internal links #

• The variogram explained — measuring whether spatial structure exists, and how far it reaches
• IDW, kriging, splines or TIN? Choosing an interpolator — picking between the methods
• How to interpolate points to a grid with IDW in Python — the practical version
• How to fit a variogram and krige a surface in Python — the method that won here
• How to cross-validate an interpolated surface — why LOO reported 60.8 m for a 54.0 m surface
• How to choose a cell size for an interpolated surface — the 457 m spacing rule
• My interpolated surface extends far beyond the data — masking on support distance
• Sample design explained: where to measure — why more samples beat a better method

FAQ #

What is spatial interpolation? #

Estimating a value at unsampled locations from measurements at sampled ones, assuming nearby places are more alike than distant ones.

Which interpolation method is best? #

On the data measured here, kriging (40.2 m RMSE) beat linear TIN (46.7 m), IDW (54.0 m) and nearest neighbour (65.1 m). The gap between methods was about the same as the gap from doubling the number of samples.

How accurate is an interpolated surface? #

It depends on distance to the nearest sample, and it varies hugely across one map — from 18.9 m within 100 m of a sample to 334.3 m beyond a kilometre in the measurement above.

Should I clip the surface to the convex hull? #

The hull is a weak proxy for support. On a clustered sample it covered 59% of the area with barely lower error inside than outside. Mask on distance to the nearest sample instead.

What cell size should I interpolate to? #

Around the mean sample spacing, sqrt(area/n). Finer grids cost storage and give no accuracy — 90× more cells changed RMSE by 0.1 m here.

Does interpolation work for any variable? #

Only for variables with spatial continuity. Fit a variogram first; if variance does not rise with distance, there is no structure and the mean is your best estimate everywhere.

How many samples do I need? #

More than you think, and more helps more than a better algorithm: 500 → 1,000 → 2,000 points cut IDW's error from 58.4 m to 43.3 m to 28.8 m.