Computing Volumes and Stockpiles in Python

Volume is the commercial output of drone surveying. A quarry pays for a monthly figure it will settle invoices against; a contractor pays for a cut-and-fill balance that determines how many lorry movements the job needs. The arithmetic is trivial — sum the height difference over an area — and the number is routinely wrong by more than the client’s tolerance, for reasons that have nothing to do with arithmetic.

Three choices dominate the answer, and none of them is a numerical method. Where the base surface came from. Where the boundary was drawn. Whether the two surveys being compared share a vertical datum. A pipeline that automates the integration while leaving those three to a technician’s judgement has automated the easy tenth of the problem.

This page covers volume as a scripted, auditable stage: the two integration methods and when the difference matters, boundary handling that survives a disagreement with the client, and an error budget that separates the error that averages away from the error that does not. It builds directly on the classified cloud from classifying point clouds with PDAL and Python.

Audience and prerequisites. Python 3.10+, a classified point cloud or a pair of rasterized surfaces in a projected CRS with metre units, and a boundary polygon. Volumes computed from surfaces in geographic coordinates are wrong by the cosine of the latitude, silently.

# Prerequisites

# Conceptual architecture

Every volume is the same integral: the height difference between an upper surface and a lower one, integrated over a region. What varies is where each of the three inputs comes from.

The upper surface is almost always the current survey’s DSM, restricted to the pile. The lower surface — the base — is the part that is genuinely ambiguous: it can be a fitted plane through the toe of the pile, the terrain model from a survey taken before the material arrived, a fixed elevation the client specifies, or a triangulation of hand-picked toe points. These four choices routinely differ by several percent on the same pile. The region is a polygon whose edge sits somewhere on the slope between pile and ground, and moving it by half a metre on a pile with a 34° angle of repose changes the volume by roughly the perimeter times a third of a square metre.

Figure 1 — Where the disagreement actually lives. The integration is identical in all eight cases.

# Step 1: Rasterize both surfaces onto one grid

The single most common source of a quietly wrong volume is two surfaces on grids that do not align. Subtracting a 5 cm DSM from a 10 cm base by resampling one onto the other introduces an interpolation error that is small everywhere and systematically signed on slopes — which is exactly where a stockpile has all its area.

import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling


def align_to(reference_path: str, other_path: str) -> tuple[np.ndarray, dict]:
    """Resample `other` onto the reference grid exactly, returning the array.

    Bilinear is correct here: both rasters are continuous surfaces. What is
    NOT correct is resampling the finer raster onto the coarser one, which
    discards the resolution that was paid for — always resample the base up
    to the survey grid, never the survey down.
    """
    with rasterio.open(reference_path) as ref, rasterio.open(other_path) as src:
        if src.res[0] < ref.res[0]:
            raise ValueError("the base surface is finer than the survey grid; "
                             "resample the base, not the survey")
        dst = np.full((ref.height, ref.width), np.nan, dtype="float32")
        reproject(
            source=rasterio.band(src, 1), destination=dst,
            src_transform=src.transform, src_crs=src.crs,
            dst_transform=ref.transform, dst_crs=ref.crs,
            resampling=Resampling.bilinear, dst_nodata=np.nan)
        profile = ref.profile
    return dst, profile

Raising on the wrong-direction resample is worth the three lines. It is a mistake that produces a plausible number and is invisible in every subsequent check.

# Step 2: Build the boundary mask, and keep its edge explicit

A boundary polygon rasterized with default settings includes a cell when its centre falls inside the polygon. On a 10 cm grid that is a sub-decimetre decision per cell and irrelevant; on a 50 cm grid over a pile with a 60 m perimeter it is a percent of the volume.

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


def boundary_mask(polygon, transform, shape_hw: tuple[int, int],
                  *, all_touched: bool = False) -> np.ndarray:
    """Rasterize the boundary, and report the area the rasterization gained or lost.

    `all_touched=True` includes any cell the polygon touches, which
    over-counts; the default includes cells whose centre is inside, which
    under-counts. Neither is wrong — but the difference must be small
    relative to the tolerance, and that has to be checked, not assumed.
    """
    mask = rasterize([(mapping(polygon), 1)], out_shape=shape_hw,
                     transform=transform, fill=0, all_touched=all_touched,
                     dtype="uint8").astype(bool)
    cell_area = abs(transform.a * transform.e)
    raster_area = mask.sum() * cell_area
    disagreement = abs(raster_area - polygon.area) / polygon.area
    if disagreement > 0.01:
        raise ValueError(
            f"rasterized boundary differs from the polygon by {disagreement:.1%} — "
            "the grid is too coarse for this boundary")
    return mask

The check turns an invisible discretisation into an explicit failure. A pile small enough that its boundary cannot be represented on the chosen grid is a pile that needs a finer grid, and the pipeline should say so rather than returning a number.

# Step 3: Integrate, with the two error terms kept apart

import numpy as np


def volume_with_uncertainty(surface: np.ndarray, base: np.ndarray,
                            mask: np.ndarray, cell_size: float,
                            *, sigma_random: float = 0.03,
                            sigma_systematic: float = 0.02) -> dict:
    """Prism volume over the mask, with random and systematic error separated.

    The distinction is the whole point. Random per-cell error averages down
    with the square root of the cell count; a systematic offset — a datum
    error, a calibration bias, a base surface that is uniformly low — does
    not average at all and scales with the full footprint area.
    """
    dz = np.where(mask, surface - base, np.nan)
    finite = np.isfinite(dz)
    n = int(np.count_nonzero(finite))
    if n == 0:
        raise ValueError("no finite cells inside the boundary")

    cell_area = cell_size ** 2
    volume = float(np.nansum(dz) * cell_area)
    cut = float(np.nansum(np.where(dz < 0, dz, 0.0)) * cell_area)
    fill = float(np.nansum(np.where(dz > 0, dz, 0.0)) * cell_area)

    sigma_v_random = sigma_random * np.sqrt(n) * cell_area
    sigma_v_systematic = sigma_systematic * n * cell_area
    sigma_total = float(np.hypot(sigma_v_random, sigma_v_systematic))

    return {
        "volume_m3": volume, "cut_m3": cut, "fill_m3": fill,
        "cells": n, "area_m2": n * cell_area,
        "sigma_random_m3": float(sigma_v_random),
        "sigma_systematic_m3": float(sigma_v_systematic),
        "sigma_total_m3": sigma_total,
        "relative_uncertainty": sigma_total / abs(volume) if volume else float("nan"),
        "void_fraction": 1.0 - n / max(int(np.count_nonzero(mask)), 1),
    }

Reporting void_fraction alongside the volume matters more than it looks. A pile with fifteen percent NoData inside the boundary has had that fraction of its volume silently excluded, and a client comparing month to month will see the pile shrink when in fact the reconstruction got worse. Gate on it.

Figure 2 — Why a single “±3 cm accuracy” figure is not an uncertainty. The term that dominates is the one that never appears in a sensor specification.

# Step 4: Fit a base plane that does not move between surveys

Where no prior terrain model exists, the base is usually a plane through the pile’s toe. Fitting it from the current survey is convenient and wrong for monitoring: a pile that grows outward pushes its own toe onto new ground, the fitted plane tilts, and the volume changes without any material moving.

import numpy as np


def fit_toe_plane(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> tuple:
    """Least-squares plane through toe points, returned as (a, b, c) for z = ax+by+c.

    Centre the coordinates before solving. A UTM easting of 500 000 in a
    normal-equations solve costs most of the available floating-point
    precision and the resulting plane can be visibly tilted.
    """
    x0, y0 = float(np.mean(x)), float(np.mean(y))
    A = np.column_stack([x - x0, y - y0, np.ones_like(x)])
    coef, *_ = np.linalg.lstsq(A, z, rcond=None)
    a, b, c = (float(v) for v in coef)
    return a, b, c - a * x0 - b * y0


def plane_surface(a: float, b: float, c: float, transform, shape_hw) -> np.ndarray:
    """Evaluate the fitted plane over a raster grid."""
    h, w = shape_hw
    cols, rows = np.meshgrid(np.arange(w), np.arange(h))
    xs = transform.c + (cols + 0.5) * transform.a
    ys = transform.f + (rows + 0.5) * transform.e
    return (a * xs + b * ys + c).astype("float32")

Fit it once, store the coefficients with the site, and reuse them. Reconciling volume differences between flights covers what happens when this discipline is missing, and choosing a base surface for stockpile volumes covers how to pick among the four candidates in the first place.

# Step 5: Emit a volume record, not a number

A volume delivered as a single figure in an email is unreviewable. Six weeks later, when the client’s quantity surveyor disputes it, nobody can reconstruct which base surface, which boundary version or which grid produced it. The remedy costs one function: emit every input alongside the answer, and make the record the deliverable.

import hashlib
import json
from datetime import datetime, timezone
from pathlib import Path


def volume_record(result: dict, *, site: str, boundary_path: str,
                  surface_path: str, base_mode: str,
                  base_params: dict, cell_size: float) -> dict:
    """Everything needed to reproduce this volume, in one JSON document."""
    def digest(path: str) -> str:
        h = hashlib.sha256()
        with open(path, "rb") as fh:
            for chunk in iter(lambda: fh.read(1 << 20), b""):
                h.update(chunk)
        return h.hexdigest()[:16]

    return {
        "site": site,
        "computed_at": datetime.now(timezone.utc).isoformat(),
        "volume_m3": round(result["volume_m3"], 1),
        "uncertainty_m3": round(result["sigma_total_m3"], 1),
        "uncertainty_split": {
            "random_m3": round(result["sigma_random_m3"], 1),
            "systematic_m3": round(result["sigma_systematic_m3"], 1),
        },
        "area_m2": round(result["area_m2"], 1),
        "void_fraction": round(result["void_fraction"], 4),
        "inputs": {
            "surface": {"path": surface_path, "sha256_16": digest(surface_path)},
            "boundary": {"path": boundary_path, "sha256_16": digest(boundary_path)},
            "base_mode": base_mode, "base_params": base_params,
            "cell_size_m": cell_size,
        },
    }

Hashing the inputs rather than merely naming them is what makes the record worth keeping. A path is a promise that the file has not changed; a digest is a fact. When a monthly series shows a step change, comparing digests answers “did the boundary change” in a second, and that single question resolves most disputes.

Two further conventions pay for themselves on a monitoring contract. Round the volume to a precision the uncertainty supports — reporting 18,472.3 m³ next to an uncertainty of ±210 m³ invites a false sense of precision, and rounding to the nearest ten is both honest and easier to read. Keep the boundary under version control with the rest of the site configuration, so a change to it appears in a diff rather than in a file share.

The same record feeds the accuracy statement in the deliverable. A client tolerance expressed as a percentage — “within two percent” — can be tested directly against uncertainty_m3 / volume_m3, and a run that cannot meet it should fail loudly at computation time rather than being discovered during settlement.

# Parameter deep-dive

# Verification and output inspection

Two checks catch nearly everything: a closed-form comparison on a synthetic shape, and a resolution sweep on the real one.

import numpy as np


def verify_against_cone(cell_size: float = 0.10, radius: float = 12.0,
                        height: float = 6.0, tol: float = 0.005) -> None:
    """Integrate a synthetic cone and compare against its analytic volume.

    This is the only test that validates the integration itself rather than
    the data. A pipeline that cannot reproduce a cone to half a percent has
    a bug in the mask, the cell area, or the NoData handling.
    """
    n = int(2 * radius / cell_size) + 4
    xs = (np.arange(n) - n / 2) * cell_size
    gx, gy = np.meshgrid(xs, xs)
    r = np.hypot(gx, gy)
    surface = np.where(r <= radius, height * (1 - r / radius), 0.0)
    base = np.zeros_like(surface)
    mask = r <= radius

    got = float(np.nansum(np.where(mask, surface - base, 0.0)) * cell_size ** 2)
    expected = np.pi * radius ** 2 * height / 3
    rel = abs(got - expected) / expected
    assert rel < tol, f"cone volume off by {rel:.3%} — integration is wrong"


def resolution_sweep(compute, sizes=(0.05, 0.10, 0.20, 0.40)) -> dict:
    """Volume at several grid resolutions; a stable pipeline varies little."""
    out = {s: compute(s)["volume_m3"] for s in sizes}
    spread = (max(out.values()) - min(out.values())) / abs(np.mean(list(out.values())))
    out["relative_spread"] = float(spread)
    return out

A resolution sweep spreading more than about one percent means the grid is interacting with the boundary or with voids, and the reported volume is a function of an arbitrary parameter. That is worth knowing before a client finds it.

Figure 3 — Five stages, five choices, each of which moves the number.

# Troubleshooting

The volume changed between two flights and nothing moved on site. The base surface was re-derived from each survey. Fix it once and reuse it. If the base is a prior DTM, confirm both surveys share a vertical datum first.

The volume is negative. The surfaces are the wrong way round, or the base sits above the pile because it was fitted to points on higher ground outside the toe. Print the mean of the difference before integrating.

Volumes differ by several percent between two operators. Almost always the boundary. Two people tracing the toe of a pile disagree by a metre routinely, and on a typical pile that is two to four percent. Store the polygon, version it, and reuse it rather than re-tracing.

The number is stable but the client’s ground survey disagrees by eight percent. Check for a systematic vertical offset between the two datums, and check what the client’s method measures — a truck-count estimate and a surface integral measure genuinely different things, and neither is wrong.

Fifteen percent of the pile is NoData. The reconstruction failed on the pile surface, usually because dark wet material has no texture. The volume is understated by roughly that fraction. Re-fly with different lighting rather than interpolating across it.

Cut and fill both look implausibly large on a flat site. The two surfaces are misaligned horizontally, so every slope produces a paired cut and fill. Co-register before differencing, as in aligning two epochs with ICP before differencing.

# Related

• Choosing a Base Surface for Stockpile Volumes
• Estimating Volume Uncertainty from Surface Error
• Clipping Stockpile Boundaries from Vector Polygons
• Reconciling Volume Differences Between Flights
• Generating DSM and DTM from Point Clouds with PDAL

← Point Cloud Processing & 3D Deliverables