Declustering Weights with Cell and Polygon Methods

TL;DR: Weight each observation by 1 / (n_occupied_cells * n_in_its_cell), average that over eight shifted grid origins, and sweep the cell size to find the extremum — minimum when the dense sampling sits in the highs. Or skip the tuning parameter entirely and use shapely.voronoi_polygons(...) areas clipped to the domain. On the worked block below both turn a 64% overstatement into a sub-1% error.

Why This Matters #

Sampling is almost never a random draw from the population you want to describe. Drill holes follow the ore, monitoring stations follow the pollution, soil pits follow the interesting soil. The result is a sample whose composition is wrong even when every individual measurement is perfect: an arithmetic mean gives one vote per observation, and the well-sampled ground has thousands of votes while the sparsely sampled ground has ten. Declustering is the correction — a set of weights $w_i$ summing to one that give each observation influence in proportion to the ground it represents rather than in proportion to how often somebody visited it.

Two things are worth stating before any code. First, declustering fixes the statistic, not the map. A declustered mean, variance and histogram are what feed simple kriging, normal score transforms and sequential simulation; the kriging weights themselves come from the variogram and are untouched. Second, declustering is one of several tools in Sampling Bias Mitigation, and it addresses only spatial over-representation. It cannot repair a value-dependent sampling rule such as reporting assays only above a cut-off. The wider treatment of what areal supports do to statistics sits in Core Concepts of Spatial Statistics & Geostatistics, and the closely related question of how the units themselves change an answer is covered in Spatial Scale & the Modifiable Areal Unit Problem.

The Two Weight Rules #

Both classical methods answer the same question — how much ground does observation $i$ speak for? — and differ only in how they measure it.

Cell declustering lays a square grid of side $c$ over the domain. Let $n_c$ be the number of cells containing at least one observation and let $n_{k(i)}$ be the number of observations sharing the cell that contains observation $i$. Then

$$w_i = \frac{1}{n_c \, n_{k(i)}}, \qquad \sum_{i=1}^{n} w_i = 1 .$$

Every occupied cell receives exactly $1/n_c$ of the total weight and splits it evenly among its occupants, so a cell holding twenty infill holes gives each of them a twentieth of what a lone hole in an empty quarter of the map receives. The declustered mean is $\bar{z}_D = \sum_i w_i z_i$.

Polygonal declustering replaces the grid with the Voronoi tessellation. Let $V_i$ be the set of points closer to observation $i$ than to any other, let $D$ be the study domain, and set

$$w_i = \frac{A_i}{\sum_j A_j}, \qquad A_i = \operatorname{area}(V_i \cap D) .$$

This is the nearest-neighbour or “polygons of influence” estimate of the mean. It has no tuning parameter, which is its great attraction, and it depends entirely on $D$, which is its great weakness — the outermost observations own whatever ground you decide to clip to, and there is no data-driven way to decide that.

Environment and Version Pinning #

Nothing here needs a geostatistics package: the weights are arithmetic over a spatial index. The Voronoi construction comes from shapely 2.x, which wraps the GEOS implementation directly.

pip install "numpy==2.1.3" "pandas==2.2.3" "shapely==2.0.6" \
            "geopandas==1.0.1" "scipy==1.14.1" "pyproj==3.7.0"

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely import voronoi_polygons
from shapely.geometry import MultiPoint, box

Step-by-Step Implementation #

1. Build a preferentially sampled block #

The example is a 1 km square gold block with a smooth background trend and one Gaussian high-grade anomaly. Sixty-four exploration holes sit on a nominal 125 m grid; eighty-five infill holes are drilled into the anomaly. Because the field is synthetic, the true area mean is available as the reference every method is judged against.

rng = np.random.default_rng(2024)
SIDE = 1000.0
BLOCK = box(0, 0, SIDE, SIDE)

def true_grade(x, y):
    """Background trend plus one Gaussian high-grade anomaly, in g/t."""
    return 1.10 + 0.0009 * y + 4.20 * np.exp(
        -(((x - 720.0) ** 2 + (y - 680.0) ** 2) / (2 * 165.0 ** 2))
    )

# 64 exploration holes on a jittered 125 m grid
g = (np.arange(8) + 0.5) * (SIDE / 8)
gx, gy = np.meshgrid(g, g)
ex = gx.ravel() + rng.uniform(-18, 18, 64)
ey = gy.ravel() + rng.uniform(-18, 18, 64)

# infill holes drawn around the anomaly, kept inside the block
ix = rng.normal(720, 105, 86)
iy = rng.normal(680, 105, 86)
keep = (ix > 5) & (ix < SIDE - 5) & (iy > 5) & (iy < SIDE - 5)
ix, iy = ix[keep], iy[keep]

x = np.concatenate([ex, ix])
y = np.concatenate([ey, iy])
z = np.maximum(true_grade(x, y) + rng.normal(0, 0.22, x.size), 0.05)

holes = gpd.GeoDataFrame({"au_gt": z},
                         geometry=gpd.points_from_xy(x, y),
                         crs="EPSG:32633")

gg = np.linspace(0.5, SIDE - 0.5, 1000)
GX, GY = np.meshgrid(gg, gg)
truth = float(true_grade(GX, GY).mean())

naive = float(holes["au_gt"].mean())
print(f"n = {len(holes)}   naive mean = {naive:.4f} g/t")
print(f"true block mean = {truth:.4f} g/t  ({100 * (naive - truth) / truth:+.1f}%)")
print(f"grid holes {holes['au_gt'][:64].mean():.3f} g/t   "
      f"infill holes {holes['au_gt'][64:].mean():.3f} g/t")

n = 149   naive mean = 3.6454 g/t
true block mean = 2.2182 g/t  (+64.3%)
grid holes 2.143 g/t   infill holes 4.776 g/t

A 64% overstatement is not an exotic failure; it is the ordinary consequence of drilling where the grade is. Note the diagnostic in the third line — the background survey already estimates the block correctly on its own, and the naive mean is simply the infill dragging it upward.

2. Cell declustering weights, averaged over grid origins #

The one-line rule is $w_i = 1/(n_c\,n_{k(i)})$. The one refinement that matters is the grid origin: a single fixed origin makes the answer depend on an arbitrary choice, so compute the weights for several origins shifted by fractions of the cell size and average the weight vectors.

def cell_weights(x, y, cell_size, n_offsets=8):
    """Cell declustering weights averaged over `n_offsets` grid origins.

    Each observation gets 1 / (occupied cells * observations in its cell).
    Returns a vector summing to exactly 1.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    w = np.zeros(x.size)
    for k in range(n_offsets):
        shift = cell_size * k / n_offsets
        col = np.floor((x + shift) / cell_size).astype(np.int64)
        row = np.floor((y + shift) / cell_size).astype(np.int64)
        key = col * 1_000_003 + row          # safe while |row| < 500_000
        _, inverse, counts = np.unique(key, return_inverse=True,
                                       return_counts=True)
        w += 1.0 / (counts.size * counts[inverse])
    return w / w.sum()


w170 = cell_weights(x, y, 170.0)
print(f"weights sum to {w170.sum():.6f}")
print(f"min w = {w170.min():.5f}   max w = {w170.max():.5f}"
      f"   ratio = {w170.max() / w170.min():.1f}")
print(f"declustered mean at 170 m = {float(w170 @ z):.4f} g/t")

weights sum to 1.000000
min w = 0.00091   max w = 0.01968   ratio = 21.6
declustered mean at 170 m = 2.2187 g/t

The integer key trick replaces a slow np.unique(..., axis=0) on stacked columns with a single pass over a 1-D array, which matters once you are sweeping fourteen cell sizes across eight origins on a hundred thousand holes.

3. Sweep the cell size and read the extremum #

No cell size is correct a priori, so compute the declustered mean across a range spanning from far below to far above the sampling spacing and look at the shape of the curve. At very fine cells almost every observation occupies its own cell and the weights return to equal; at very coarse cells everything shares a handful of cells and the weights return to equal again. Between the two the curve turns, and the turning point is the declustered estimate.

sizes = [25, 50, 75, 100, 125, 150, 170, 200, 250, 300, 400, 500, 700, 1000]
rows = []
for cs in sizes:
    w = cell_weights(x, y, cs)
    key0 = (np.floor(x / cs).astype(np.int64) * 1_000_003
            + np.floor(y / cs).astype(np.int64))
    rows.append({"cell_m": cs,
                 "n_cells": int(np.unique(key0).size),
                 "dec_mean": float(w @ z),
                 "n_eff": float(1.0 / np.sum(w ** 2))})

sweep = pd.DataFrame(rows)
print(sweep.to_string(index=False, formatters={
    "dec_mean": "{:.4f}".format, "n_eff": "{:.1f}".format}))

best = sweep.loc[sweep["dec_mean"].idxmin()]
print(f"\nminimum {best.dec_mean:.4f} g/t at a {best.cell_m:.0f} m cell")

 cell_m  n_cells  dec_mean  n_eff
     25      132    3.4139  142.8
     50      102    3.0393  126.7
     75       85    2.6615  108.0
    100       74    2.4102   93.3
    125       64    2.2562   83.5
    150       49    2.2268   80.8
    170       36    2.2187   79.6
    200       25    2.2360   77.8
    250       16    2.2722   78.0
    300       16    2.2673   73.2
    400        9    2.3537   72.4
    500        4    2.3815   72.3
    700        4    2.5741   81.3
   1000        1    2.7652   88.0

minimum 2.2187 g/t at a 170 m cell

The minimum lands at 2.2187 g/t against a truth of 2.2182 — a bias of two hundredths of a per cent, from a starting point of plus sixty-four per cent. The cell size that finds it, 170 m, is close to the 125 m spacing of the background survey, which is the standard heuristic: the right cell is one that holds roughly one background observation, so that the sparse ground is represented once and the dense ground is represented once as well.

4. Polygonal weights from clipped Voronoi cells #

The polygonal method needs no sweep. It needs a domain polygon, and it needs the pairing between Voronoi cells and observations to be recovered by geometry rather than by position, because voronoi_polygons does not return its output in input order.

cells = voronoi_polygons(MultiPoint(list(zip(x, y))), extend_to=BLOCK)
vor = gpd.GeoDataFrame(geometry=list(cells.geoms), crs=holes.crs)
vor["cell_id"] = np.arange(len(vor))

# extend_to only guarantees the diagram covers the block; the outer cells
# still run past it, so clip explicitly.
vor["geometry"] = vor.geometry.intersection(BLOCK)

pair = gpd.sjoin(holes, vor[["cell_id", "geometry"]],
                 how="left", predicate="within")["cell_id"].to_numpy()
assert not np.isnan(pair).any(), "a hole fell outside every polygon"
pair = pair.astype(int)

areas = vor.geometry.area.to_numpy()[pair]
w_poly = areas / areas.sum()

print(f"clipped area total = {areas.sum():,.0f} m2 "
      f"(block = {SIDE ** 2:,.0f} m2)")
print(f"areas: min {areas.min():,.0f}   median {np.median(areas):,.0f}   "
      f"max {areas.max():,.0f} m2")
print(f"polygonal declustered mean = {float(w_poly @ z):.4f} g/t")

clipped area total = 1,000,000 m2 (block = 1,000,000 m2)
areas: min 51   median 3,847   max 18,525 m2
polygonal declustered mean = 2.2031 g/t

The total area recovering the block area exactly is the check that the clip worked: Voronoi cells tile the plane, so their clipped areas must tile the domain. If that total comes back short, a polygon failed to intersect; if it comes back long, you clipped to the wrong thing.

5. Show what the clip is worth #

The single most important number in polygonal declustering is not the mean but the sensitivity of the mean to the domain you clipped to. Compute it by re-clipping the same tessellation to the convex hull of the observations, which is what you implicitly get when nobody supplies a boundary.

hull_cells = vor.geometry.intersection(holes.union_all().convex_hull)
a_hull = hull_cells.area.to_numpy()[pair]
w_hull = a_hull / a_hull.sum()

def report(name, w):
    mu = float(w @ z)
    sd = float(np.sqrt(w @ (z - mu) ** 2))
    print(f"{name:<26} {mu:6.4f}  {100 * (mu - truth) / truth:+6.2f}%  "
          f"sd {sd:.4f}  n_eff {1 / np.sum(w ** 2):5.1f}")

report("naive (equal weights)", np.full(len(z), 1 / len(z)))
report("cell, 170 m, 8 origins", w170)
report("polygonal, block clip", w_poly)
report("polygonal, hull clip", w_hull)

naive (equal weights)      3.6454  +64.34%  sd 1.6756  n_eff 149.0
cell, 170 m, 8 origins     2.2187   +0.02%  sd 1.2661  n_eff  79.6
polygonal, block clip      2.2031   -0.68%  sd 1.2119  n_eff  80.0
polygonal, hull clip       2.3087   +4.08%  sd 1.2823  n_eff  77.9

Interpreting the Output #

Three quantities carry the result. The declustered mean is the headline: 2.2187 or 2.2031 g/t depending on method, against a naive 3.6454. The declustered standard deviation, $s_D^2 = \sum_i w_i (z_i - \bar{z}_D)^2$, falls from 1.676 to about 1.27 because the over-represented high values stop inflating the spread; this matters as much as the mean, since the declustered histogram is what a normal score transform and the reference distribution for simulation are built from. The effective sample size, $n_{\text{eff}} = 1/\sum_i w_i^2$, drops from 149 to about 80. That last number is the honest statement of what the survey bought: roughly eighty independent pieces of information about the block mean, not a hundred and forty-nine.

What good looks like is a sweep curve with a clear interior turning point at a cell size within a factor of two of the background sampling spacing, and two methods agreeing to within a couple of per cent. Here they differ by 0.7%, which is noise. What is not noise is the fourth row: clipping the same Voronoi cells to the sample convex hull instead of the block moves the answer by 4.8 percentage points relative to the block clip, from 2.2031 to 2.3087. Nothing about the data changed. The hull excludes the barren edges of the block that no hole quite reached, so the low-grade perimeter loses the area it deserves.

Warning signs are worth naming. A sweep with no interior extremum — monotone across the whole range — usually means the sampling is not actually clustered at any resolvable scale, or that the range you swept does not reach past the spacing of the sparse data. A minimum that lands at the smallest cell you tried is the same complaint from the other side. And a very jagged sweep, one that jumps by several per cent between adjacent cell sizes, means you are seeing origin sensitivity rather than the declustering signal; raise n_offsets.

Critical Best Practices #

Decide the direction of the sweep before you run it #

The rule is that the extremum you take must match where the extra sampling went: minimum when the dense groups sit in high values, maximum when they sit in low values. This is a statement about the sampling campaign, not about the data, and it has to be settled from the drilling history, the monitoring programme design or a simple map before the sweep is computed. Reading the sweep first and then choosing the extremum that produces the number you expected is how declustering becomes a way of manufacturing an answer.

Always average over shifted grid origins #

With a single fixed origin the declustered mean is a function of where you happened to put the corner of the grid. On this dataset the single-origin result at a 150 m cell is 2.297 g/t and at 170 m it is 2.258, while the eight-origin averages are 2.227 and 2.219 — the jitter is comparable to the effect being measured. Eight offsets is the usual default and costs eight passes over an integer array, which is nothing.

Clip Voronoi cells to the real domain, never to the sample hull #

The convex hull of the observations is not the study area; it is a summary of where somebody drilled, and using it deletes exactly the sparsely sampled perimeter that declustering exists to protect. The 4.1% shift shown above is on a well-behaved rectangular block, and the effect is far larger when the sampling is concentrated in one corner of an irregular licence. If you genuinely do not have a boundary polygon, prefer cell declustering, which never asks the question.

Report the effective sample size with the mean #

Weighted means have larger standard errors than equally weighted ones, and $n_{\text{eff}} = 1/\sum w_i^2$ quantifies it. Falling from 149 to 79.6 means a confidence interval built from the raw count is roughly $\sqrt{149/79.6} \approx 1.37$ times too narrow. A declustered mean quoted without $n_{\text{eff}}$ invites exactly the over-confidence the correction was supposed to remove.

Decluster before the variogram, not after #

The declustered distribution is what the normal score transform maps to Gaussian, and a transform built from an inflated histogram distorts every subsequent structural analysis. Fit structure on declustered, transformed values — see Empirical Variogram Estimation for the estimator that consumes them. The weights themselves do not enter the variogram estimator, which pairs observations by lag, but a clustered dataset still biases short-lag bins towards the dense area, so read the lag counts before trusting a nugget.

Troubleshooting #

Next Steps #

Fold the weights into the wider pipeline described in Correcting Spatial Sampling Bias with GeoPandas, then carry the declustered distribution forward into Empirical Variogram Estimation so the structural model is fitted to a histogram that represents the ground rather than the drill plan.

Frequently Asked Questions #

Should I take the minimum or the maximum of the cell-size sweep?

It depends on where the extra sampling went. If the tight groups of observations sit in the high values, as they do when infill drilling follows a high-grade anomaly, the naive mean is inflated and the declustered mean dips to a minimum before rising again, so take the minimum. If the tight groups sit in low values, which happens when monitoring concentrates on clean reference sites, take the maximum. Decide this from the sampling history and a map before you look at the sweep, never by picking whichever extreme is more convenient.

Does declustering change the kriged map as well as the mean?

Not directly. Kriging weights come from the variogram and the local configuration of neighbours, so a declustering weight never enters the kriging system. Declustering fixes global statistics: the mean, the variance, the histogram used for normal score transforms and for simple kriging. Because simple kriging and sequential Gaussian simulation both take that global mean and that reference distribution as inputs, a biased histogram propagates into the map through the back door even though the kriging equations never saw the weights.

Which method should I use when the two disagree?

Prefer cell declustering when the study domain has a hard, well-defined edge and the sampling is a background survey with dense infill, because the cell weights ignore the boundary entirely. Prefer polygonal declustering when the sampling is irregular at every scale, so no single cell size describes it, and when the domain boundary is genuinely known. A disagreement of more than a few per cent is almost always the boundary: the polygonal weights are absorbing whatever you clipped to, so check that polygon before trusting either number.

How many observations do I lose by declustering?

Report the Kish effective sample size, one divided by the sum of the squared weights, which equals the count only when the weights are equal. In the worked example 149 observations fall to an effective 79.6 under cell declustering and 80.0 under polygonal weights, so roughly half the nominal precision was never real. That number belongs in the report next to the declustered mean, because a confidence interval computed from the raw count will be about forty per cent too narrow.

Related

• Correcting Spatial Sampling Bias with GeoPandas — the surrounding workflow these weights slot into
• Spatial Scale & the Modifiable Areal Unit Problem — why the cell size changes the answer at all
• Empirical Variogram Estimation — what the declustered distribution feeds next

← Back to Sampling Bias Mitigation