Sample Design Explained: Where to Measure

Problem statement

Sampling design decides the accuracy of an interpolated surface more than the interpolation method does, and it is fixed before any code is written.

Three designs, the same 500 samples, the same IDW interpolation, the same real elevation surface as ground truth:

design       RMSE     mean gap   95th pct gap   worst gap
random      49.2 m      227 m         450 m        819 m
regular     44.8 m      174 m         272 m        671 m
clustered  257.1 m    2,365 m       6,244 m      8,385 m

A regular grid is 9% more accurate than random placement at the same cost, and a stratified-random design sits between them at 47.5 m (Example 4). Clustered sampling — which is what a survey following roads, rivers or convenient access produces — is more than five times worse.

Quick answer

Two numbers tell you whether a design will support a surface, and both are computable before you interpolate anything:

import numpy as np
from scipy.spatial import cKDTree


def design_quality(points, bounds):
    """Coverage, not count, is what determines interpolation error."""
    left, bottom, right, top = bounds
    area = (right - left) * (top - bottom)

    ideal = np.sqrt(area / len(points))       # spacing of a perfect grid

    # how far is a typical location from its nearest sample?
    xs = np.random.default_rng(0).uniform(left, right, 20000)
    ys = np.random.default_rng(1).uniform(bottom, top, 20000)
    gap, _ = cKDTree(points).query(np.column_stack([xs, ys]), k=1)

    print(f"  {len(points)} samples over {area / 1e6:.1f} km²")
    print(f"  ideal grid spacing   {ideal:8.0f} m")
    print(f"  mean gap             {gap.mean():8.0f} m")
    print(f"  95th percentile gap  {np.percentile(gap, 95):8.0f} m")
    print(f"  worst gap            {gap.max():8.0f} m")
    return {"ideal": float(ideal), "mean_gap": float(gap.mean()),
            "p95_gap": float(np.percentile(gap, 95)),
            "max_gap": float(gap.max())}

The 95th percentile gap is the number to design against. It is what determines the error in the worst parts of the map, which is where every complaint comes from.

Design against the gap: grid, stratified, pairs, boundary

1. Design against the gap, not the count

"500 samples" says nothing. 500 samples in five clusters left a 95th-percentile gap of 6,244 m and produced an RMSE of 257.1 m — worse than 100 random samples, which scored 102.3 m in Spatial interpolation explained.

Coverage is the quantity that matters, and coverage is a distribution, not a number. Report the mean, the 95th percentile and the maximum gap.

2. Prefer a regular grid, with one caveat

The grid was the most accurate design here: 44.8 m against random's 49.2 m, because it has no gaps larger than its spacing. On this rectangular window it also had the smallest worst gap, 671 m against random's 819 m. On an irregular study area that changes: a grid laid over a ragged boundary leaves strips along the edges that no sample reaches.

That is the standard trade-off. Grids minimise typical gaps and can leave systematic edge gaps. Check the boundary explicitly.

Grids have one real failure mode: if the phenomenon is periodic at the grid spacing — plough furrows, tree rows, street blocks — every sample lands on the same phase and the survey is systematically biased. Offsetting alternate rows or using a stratified-random design inside grid cells removes this.

3. Use stratified random when you need both

Divide the area into equal cells and place one random point in each. You get the coverage guarantee of a grid without the periodicity risk, and the random offsets provide the short-distance pairs a variogram needs.

4. Add short-distance pairs on purpose

A variogram needs pairs at every distance, including short ones. A pure grid at 400 m spacing has no pairs below 400 m, so the nugget is unmeasurable and the near-origin part of the curve — the part that matters most for interpolation — is pure extrapolation.

Add a handful of nested clusters: a few sites with additional samples at 10 m, 30 m and 100 m. It costs a small fraction of the budget and is the only way to measure the nugget.

5. Sample the boundary, not just the interior

Interpolation degrades at edges because targets there have samples on one side only. Measured across distance bands on a clustered design in cross-validating an interpolation, error rose from 21.2 m within 100 m of a sample to 258.7 m beyond a kilometre — and the far band is disproportionately edges and gaps.

If the study area has a boundary you care about, put samples on it.

Three designs and where the next sample goes

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, used to interpolate each design in Example 4, where every design's samples read the DEM through truth_fn.

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 — generating the three designs

import numpy as np


def sample_design(bounds, n, kind="grid", seed=0, clusters=5, spread=400):
    """Random, regular grid or clustered sample points within bounds."""
    rng = np.random.default_rng(seed)
    left, bottom, right, top = bounds

    if kind == "random":
        return np.column_stack([rng.uniform(left, right, n),
                                rng.uniform(bottom, top, n)])

    if kind == "grid":
        aspect = (right - left) / (top - bottom)
        cols = int(round(np.sqrt(n * aspect)))
        rows = int(np.ceil(n / cols))
        step_x = (right - left) / cols
        step_y = (top - bottom) / rows
        xs = left + (np.arange(cols) + 0.5) * step_x
        ys = bottom + (np.arange(rows) + 0.5) * step_y
        grid = np.array([(x, y) for y in ys for x in xs])
        return grid[:n]

    if kind == "stratified":
        aspect = (right - left) / (top - bottom)
        cols = int(round(np.sqrt(n * aspect)))
        rows = int(np.ceil(n / cols))
        step_x = (right - left) / cols
        step_y = (top - bottom) / rows
        pts = [(left + (i + rng.random()) * step_x,
                bottom + (j + rng.random()) * step_y)
               for j in range(rows) for i in range(cols)]
        return np.array(pts)[:n]

    if kind == "clustered":
        centres = np.column_stack([rng.uniform(left, right, clusters),
                                   rng.uniform(bottom, top, clusters)])
        which = rng.integers(0, clusters, n)
        pts = centres[which] + rng.normal(0, spread, (n, 2))
        return np.clip(pts, [left, bottom], [right, top])

    raise ValueError(f"unknown design {kind!r}")

Example 2 — nested clusters for the variogram

import numpy as np


def add_nested_pairs(points, n_sites=6, offsets=(10, 30, 100, 300), seed=0):
    """Extra samples at short lags, so the near-origin variogram is measured."""
    rng = np.random.default_rng(seed)
    chosen = rng.choice(len(points), size=min(n_sites, len(points)),
                        replace=False)

    extra = []
    for idx in chosen:
        for distance in offsets:
            angle = rng.uniform(0, 2 * np.pi)
            extra.append(points[idx] + distance * np.array([np.cos(angle),
                                                            np.sin(angle)]))
    extra = np.array(extra)
    print(f"  {len(extra)} extra samples at {offsets} m from {n_sites} sites")
    print(f"  cost: {len(extra) / len(points):.0%} of the original budget")
    return np.vstack([points, extra])


_ = add_nested_pairs(sample_design(bounds, 500, "random"))

  24 extra samples at (10, 30, 100, 300) m from 6 sites
  cost: 5% of the original budget

Five percent of the budget buys the entire short-lag half of the variogram. Without it the nugget is a free parameter the optimiser will choose to fit noise.

Example 3 — where to put the next sample

import numpy as np
from scipy.spatial import cKDTree


def next_sample_locations(existing, bounds, n_new=10, candidates=40000, seed=0):
    """Greedily place new samples where the coverage gap is largest."""
    rng = np.random.default_rng(seed)
    left, bottom, right, top = bounds
    grid = np.column_stack([rng.uniform(left, right, candidates),
                            rng.uniform(bottom, top, candidates)])

    points = np.array(existing, dtype=float)
    chosen = []
    for _ in range(n_new):
        gap, _ = cKDTree(points).query(grid, k=1)
        pick = int(np.argmax(gap))
        chosen.append(grid[pick])
        print(f"  new sample at ({grid[pick][0]:.0f}, {grid[pick][1]:.0f}), "
              f"filling a {gap[pick]:.0f} m gap")
        points = np.vstack([points, grid[pick]])

    return np.array(chosen)


_ = next_sample_locations(sample_design(bounds, 500, "clustered"), bounds, n_new=3)

  new sample at (432859, 5875449), filling a 8403 m gap
  new sample at (427408, 5875522), filling a 5448 m gap
  new sample at (432834, 5879932), filling a 4482 m gap

This is the practical version of "design against the gap": for a second field season, the highest-value locations are exactly the ones furthest from anything already measured. Greedy maximin placement is simple, fast, and hard to beat.

Example 4 — the design comparison at the top of the page

print("design       RMSE     mean gap   95th pct gap   worst gap")
for kind in ("random", "grid", "stratified", "clustered"):
    pts = sample_design(bounds, 500, kind)
    vals = truth_fn(pts)
    rmse = float(np.sqrt(np.mean((idw(pts, vals, grid_points)[0] - truth) ** 2)))
    import contextlib, io
    with contextlib.redirect_stdout(io.StringIO()):
        q = design_quality(pts, bounds)
    print(f"{kind:10} {rmse:6.1f} m {q['mean_gap']:8,.0f} m {q['p95_gap']:10,.0f} m {q['max_gap']:10,.0f} m")

design       RMSE     mean gap   95th pct gap   worst gap
random       49.2 m      227 m        450 m        819 m
grid         44.8 m      174 m        272 m        671 m
stratified   47.5 m      196 m        347 m        789 m
clustered   257.1 m    2,365 m      6,244 m      8,385 m

Why clustered sampling fails twice

Why clustered sampling fails so badly

The clustered design had 500 samples and a mean gap of 2,365 m, against random's 227 m. Every one of those 500 measurements is real, and most of them are nearly redundant: within a cluster, samples are close enough to be strongly correlated, so the tenth sample in a cluster adds a fraction of the information the first one did.

The effective sample size is therefore much smaller than 500 — closer to the number of clusters than to the number of points. Five clusters is five pieces of information about a landscape 10 km across.

This is the design that real surveys drift into. Access is easier along roads, rivers and existing plots, and each individually reasonable decision moves the design towards clustering.

Why clustering also breaks your validation

Clustered sampling is doubly dangerous because it corrupts the check that would reveal it. Leave-one-out cross-validation removes one point and predicts it from the rest — but in a cluster, the removed point still has near neighbours a few metres away, so it is predicted almost perfectly.

The result is an optimistic cross-validation score on a surface that is genuinely bad. On a well-spread random design of the same surface, LOO matched the true error within half a percent (55.8 m against 56.0 m); under clustering it becomes optimistic, and the error grows with the clustering. See How to cross-validate an interpolated surface.

Why more samples in the same place do not help

Doubling a well-spread sample from 500 to 1,000 improved IDW's RMSE from 55.4 m to 40.6 m. Doubling a clustered sample mostly adds more points inside the existing clusters, filling in distances where the surface was already accurate.

The value of a sample is roughly proportional to how much new area it brings within the variogram range. Once a location is within range of an existing sample, another one there is nearly free of information.

Why the range sets the spacing

The variogram range is the distance beyond which samples tell you nothing about each other. Sample spacing well above the range means every point is isolated and interpolation is impossible. Spacing well below it means neighbouring samples are nearly redundant.

A reasonable target is a spacing of roughly a third to a half of the range, which keeps every location within range of several samples without paying for redundancy. With the fitted range of 3,959 m in The variogram explained, that suggests 1.3–2 km spacing — and the 500-point designs above, at 447 m ideal spacing, are comfortably denser than needed for that structure. Their remaining error is dominated by variation at scales shorter than the range.

Gap percentiles, periodic features and declustering

• Report gap percentiles, not sample counts. 500 clustered samples had a 95th-percentile gap of 3.6 km.
• Grids can miss the boundary. Check the maximum gap; the grid here had the worst single gap of the three designs.
• Avoid a grid spacing that matches a periodic feature — furrows, street blocks, tree rows.
• Stratified random gets most of the grid's coverage without the periodicity risk.
• Budget about 5% for nested short-lag pairs, or the nugget is unmeasurable.
• Sample the boundary if you care about the edges of the map.
• Preferential sampling biases the mean too, not just the surface: if high values are easier to reach, the map is high everywhere.
• Declustering weights can partly rescue an existing clustered dataset, but they cannot create coverage that was never measured.

Internal links

• Spatial interpolation explained — how error grows with the gap
• The variogram explained — the range that sets the spacing
• How to cross-validate an interpolated surface — why clustering corrupts validation
• My interpolated surface extends far beyond the data — what a large gap looks like in the output
• IDW, kriging, splines or TIN? Choosing an interpolator — the smaller lever
• How to choose a cell size for an interpolated surface — spacing on the output side
• The modifiable areal unit problem explained — the same "your units decide your answer" trap
• How to find and remove spatial outliers in a point dataset — cleaning before designing round the data you have
• Sampling bias in occurrence data explained — what opportunistic data looks like against a designed sample

FAQ

How many samples do I need?

Enough that the 95th-percentile distance from any location to its nearest sample is comfortably inside the variogram range. Count alone is meaningless — 500 clustered samples left a 6.2 km gap at the 95th percentile.

Is a regular grid better than random sampling?

For coverage, yes: 44.8 m RMSE against 49.2 m for the same 500 points here. But grids can leave systematic edge gaps, and they are vulnerable to periodic features in the landscape.

What is wrong with clustered sampling?

It wastes budget on redundant nearby measurements and leaves large unsampled gaps. The same 500 points gave 257.1 m RMSE against 49.2 m for random placement — and cross-validation will not warn you.

What is stratified random sampling?

One random point inside each cell of a regular grid. It gets the grid's coverage guarantee without the risk of aligning with a periodic feature, and provides varied short-distance pairs.

Why do I need short-distance sample pairs?

Because the variogram's behaviour near the origin — the nugget and the initial slope — determines the interpolation weights, and a uniform grid contains no pairs closer than its spacing.

Should I sample on the boundary?

Yes, if the boundary is part of the study area. Locations there have samples on one side only, so their error is systematically higher.

Can I fix a clustered dataset afterwards?

Partly. Declustering weights reduce the bias in summary statistics and kriging handles redundancy natively, but neither creates coverage where nobody measured.