My IDW Surface Is Full of Bullseyes

On this page

Problem statement

An inverse-distance-weighted surface comes out covered in concentric rings, one around each sample point, with a flat plain between them. It looks like a bad radar image rather than a landscape.

This is not a bug and it is not a parameter you set wrongly. It is a structural property of the method: IDW is an exact interpolator, so every sample is a local extremum of the surface.

Measured on a 500-point IDW surface derived from a real elevation model:

sample points that are a local extremum of the IDW surface   81.8%
cells that are a local extremum in the real DEM               0.38%

Four out of five samples become their own little peak or pit. Those are the bullseyes.

Quick answer

You cannot remove bullseyes from IDW. You can make them smaller, or use a method that does not have them:

# 1. lower the power — flattens the spike at each sample
idw(points, values, targets, power=1.5, k=12)

# 2. use more neighbours — but not many more
idw(points, values, targets, power=2, k=16)

# 3. use a smoothing distance so no weight is ever infinite
def idw_smoothed(points, values, targets, power=2, k=12, smoothing=50.0):
    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 / (d ** 2 + smoothing ** 2) ** (power / 2)   # never infinite
    return (w * values[i]).sum(1) / w.sum(1)

The third is the real fix. Adding a smoothing distance makes IDW an inexact interpolator: it no longer passes exactly through the samples, so the samples stop being extrema.

The weight at a sample is infinite relative to all others, so the surface must pass exactly through it — and then fall away in every direction.

Confirm, then lower the power, smooth or switch

1. Confirm it is bullseyes and not something else

import numpy as np


def local_extrema(grid):
    """Cells higher or lower than all eight neighbours."""
    centre = grid[1:-1, 1:-1]
    stack = np.stack([grid[:-2, 1:-1], grid[2:, 1:-1],
                      grid[1:-1, :-2], grid[1:-1, 2:],
                      grid[:-2, :-2], grid[:-2, 2:],
                      grid[2:, :-2], grid[2:, 2:]])
    return (centre > stack.max(0)) | (centre < stack.min(0))


for name, grid in (("IDW surface", surface), ("real DEM   ", dem)):
    extrema = local_extrema(grid)
    print(f"{name}: {extrema.sum():,} local extrema ({extrema.mean():.2%} of cells)")

IDW surface: 670 local extrema (0.61% of cells)
real DEM   : 418 local extrema (0.38% of cells)

Then check whether those extrema sit on your sample points. If most of them do, it is bullseyes. If they are scattered, it is noise in the data.

2. Lower the power

A higher power concentrates weight on the nearest sample, which makes the spike at each sample sharper and the plain between them flatter. Power 4 gives dramatic bullseyes; power 1.5 gives soft ones.

The cost is small: measured RMSE went from 52.5 m at power 3 to 56.1 m at power 2 to 67.2 m at power 1 (in Spatial interpolation explained, on the same samples). Dropping from 2 to 1.5 is a modest accuracy price for a visibly better surface.

3. Add a smoothing distance

The mathematical cause is $1/d^p$ going to infinity as $d$ goes to zero. Replace it with $1/(d^2 + s^2)^{p/2}$ and the weight at zero distance is finite, so no sample dominates absolutely.

Choose s as a fraction of the mean sample spacing — a tenth is a reasonable start. With 447 m spacing that is about 45 m; measured, s = 50 m halved the share of samples that are extrema (81.8% to 49.1%) for 0.3 m of RMSE, and s = 100 m brought it to 13.1% for 1.4 m (Example 2).

This is the standard trade: exactness for smoothness. If your samples have measurement error, you want an inexact interpolator, because honouring a noisy measurement exactly is honouring the noise.

4. Or switch method

Bullseyes are specific to distance-weighted exact interpolators. The alternatives behave differently:

Kriging with a nugget is inexact by construction: a non-zero nugget means the surface does not pass through the data, and there are no rings. With a zero nugget it is exact and does have them, though far weaker than IDW's because the weights come from structure.
Linear TIN is exact but produces flat triangular facets, not rings. Each sample is a vertex, so it is still a local extremum — but the artefact reads as faceting.
Splines are exact and smooth, and swap rings for overshoot near steep changes.

5. Consider that the bullseyes are telling you something

Prominent rings mean the samples are far apart relative to the variation between them. The surface is honestly reporting that it knows the value at 500 places and is guessing everywhere else.

Smoothing them away makes the map prettier without adding information. If the rings are severe, the fix that actually helps is more samples: doubling from 500 to 1,000 cut IDW's RMSE from 55.4 m to 40.6 m and roughly halves the visual prominence of each ring.

Exactness is the cause. A smoothing distance buys a surface that passes near the data instead of through it.

Measuring and smoothing the bullseyes

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.

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 — measuring how bad the bullseyes are

import numpy as np
from scipy.spatial import cKDTree


def bullseye_score(surface, transform, points):
    """What fraction of sample locations are local extrema of the surface?"""
    a, _, c, _, e, f = transform[:6]
    cols = ((points[:, 0] - c) / a).astype(int)
    rows = ((points[:, 1] - f) / e).astype(int)

    height, width = surface.shape
    inside = (rows > 0) & (rows < height - 1) & (cols > 0) & (cols < width - 1)

    extrema = local_extrema(surface)
    hits = extrema[rows[inside] - 1, cols[inside] - 1]

    print(f"  {int(inside.sum())} samples inside the grid")
    print(f"  {hits.mean():.1%} of them are local extrema of the surface")
    print(f"  {extrema.mean():.2%} of all cells are local extrema")
    return float(hits.mean())


_ = bullseye_score(surface, transform, points)

  495 samples inside the grid
  81.8% of them are local extrema of the surface
  0.61% of all cells are local extrema

A score near 100% means a pure exact interpolator. Below about 20% the artefact is no longer visually obvious.

Example 2 — comparing smoothing distances

import numpy as np
from scipy.spatial import cKDTree


def smoothed_idw(points, values, targets, power=2, k=12, smoothing=0.0):
    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)

    if smoothing > 0:
        w = 1.0 / (d ** 2 + smoothing ** 2) ** (power / 2)
    else:
        w = 1.0 / np.maximum(d, 1e-12) ** power
        exact = d[:, 0] < 1e-9
    out = (w * values[i]).sum(1) / w.sum(1)
    if smoothing <= 0:
        out[exact] = values[i[exact, 0]]
    return out


def smoothing_sweep(points, values, grid_points, truth, shape, transform,
                    candidates=(0, 25, 50, 100, 200)):
    for s in candidates:
        surface = smoothed_idw(points, values, grid_points,
                               smoothing=s).reshape(shape)
        rmse = float(np.sqrt(np.mean((surface.ravel() - truth) ** 2)))
        score = bullseye_score(surface, transform, points)
        print(f"  smoothing {s:4d} m: RMSE {rmse:6.2f}  bullseye score {score:5.1%}")


smoothing_sweep(points, values, grid_points, truth, shape, transform)

  495 samples inside the grid
  81.8% of them are local extrema of the surface
  0.61% of all cells are local extrema
  smoothing    0 m: RMSE  55.99  bullseye score 81.8%
  495 samples inside the grid
  71.1% of them are local extrema of the surface
  0.56% of all cells are local extrema
  smoothing   25 m: RMSE  56.07  bullseye score 71.1%
  495 samples inside the grid
  49.1% of them are local extrema of the surface
  0.51% of all cells are local extrema
  smoothing   50 m: RMSE  56.33  bullseye score 49.1%
  495 samples inside the grid
  13.1% of them are local extrema of the surface
  0.45% of all cells are local extrema
  smoothing  100 m: RMSE  57.39  bullseye score 13.1%
  495 samples inside the grid
  1.0% of them are local extrema of the surface
  0.43% of all cells are local extrema
  smoothing  200 m: RMSE  60.84  bullseye score  1.0%

Run this on your own data and pick the smallest smoothing that gets the score below about 20%. Beyond that you are trading accuracy for cosmetics.

A tenth of the sample spacing halves the rings; about a fifth removes them.

Example 3 — the honest alternative: show the support

import numpy as np
from scipy.spatial import cKDTree


def surface_with_support(points, values, targets, shape, power=2, k=12,
                         range_m=None):
    """Interpolate, and mask where no sample is within the variogram range."""
    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-12) ** power
    surface = ((w * values[i]).sum(1) / w.sum(1)).reshape(shape)

    nearest = d[:, 0].reshape(shape)
    if range_m:
        surface = np.where(nearest <= range_m, surface, np.nan)
        print(f"  masked {np.isnan(surface).mean():.1%} beyond {range_m} m")
    return surface, nearest


_ = surface_with_support(points, values, grid_points, shape, range_m=3762)

  masked 0.0% beyond 3762 m

Rings are a symptom of sparse sampling, and sparse sampling is what the support raster shows directly. Publishing both is more useful than smoothing the surface until the sparseness is invisible.

Exactness creates the rings

Why exactness causes extrema

At a sample location, $d = 0$ for that sample and positive for all others. The weight $1/0^p$ is infinite, so the weighted average equals that sample's value exactly.

Move a few metres away and the sample's weight becomes large but finite, and the other samples begin to contribute. Their values pull the estimate towards a local average. Whether the sample is above or below that local average decides whether the surface falls or rises away from it — and either way, the sample is a local extremum.

Since every sample is above or below its local average (a sample exactly equal to it is measure zero), essentially every sample becomes a local extremum. The measurement of 83.1% is the practical version of "essentially every"; the shortfall is samples whose immediate neighbourhood is a plateau or whose closest neighbour is within a cell.

Why the rings are circular

The weight depends only on distance, so the surface around an isolated sample has circular symmetry. Contour a spike with circular symmetry and you get concentric circles — the bullseye.

That is also why the artefact is worse when samples are isolated: near several samples, the influences overlap and the symmetry breaks. Bullseyes are most visible exactly where they matter most, in the sparsely sampled parts of the map.

Why a nugget removes them in kriging

A kriging system with a non-zero nugget has a discontinuity at zero distance: γ(0) = 0 but γ(h) → nugget as h → 0⁺. The consequence is that the surface does not pass through the data; it passes near it, with the gap controlled by the nugget-to-sill ratio.

That is the same trick as IDW's smoothing distance, arrived at from a statistical model rather than by patching the weight function — and with the advantage that the amount of smoothing is estimated from the data rather than chosen.

Why smoothing is not always the right answer

If your measurements are accurate, honouring them exactly is correct and the bullseyes are an honest depiction of a sparse survey. Smoothing then hides real information: the surface no longer tells you where the observations were.

If your measurements have error — most field data — then exact interpolation reproduces that error as topography, and smoothing is a genuine improvement rather than cosmetics. The nugget is the formal way to say how much error you believe there is.

Post-smoothing is the popular one. It changes the picture and not the numbers.

Power, neighbour count and contour intervals

Bullseyes get worse as the power rises. Power 4 concentrates almost all weight on the nearest sample.
They get worse with fewer neighbours too. k=1 is pure nearest neighbour, which is all plateau and no ring.
A smoothing distance makes IDW inexact. That is the point, and it should be recorded in the output metadata.
Kriging with a zero nugget is also exact and also produces rings, weaker than IDW's.
Linear TIN is exact but shows faceting, not rings.
Duplicate points at one location produce a doubled spike; aggregate first.
Rings around a single sample with none elsewhere usually means that sample is an outlier — check it before smoothing it away.
Contour intervals amplify the artefact. The same surface with fewer contours can look fine.

Internal links

How to interpolate points to a grid with IDW in Python — the parameters that control severity
Spatial interpolation explained — why sparse sampling is the underlying cause
IDW, kriging, splines or TIN? Choosing an interpolator — methods without rings
How to fit a variogram and krige a surface in Python — the nugget as principled smoothing
The variogram explained — what the nugget actually measures
Sample design explained: where to measure — the fix that adds information
My interpolated surface extends far beyond the data — the other sparse-sampling artefact
How to find and remove spatial outliers in a point dataset — when one ring is a bad measurement

FAQ

Why does my IDW surface have circular rings around each point?

Because IDW is an exact interpolator: the weight at a sample is infinite relative to the others, so the surface must pass through it and then fall away in every direction. 81.8% of samples became local extrema in the surface measured here.

How do I get rid of IDW bullseyes?

Add a smoothing distance so the weight at zero distance is finite: $1/(d^2 + s^2)^{p/2}$ . Start with about a tenth of the mean sample spacing.

Does lowering the power help?

Yes, somewhat. It flattens the spike at each sample, at a modest accuracy cost — RMSE went from 52.5 m at power 3 to 67.2 m at power 1 here.

Does kriging have bullseyes?

Only with a zero nugget, and much weaker. A non-zero nugget makes kriging inexact, and the rings disappear.

Are bullseyes always bad?

No. With accurate measurements they honestly show where the observations are. They are a problem when the data have measurement error, because then the surface is reproducing noise as topography.

What if only one point has a ring around it?

That point is probably an outlier — a transcription error or a genuinely unusual measurement. Investigate it rather than smoothing the whole surface.

Will more samples fix it?

They will make each ring smaller and less prominent, and they improve accuracy for real: 500 to 1,000 samples cut RMSE from 55.4 m to 40.6 m here.