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 of Snowdonia as ground truth, sampling 377 points from it in eight clusters, and interpolating back with IDW (Example 4):

distance to nearest sample      cells   RMSE
        0 -   100 m             9,680   21.2 m
      100 -   250 m            14,020   50.3 m
      250 -   500 m            15,753  101.6 m
      500 - 1,000 m            25,193  142.2 m
    over   1,000 m             46,243  258.7 m

Same surface, same method, same colour ramp. Error varies by a factor of twelve 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.

Is it continuous, which method, which cell size

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)   41.3 m
TIN, linear                        49.1 m
IDW, power 2                       56.1 m
nearest neighbour                  65.4 m

Kriging is the best here, by about a quarter over IDW, and that gain is no bigger than the one from adding samples: going from 500 to 1,000 random points cut IDW's error from 55.4 m to 40.6 m (Example 4). How the methods differ, and when each is the right one, is choosing an interpolator.

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

And its neighbour count:

k      1      3      6     12     24     48    500
RMSE  65.4  53.7   52.5   56.1   62.2   69.9   97.1

Both have a broad minimum. Power 3 beats power 2 by 6%; using all 500 points instead of the nearest 6 costs 85%. 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 99.8 km², the mean spacing is sqrt(area/n) = 447 m. Interpolating to finer and finer grids, in the cell-size guide:

cell size    output cells    RMSE
     10 m        998,001    55.86 m
     30 m        110,889    55.99 m
    120 m          6,889    55.68 m
    250 m          1,521    56.06 m

A 10 m grid has 145 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 55.8 m. The true error across every cell of the surface was 56.0 m.

For a random sample LOO is a fair estimate — here within half a percent. Under clustered sampling it becomes optimistic, because each held-out point still has near neighbours from its own cluster while most of the map does not.

Interpolating Snowdon elevations with support reported

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. The clustered sample in Example 4 is stated there.

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]]))
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 — 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


_ = interpolate_with_support(points, values, grid_points)

      0-  100 m:  16,824 cells
    100-  250 m:  51,162 cells
    250-  500 m:  40,026 cells
    500- 1000 m:   2,877 cells

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


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

  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%

Note the coverage column. Linear TIN interpolation only produces values inside the convex hull of the samples — 97.8% of held-out cells at 500 points here, and 94.7% at 100 points (Example 4). 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
    random_gap = 0.5 * ideal                      # mean nearest-neighbour gap of a random pattern
    tree = cKDTree(points)
    nn, _ = tree.query(points, k=2)
    observed = float(np.mean(nn[:, 1]))           # mean sample-to-sample gap

    r = observed / random_gap                     # Clark-Evans ratio
    verdict = "clustered" if r < 0.8 else "evenly spread" if r > 1.2 else "about random"
    print(f"  {n} samples over {area / 1e6:.1f} km²")
    print(f"  ideal grid spacing        {ideal:7.0f} m")
    print(f"  mean nearest-neighbour    {observed:7.0f} m  (random would be {random_gap:.0f} m)")
    print(f"  Clark-Evans ratio         {r:7.2f}  ({verdict})")
    return {"ideal_spacing": ideal, "observed_spacing": observed, "clark_evans": r}


_ = sampling_report(points, bounds)

  500 samples over 99.8 km²
  ideal grid spacing            447 m
  mean nearest-neighbour        225 m  (random would be 223 m)
  Clark-Evans ratio            1.01  (about random)

A random pattern's nearest neighbours are, on average, half the ideal grid spacing apart, so the ratio compares against that rather than against the grid: about 1 for a random sample, above 1 for a planned, evenly spread one, and well below 1 for a clustered one. A clustered sample covers less of the area than its count suggests, and cross-validation on it will be optimistic. (An earlier version of this function compared against the grid spacing itself, and called every random sample "clustered".)

Example 4 — every number quoted on this page

from pykrige.ok import OrdinaryKriging
from scipy.spatial import ConvexHull, Delaunay


def rmse(pred, obs):
    ok = np.isfinite(pred)
    return float(np.sqrt(np.mean((pred[ok] - obs[ok]) ** 2)))


# a clustered design: 8 centres at random, 50 samples round each (sigma 400 m), snapped to cells
crng = np.random.default_rng(1)
centres = np.column_stack([crng.uniform(bounds[0], bounds[2], 8), crng.uniform(bounds[1], bounds[3], 8)])
raw = (centres[:, None, :] + crng.normal(0, 400, (8, 50, 2))).reshape(-1, 2)
raw = raw[(raw[:, 0] > bounds[0]) & (raw[:, 0] < bounds[2]) & (raw[:, 1] > bounds[1]) & (raw[:, 1] < bounds[3])]
cols = ((raw[:, 0] - bounds[0]) // 30).astype(int)
rows = ((bounds[3] - raw[:, 1]) // 30).astype(int)
cidx = np.unique(rows * shape[1] + cols)
cpoints, cvalues = grid_points[cidx], truth[cidx]
pred, near = idw(cpoints, cvalues, grid_points)
print(f"clustered sample: {len(cpoints)} points")
for lo, hi in ((0, 100), (100, 250), (250, 500), (500, 1000), (1000, np.inf)):
    band = (near >= lo) & (near < hi)
    print(f"  {lo:>5.0f}-{hi:>6.0f} m  {int(band.sum()):7,} cells  RMSE {rmse(pred[band], truth[band]):6.1f} m")
hull = Delaunay(cpoints[ConvexHull(cpoints).vertices])
inside = hull.find_simplex(grid_points) >= 0
print(f"  hull covers {inside.mean():.1%}; RMSE inside {rmse(pred[inside], truth[inside]):.0f} m, "
      f"outside {rmse(pred[~inside], truth[~inside]):.0f} m")
_ = sampling_report(cpoints, bounds)

from scipy.optimize import curve_fit
from scipy.spatial.distance import pdist

d, g = pdist(points), 0.5 * pdist(values[:, None], "sqeuclidean")      # the variogram guide's fit
edges = np.linspace(0, d.max() / 2, 16)
lag = [d[(d >= a) & (d < b)].mean() for a, b in zip(edges[:-1], edges[1:])]
gam = [g[(d >= a) & (d < b)].mean() for a, b in zip(edges[:-1], edges[1:])]
spherical = lambda h, n, s, r: n + s * np.where(h < r, 1.5 * h / r - 0.5 * (h / r) ** 3, 1.0)
(nug, psill, rng_m), _ = curve_fit(spherical, lag, gam, p0=[0, max(gam), d.max() / 4],
                                   bounds=([0, 0, 1], [max(gam), 3 * max(gam), d.max()]), maxfev=20000)
for label, kw in (("fitted variogram", {"variogram_parameters": {"sill": psill + nug, "range": rng_m, "nugget": nug}}),
                  ("pykrige's own fit", {})):
    ok = OrdinaryKriging(points[:, 0], points[:, 1], values, variogram_model="spherical", **kw)
    kriged, _ = ok.execute("points", test_points[:, 0], test_points[:, 1])
    psill_, range_, nugget_ = ok.variogram_model_parameters
    print(f"kriging, {label:18} nugget {nugget_:8,.0f} m²  range {range_:6,.0f} m  "
          f"RMSE {rmse(np.asarray(kriged), test_values):.1f} m")

print("power:", {p: round(rmse(idw(points, values, test_points, power=p)[0], test_values), 1)
                 for p in (0.5, 1, 2, 3, 4, 6)})
print("k:    ", {k: round(rmse(idw(points, values, test_points, k=k)[0], test_values), 1)
                 for k in (1, 3, 6, 12, 24, 48, 500)})

loo = []
for j in range(len(points)):
    keep = np.arange(len(points)) != j
    loo.append(idw(points[keep], values[keep], points[j:j + 1])[0][0] - values[j])
print(f"leave-one-out RMSE {np.sqrt(np.mean(np.square(loo))):.1f} m; "
      f"true error over every cell {rmse(idw(points, values, grid_points)[0], truth):.1f} m")

pool = np.setdiff1d(np.arange(truth.size), test_idx)
for n in (100, 250, 500, 1000, 2000):
    pick = np.random.default_rng(n).choice(pool, n, replace=False)
    p_, v_ = grid_points[pick], truth[pick]
    tin = griddata(p_, v_, test_points, method="linear")
    print(f"{n:5,} samples: IDW {rmse(idw(p_, v_, test_points)[0], test_values):5.1f} m, "
          f"TIN coverage {np.isfinite(tin).mean():.1%}")

clustered sample: 377 points
      0-   100 m    9,680 cells  RMSE   21.2 m
    100-   250 m   14,020 cells  RMSE   50.3 m
    250-   500 m   15,753 cells  RMSE  101.6 m
    500-  1000 m   25,193 cells  RMSE  142.2 m
   1000-   inf m   46,243 cells  RMSE  258.7 m
  hull covers 54.0%; RMSE inside 123 m, outside 238 m
  377 samples over 99.8 km²
  ideal grid spacing            515 m
  mean nearest-neighbour        125 m  (random would be 257 m)
  Clark-Evans ratio            0.48  (clustered)
kriging, fitted variogram   nugget        0 m²  range  3,762 m  RMSE 41.3 m
kriging, pykrige's own fit  nugget   16,303 m²  range 11,902 m  RMSE 78.6 m
power: {0.5: 75.3, 1: 67.2, 2: 56.1, 3: 52.5, 4: 52.4, 6: 54.7}
k:     {1: 65.4, 3: 53.7, 6: 52.5, 12: 56.1, 24: 62.2, 48: 69.9, 500: 97.1}
leave-one-out RMSE 55.8 m; true error over every cell 56.0 m
  100 samples: IDW 102.3 m, TIN coverage 94.7%
  250 samples: IDW  80.0 m, TIN coverage 94.5%
  500 samples: IDW  55.4 m, TIN coverage 97.4%
1,000 samples: IDW  40.6 m, TIN coverage 99.0%
2,000 samples: IDW  29.1 m, TIN coverage 99.3%

Every interpolator is a weighted average

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: 21.2 m within 100 m of a sample, 258.7 m beyond a kilometre.

A common heuristic is to clip the surface to the convex hull of the samples. On the clustered sample, the hull covered 54.0% of the study area, and the error inside it was 123 m against 238 m outside. That is a real difference, and it is still a blunt one: inside the hull a location can sit in the gap between two clusters, far from any sample, and the error bands run from 21 m to 259 m according to that distance alone.

Distance to the nearest sample is what predicts error at a location, and it is what you should mask on. The hull is a coarse summary of it.

Why more data beats a better method

Doubling the sample from 500 to 1,000 improved IDW from 55.4 m to 40.6 m — 27%. Switching from IDW to kriging on 500 points improved it from 56.1 m to 41.3 m — 26%.

They are comparable in size, and only one of them is free. But the sampling gain compounds: 2,000 points brought IDW to 29.1 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.

Duplicates, anisotropy, rates and barriers

• 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. 145 times more cells changed RMSE by 0.2 m 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 — when leave-one-out can be trusted
• How to choose a cell size for an interpolated surface — the 447 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 with a fitted variogram (41.3 m RMSE) beat linear TIN (49.1 m), IDW (56.1 m) and nearest neighbour (65.4 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 21.2 m within 100 m of a sample to 258.7 m beyond a kilometre in the measurement above.

Should I clip the surface to the convex hull?

The hull is a coarse proxy for support. On a clustered sample it covered 54% of the area, with error inside about half of that outside — while across the map error ranged twelve-fold with distance to the nearest sample. Mask on that distance instead.

What cell size should I interpolate to?

Around the mean sample spacing, sqrt(area/n). Finer grids cost storage and give no accuracy — 145× more cells changed RMSE by 0.2 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 55.4 m to 40.6 m to 29.1 m.