How to Interpolate Points to a Grid with IDW in Python

On this page

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/\text{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     52.5 m RMSE
nearest 12 points    56.1 m
nearest 24 points    62.2 m
all 500 points       97.1 m

Using every sample instead of the nearest six is 85% 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

The power controls how quickly influence decays. The neighbour count controls whether distant samples get a vote at all.

What you get — from The data used on this page

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

Project, small k, any power, a search cap, a matched cell

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 65.4  53.7  52.5  56.1  62.2  69.9  97.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 85%; the power is worth 7% (Example 4).

3. Choose the power, then stop worrying about it

power  0.5     1      2      3      4      6
RMSE  75.3   67.2   56.1   52.5   52.4   54.7

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, on a clustered design where the far bands are large (My interpolated surface extends far beyond the data):

   0 -   100 m    21.2 m
 100 -   250 m    50.3 m
 250 -   500 m   101.6 m
 500 - 1,000 m   142.2 m
over 1,000 m     258.7 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 99.8 km², the mean spacing is 447 m. Output grids:

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

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

The default of "use every sample" is the single most expensive choice in IDW.

IDW to GeoTIFF with k tuned by leave-one-out

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.

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


_ = idw_raster(points, values, bounds, 30, "EPSG:32630", max_distance=3959, out_path="idw.tif")

  333 x 333 at 30 m, 0.0% masked

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


_ = tune_k(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
  best k = 6

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)


rmse = lambda pred: float(np.sqrt(np.mean((pred - test_values) ** 2)))
for angle in (0, 45, 90, 135):
    print(f"  major axis {angle:3d}°, ratio 0.5: RMSE {rmse(idw_anisotropic(points, values, test_points, angle, 0.5)):.1f}")
print(f"  isotropic:                  RMSE {rmse(idw(points, values, test_points)):.1f}")

  major axis   0°, ratio 0.5: RMSE 59.8
  major axis  45°, ratio 0.5: RMSE 61.6
  major axis  90°, ratio 0.5: RMSE 62.3
  major axis 135°, ratio 0.5: RMSE 60.6
  isotropic:                  RMSE 56.1

On this terrain every stretched version was worse than isotropic IDW. Use it only when a directional variogram shows a longer range one way.

Stretching the search along the wrong axis costs accuracy; check a directional variogram first.

Example 4 — the parameter sweeps in the steps

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

k     {1: 65.4, 3: 53.7, 6: 52.5, 12: 56.1, 24: 62.2, 48: 69.9, 500: 97.1}
power {0.5: 75.3, 1: 67.2, 2: 56.1, 3: 52.5, 4: 52.4, 6: 54.7}

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.

Flattening, bullseyes and no uncertainty

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 52.5 m at k=6 to 97.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: 81.8% 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 56.1 m to 52.4 m — 7%. 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 (75.3 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 63% in How to fit a variogram and krige a surface.

Of the four, only power=2 is a default worth keeping.

Zero distances, 1-D arrays and ties

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.

Internal links

Spatial interpolation explained — how error grows with distance
My IDW surface is full of bullseyes — the artefact this method always produces
IDW, kriging, splines or TIN? Choosing an interpolator — when to use something else
How to fit a variogram and krige a surface in Python — the accurate alternative
How to choose a cell size for an interpolated surface — sizing the output grid
How to cross-validate an interpolated surface — tuning k honestly
The variogram explained — where max_distance comes from
My interpolated surface extends far beyond the data — why to cap the search radius

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/\text{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 52.5 m to 97.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 — 81.8% 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 — 21.2 m within 100 m of a sample, 258.7 m beyond a kilometre, on a clustered design.

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 26% more accurate on these samples.