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      58.7 m      234 m         469 m        811 m
regular     44.6 m      188 m         297 m      1,100 m
clustered  215.7 m    1,672 m       3,644 m      5,021 m

A regular grid is 24% more accurate than random placement at the same cost. Clustered sampling β€” which is what a survey following roads, rivers or convenient access produces β€” is nearly four 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.

Random, regular grid and clustered sampling of the same 500 points, with RMSE of 58.7, 44.6 and 215.7 metres.
Same budget, same interpolator, same truth. The placement is worth a factor of five.

Step-by-step solution

1. Design against the gap, not the count

"500 samples" says nothing. 500 samples in five clusters left a 95th-percentile gap of 3,644 m and produced an RMSE of 215.7 m β€” worse than 100 well-spread samples would have.

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.6 m against random's 58.7 m, because it has no gaps larger than its spacing. It also had the largest single gap (1,100 m against random's 811 m) β€” an artefact of the grid not fitting the irregular study boundary, leaving a strip along one edge.

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, error rose from 18.9 m within 100 m of a sample to 334.3 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.

Distribution of distance to the nearest sample for random, grid and clustered designs, with the clustered design's 95th percentile at 3,644 metres.
The clustered design's typical location is 1.7 km from any observation. That is what an RMSE of 215 m is made of.

Code examples

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

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.

Explanation

Why clustered sampling fails so badly

The clustered design had 500 samples and a mean gap of 1,672 m, against random's 234 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 the well-spread random design measured here, LOO was 13% pessimistic (60.8 m against a true 54.0 m); under clustering the sign flips and the magnitude grows. 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 58.4 m to 43.3 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 a range of 4,537 m, that suggests 1.5–2 km spacing β€” and the 500-point designs above, at 457 m ideal spacing, are comfortably denser than needed for that structure. Their remaining error is dominated by variation at scales shorter than the range.

A grid design supplemented with nested short-distance pairs at six sites, costing five percent of the budget.
Five percent of the budget buys the entire short-lag half of the variogram.

Edge cases or notes

  • 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.

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 3.6 km gap.

Is a regular grid better than random sampling?

For coverage, yes: 44.6 m RMSE against 58.7 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 215.7 m RMSE against 58.7 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.