Clipping Stockpile Boundaries from Vector Polygons

The volume came back 3.2 % higher than last month and nobody delivered any material. The surfaces are identical to within a centimetre on the hardstanding. The difference is that somebody re-traced the pile boundary, and their line sits forty centimetres further out than the previous one all the way round.

Boundary geometry is the second-largest discretionary input to a volume, after the base surface, and unlike the base it changes every time a human touches it. This page covers rasterizing a boundary so the result is reproducible, deciding what to do with the cells the polygon only partly covers, and the version discipline that stops the boundary being an unrecorded variable.

# Why a boundary offset costs so much volume

A stockpile does not end at a line; it tapers. The material sits at its angle of repose, typically 32–38° for crushed aggregate, so the surface near the toe rises roughly 0.7 m for every metre inward. Moving the boundary outward by $d$ metres adds a wedge whose cross-section is about $\tfrac{1}{2} d^2 \tan\theta$, multiplied by the perimeter.

For a pile with a 160 m perimeter and a 34° repose angle, a 0.4 m outward shift adds around 8.6 m³ — small. But the same shift also includes ground that sits below the base surface plane near the toe, and on a pile where the boundary was previously drawn inside the toe, the recovered slope area is the dominant term. In practice, boundary disagreements of half a metre routinely move real piles by two to four percent, which is larger than most clients’ tolerance.

Figure 1 — Why “it’s only forty centimetres” is not a defence. The wedge is small; what the boundary excludes is not.

# Minimal reproducible solution

Rasterize the polygon deterministically, and make the partial-cell decision explicit rather than inheriting a library default.

import numpy as np
import rasterio
from rasterio.features import rasterize
from shapely.geometry import mapping


def fractional_mask(polygon, transform, shape_hw: tuple[int, int],
                    subsample: int = 4) -> np.ndarray:
    """Per-cell coverage fraction of the polygon, in [0, 1].

    A binary mask forces every edge cell to be fully in or fully out, which
    on a coarse grid is the difference between two operators' volumes. A
    fractional mask computed by supersampling costs a few seconds and removes
    the discretisation from the answer entirely.
    """
    h, w = shape_hw
    fine = rasterize(
        [(mapping(polygon), 1)],
        out_shape=(h * subsample, w * subsample),
        transform=transform * transform.scale(1 / subsample, 1 / subsample),
        fill=0, all_touched=False, dtype="uint8")
    # Average each subsample block back down to the target grid.
    return (fine.reshape(h, subsample, w, subsample)
                .mean(axis=(1, 3))
                .astype("float32"))


def volume_with_fractional_mask(surface: np.ndarray, base: np.ndarray,
                                coverage: np.ndarray, cell: float) -> float:
    """Integrate with partial cells weighted by their coverage."""
    dz = surface - base
    valid = np.isfinite(dz) & (coverage > 0)
    return float(np.nansum(np.where(valid, dz * coverage, 0.0)) * cell ** 2)

Supersampling at 4× means each cell is resolved into sixteen sub-cells, which puts the residual discretisation error below a quarter of a percent of a cell — far below anything else in the budget. Going to 8× costs four times the memory and buys nothing measurable.

The alternative approaches both have a bias. all_touched=False includes only cells whose centre is inside, which systematically under-counts by about half a cell all the way round. all_touched=True includes any cell the polygon touches, which over-counts by a similar amount. On a 10 cm grid over a 160 m perimeter, either bias is worth roughly 8 m² of footprint — small on a tall pile, significant on a shallow one.

# Where to draw the line

The polygon is the single largest lever on a stockpile volume, and the best place to draw it depends on what the number is for.

Figure 3 — The right boundary depends on whether the deliverable is a volume or a difference.

# Edge-case matrix

from shapely.validation import make_valid


def sanitise_boundary(polygon, raster_crs, polygon_crs, raster_bounds):
    """Every check that has to pass before a polygon may define a volume."""
    if str(polygon_crs) != str(raster_crs):
        raise ValueError(f"boundary is {polygon_crs}, raster is {raster_crs}")
    if not polygon.is_valid:
        fixed = make_valid(polygon)
        if abs(fixed.area - polygon.area) / max(polygon.area, 1e-9) > 0.001:
            raise ValueError("repairing the boundary changed its area by >0.1 %")
        polygon = fixed
    minx, miny, maxx, maxy = polygon.bounds
    rminx, rminy, rmaxx, rmaxy = raster_bounds
    if minx < rminx or miny < rminy or maxx > rmaxx or maxy > rmaxy:
        raise ValueError("boundary extends beyond the surface raster")
    return polygon

# Verification snippet

The boundary’s contribution to the answer can be measured directly by perturbing it, which turns “is the boundary good enough” into a number in the same units as the deliverable.

import numpy as np


def boundary_sensitivity(polygon, surface, base, transform, shape_hw,
                         cell: float, offsets=(-0.5, -0.25, 0.0, 0.25, 0.5)) -> dict:
    """Volume as a function of buffering the boundary in and out."""
    out = {}
    for d in offsets:
        buffered = polygon.buffer(d)
        cov = fractional_mask(buffered, transform, shape_hw)
        out[d] = volume_with_fractional_mask(surface, base, cov, cell)

    v0 = out[0.0]
    per_metre = (out[0.5] - out[-0.5]) / 1.0
    return {"volumes_m3": out,
            "m3_per_metre_of_offset": float(per_metre),
            "percent_per_0_25m": float(abs(per_metre * 0.25 / v0) * 100)}

Reporting “each 25 cm of boundary uncertainty is worth 1.8 % of this volume” does two useful things. It tells the operator how carefully to trace, and it gives the client a concrete reason why the boundary is versioned rather than re-drawn each month.

Figure 2 — The discretisation, and the cheapest way to stop caring about it.

# Versioning the boundary so it stops being a variable

The technical fix above removes the grid from the answer. It does not remove the operator, and the operator is the larger source of variance. The organisational fix is to treat the boundary as a versioned asset with the same seriousness as the base surface.

Three conventions work. One polygon per pile, stored with the site, not per survey — a new survey reads the existing geometry rather than producing its own. A valid-from date on each version, so that when a pile genuinely changes footprint the change is a new version with a reason attached rather than an untracked edit. A diff in the run record: if the boundary used differs from the previous run’s, the volume record says so and reports the sensitivity figure alongside, so a reviewer sees immediately that the change is explained.

import hashlib
import json
from pathlib import Path


def boundary_version(path: str) -> dict:
    """Identify the exact boundary geometry a volume was computed against."""
    raw = Path(path).read_bytes()
    return {"path": path,
            "sha256_16": hashlib.sha256(raw).hexdigest()[:16],
            "feature_count": len(json.loads(raw)["features"])}


def compare_to_previous(current: dict, previous: dict | None) -> str | None:
    """Return a human-readable note when the boundary changed between runs."""
    if previous is None:
        return "first run for this pile — boundary established"
    if current["sha256_16"] != previous["sha256_16"]:
        return (f"boundary changed since the previous run "
                f"({previous['sha256_16']} → {current['sha256_16']}): "
                "volume comparison is not like-for-like")
    return None

The note is what makes the discipline stick. A monitoring report that silently used a different boundary is indistinguishable from one that did not; a report carrying that sentence in bold is one where nobody has to reconstruct what happened.

# When to escalate

• The pile genuinely changed shape and the old boundary no longer contains it. Cut a new version with a reason, and expect the series to step. Do not quietly enlarge the old one.
• Two piles have merged. The geometry question has become a commercial one — which material belongs to which stock — and no automated answer is appropriate. Escalate before processing.
• The boundary sensitivity exceeds the client’s tolerance on its own. That is a shallow pile with a long perimeter, where a volume is a weak measurement. Say so, and consider whether an area-and-depth product answers the question better.

# Related

• Computing Volumes and Stockpiles in Python
• Reconciling Volume Differences Between Flights
• Estimating Volume Uncertainty from Surface Error

← Computing Volumes and Stockpiles in Python