⚙️ How-To Interpolation & surfaces Raster & rasterio Interpolation & Surfaces · step 3 of 13 Updated

How to Interpolate Points to a Grid with IDW in Python

Problem statement

Inverse distance weighting is the interpolator you can write from memory: predict each location as a weighted average of nearby samples, with weights falling off as 1/distance^p.

It is also the interpolator with the most ways to be quietly wrong. The defaults in most implementations — all points, power 2, no distance limit — cost real accuracy. Measured on 500 samples from a real elevation model:

nearest 6 points     51.5 m RMSE
nearest 12 points    54.0 m
nearest 24 points    59.6 m
all 500 points       90.1 m

Using every sample instead of the nearest six is 75% worse. That is the single largest parameter effect in the method, and it is the one most implementations leave switched on.

Quick answer 

import numpy as np
from scipy.spatial import cKDTree


def idw(points, values, targets, power=2, k=12, max_distance=None):
    """Inverse distance weighting over the k nearest samples."""
    points = np.asarray(points, dtype=float)
    values = np.asarray(values, dtype=float)
    targets = np.asarray(targets, dtype=float)

    tree = cKDTree(points)
    distance, index = tree.query(targets, k=min(k, len(points)))
    # query returns 1-D arrays when k == 1
    distance = distance.reshape(len(targets), -1)
    index = index.reshape(len(targets), -1)

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

    # a target sitting exactly on a sample takes that sample's value
    exact = distance[:, 0] < 1e-9
    predicted[exact] = values[index[exact, 0]]

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

    return predicted
Inverse distance weights falling off with distance for powers of one, two and four, showing how a higher power concentrates weight on the nearest sample.
The power controls how quickly influence decays. The neighbour count controls whether distant samples get a vote at all.

Step-by-step solution

1. Project the coordinates first

IDW weights by distance, and distance in degrees is not distance. One degree of longitude is 111 km at the equator and 71 km at 50° north; weighting in degrees stretches the east-west axis by whatever 1/cos(latitude) happens to be.

gdf = gdf.to_crs(gdf.estimate_utm_crs())
points = np.column_stack([gdf.geometry.x, gdf.geometry.y])

2. Set k, and set it small 

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

Six to twelve is the sweet spot. Below that the surface is noisy and tends towards nearest-neighbour blockiness; above it, distant samples with small but non-zero weights collectively drag every prediction towards the global mean, flattening the map.

This is the parameter to tune. It is worth 75%; the power is worth 6%.

3. Choose the power, then stop worrying about it 

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

Anything between 2 and 4 is fine. Below 2 the surface is over-smoothed towards the local mean; above 4 it collapses towards nearest neighbour, with each prediction dominated by its single closest sample.

4. Cap the search distance

Beyond the variogram range, distance carries no information about similarity — so weighting by it produces a confident-looking number with nothing behind it. Measured error by distance to the nearest sample:

   0 -   100 m    18.9 m
 100 -   250 m    54.0 m
 250 -   500 m   106.1 m
 500 - 1,000 m   169.8 m
over 1,000 m     334.3 m

Set max_distance to something defensible — the variogram range, or two to three times the mean sample spacing — and let the far areas come out as NaN.

5. Choose an output cell size that matches the sample spacing

With 500 samples over 105 km², the mean spacing is 457 m. Output grids:

cell size    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

Ninety times more cells buys nothing. See How to choose a cell size for an interpolated surface.

IDW error against neighbour count rising from 51.5 metres at six neighbours to 90.1 metres using all 500 samples.
The default of "use every sample" is the single most expensive choice in IDW.

Code examples

Example 1 — IDW to a raster, written straight to GeoTIFF 

import numpy as np
import rasterio
from rasterio.transform import from_origin
from scipy.spatial import cKDTree


def idw_raster(points, values, bounds, cell_size, crs, power=2, k=12,
               max_distance=None, out_path=None):
    """Interpolate to a regular grid and write it with its support raster."""
    left, bottom, right, top = bounds
    width = int(np.ceil((right - left) / cell_size))
    height = int(np.ceil((top - bottom) / cell_size))
    transform = from_origin(left, top, cell_size, cell_size)

    xs = left + (np.arange(width) + 0.5) * cell_size
    ys = top - (np.arange(height) + 0.5) * cell_size
    grid = np.column_stack([np.repeat(xs[None, :], height, 0).ravel(),
                            np.repeat(ys[:, None], width, 1).ravel()])

    tree = cKDTree(points)
    distance, index = tree.query(grid, k=min(k, len(points)))
    distance = distance.reshape(len(grid), -1)
    index = index.reshape(len(grid), -1)

    weights = 1.0 / np.maximum(distance, 1e-12) ** power
    predicted = (weights * np.asarray(values)[index]).sum(1) / weights.sum(1)
    nearest = distance[:, 0]
    if max_distance is not None:
        predicted = np.where(nearest <= max_distance, predicted, np.nan)

    surface = predicted.reshape(height, width).astype("float32")
    support = nearest.reshape(height, width).astype("float32")

    if out_path:
        profile = dict(driver="GTiff", height=height, width=width, count=2,
                       dtype="float32", crs=crs, transform=transform,
                       nodata=np.nan, compress="deflate", tiled=True)
        with rasterio.open(out_path, "w", **profile) as dst:
            dst.write(surface, 1)
            dst.write(support, 2)
            dst.set_band_description(1, "interpolated value")
            dst.set_band_description(2, "distance to nearest sample (m)")
            dst.update_tags(method="idw", power=power, k=k,
                            max_distance=str(max_distance))

    print(f"  {width} x {height} at {cell_size} m, "
          f"{np.isnan(surface).mean():.1%} masked")
    return surface, support, transform

Writing the support distance as a second band is the habit worth forming. It travels with the surface, it costs one extra band, and it is the only thing that tells the next person which parts of the map to believe.

Example 2 — tuning k by leave-one-out 

import numpy as np
from scipy.spatial import cKDTree


def tune_k(points, values, candidates=(3, 6, 9, 12, 18, 24), power=2):
    """Leave-one-out RMSE for each neighbour count."""
    tree = cKDTree(points)
    results = {}
    for k in candidates:
        distance, index = tree.query(points, k=min(k + 1, len(points)))
        distance, index = distance[:, 1:], index[:, 1:]     # drop self
        weights = 1.0 / np.maximum(distance, 1e-12) ** power
        predicted = (weights * values[index]).sum(1) / weights.sum(1)
        results[k] = float(np.sqrt(np.mean((predicted - values) ** 2)))
        print(f"  k={k:3d}  LOO RMSE {results[k]:7.2f}")
    best = min(results, key=results.get)
    print(f"  best k = {best}")
    return best

The [:, 1:] slice removes each point from its own prediction. Without it, every sample predicts itself exactly, LOO RMSE is zero for every k, and the tuning is meaningless.

Note that this tunes against the sample locations, not against the whole surface, so it inherits whatever bias the sampling design has — see How to cross-validate an interpolated surface.

Example 3 — anisotropic IDW, when structure has a direction 

import numpy as np
from scipy.spatial import cKDTree


def idw_anisotropic(points, values, targets, major_angle_deg, ratio,
                    power=2, k=12):
    """IDW in a stretched coordinate system, so influence reaches further
    along the major axis than across it.

    ratio: minor/major, e.g. 0.25 means influence reaches 4x further along.
    """
    theta = np.radians(major_angle_deg)
    rotate = np.array([[np.cos(theta), np.sin(theta)],
                       [-np.sin(theta), np.cos(theta)]])
    stretch = np.array([[1.0, 0.0], [0.0, 1.0 / ratio]])
    transform = stretch @ rotate

    p = points @ transform.T
    t = targets @ transform.T

    tree = cKDTree(p)
    distance, index = tree.query(t, k=min(k, len(p)))
    distance = distance.reshape(len(t), -1); index = index.reshape(len(t), -1)
    weights = 1.0 / np.maximum(distance, 1e-12) ** power
    return (weights * values[index]).sum(1) / weights.sum(1)

Rainfall along a valley, a contaminant plume down a hydraulic gradient and coastal salinity all have a preferred direction. Get the angle and ratio from directional variograms rather than by eye.

Explanation

Why too many neighbours flattens the map

With power 2, a sample at 5 km has a weight of 1/25,000,000 — tiny, but not zero. With one such sample it is irrelevant. With four hundred of them spread around the study area, their weights sum to something comparable with the handful of near samples, and they all pull towards the same place: the global mean.

That is exactly what the measurement shows. RMSE climbed from 51.5 m at k=6 to 90.1 m using all 500 samples, and the surface flattened towards the average elevation.

Restricting to the nearest k is not an approximation for speed. It is the correct model: only nearby samples carry information.

Why IDW produces bullseyes

IDW is an exact interpolator — at a sample location the weight on that sample is infinite relative to all others, so the surface passes exactly through every observation.

The consequence is that every sample is a local extremum of the surface. Measured on a 500-point IDW surface: 83.1% of sample locations were local maxima or minima of the interpolated grid, against 0.38% of cells being local extrema in the real DEM.

Each of those becomes a concentric ring in the output. See My IDW surface is full of bullseyes.

Why the power matters so little

Between power 2 and power 4 the RMSE moved from 54.0 m to 50.9 m — under 6%. The reason is that within a small neighbourhood, all reasonable weightings put most of the weight on the same two or three nearest samples.

The power only becomes decisive at the extremes: at 0.5 the weights are nearly equal and the surface is a local mean (72.2 m), and above 6 it is nearly nearest-neighbour.

Why IDW cannot give you uncertainty

There is no statistical model behind IDW. The weights come from an exponent you chose, not from measured spatial structure, so there is no expression for the variance of the prediction.

What you can and should report instead is the distance to the nearest sample, which is the geometric part of uncertainty. If you need a real variance, use kriging — while remembering that kriging variance overstated the true error by 68% in How to fit a variogram and krige a surface.

Four IDW settings to change before trusting a surface: neighbour count, maximum distance, power and CRS.
Of the four, only `power=2` is a default worth keeping.

Edge cases or notes

  • Project first. Degrees are not metres, and IDW weights by distance.
  • Guard the zero distance. A target exactly on a sample gives 1/0; take that sample's value.
  • cKDTree.query returns 1-D arrays when k=1. Reshape, or the weighting silently transposes.
  • Set k to 6–12 and cap max_distance.
  • IDW never extrapolates beyond the sample range. Predictions are bounded by the minimum and maximum observed values, so a true peak between samples is always underestimated.
  • Duplicate coordinates are harmless for IDW (unlike kriging) but double-count that location.
  • Use a projected equal-area CRS for large areas; UTM distorts distance beyond a few zones.
  • A tie at the kth neighbour is resolved arbitrarily by the tree, which makes the surface very slightly non-deterministic across scipy versions.

FAQ

How do I do IDW in Python?

Build a scipy.spatial.cKDTree over the sample points, query the k nearest for each target, weight them by 1/distance^p and take the weighted mean. Guard the zero-distance case.

What power should I use for IDW?

Between 2 and 4. The difference across that range was under 6% here, so it is not worth extensive tuning.

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 drag every prediction to the global mean.

Why does my IDW surface look like a set of rings?

Because IDW is an exact interpolator, so every sample is a local extremum — 83.1% of sample locations were extrema in the surface measured here.

Can IDW extrapolate?

It will produce values anywhere, but they are always bounded by the observed minimum and maximum, and the error grows steeply with distance — 18.9 m within 100 m of a sample, 334.3 m beyond a kilometre.

Does IDW give me an uncertainty estimate?

No. There is no statistical model behind it. Report the distance to the nearest sample instead, or use kriging.

Should I use IDW or kriging?

IDW when you need something defensible quickly and the samples are dense and even. Kriging when the variogram shows clear structure, the field is anisotropic, or you need an uncertainty map. Kriging was 25% more accurate here.