Article

Measuring Throughput of a Raster Batch Pipeline

Time each unit individually, aggregate into throughput and percentiles rather than a mean, and sample peak memory on a background thread throughout. Reporting all three together is what stops a change that doubles tiles-per-minute and triples resident memory being recorded as a win. It builds on the Benchmarking Spatial Batch Pipelines guide, part of the broader Spatial Batch Processing & Async Workflows reference.

Prerequisites

  • Python 3.10 or later
  • pip install psutil rasterio numpy
  • A corpus that spans the real size distribution, per the parent guide

Why the Mean Lies

A batch’s wall clock is set by its slowest units, and geospatial inputs are famously unequal.

Which statistic predicts the wall clock The mean is dragged down by many fast tiles and underestimates the batch. The median is worse. The 95th percentile tracks the units that set the finish time. The maximum identifies the single worst input, which is often a data problem rather than a code one. mean per-tile duration dominated by the many small tiles underestimates badly median per-tile duration worse still — ignores the tail entirely not useful here 95th percentile tracks the units that set the finish time the number to quote maximum names the single worst input often a data problem

Quoting the mean is how a change that speeds up small tiles and slows large ones gets recorded as an improvement while the batch takes longer.

Complete Working Implementation

# bench/harness.py
from __future__ import annotations

import statistics
import threading
import time
from contextlib import contextmanager
from dataclasses import dataclass, field

import psutil


@dataclass
class Result:
    durations: list[float] = field(default_factory=list)
    peak_rss: int = 0
    wall: float = 0.0

    @property
    def throughput(self) -> float:
        return len(self.durations) / self.wall * 60 if self.wall else 0.0

    def percentile(self, q: float) -> float:
        if not self.durations:
            return 0.0
        ordered = sorted(self.durations)
        idx = min(int(q * len(ordered)), len(ordered) - 1)
        return ordered[idx]

    def summary(self) -> dict:
        return {
            "units": len(self.durations),
            "wall_s": round(self.wall, 2),
            "throughput_per_min": round(self.throughput, 1),
            "p50_ms": round(self.percentile(0.50) * 1000, 1),
            "p95_ms": round(self.percentile(0.95) * 1000, 1),
            "max_ms": round(max(self.durations, default=0) * 1000, 1),
            "peak_rss_mb": round(self.peak_rss / 2**20, 1),
        }


class _PeakSampler(threading.Thread):
    """Sample resident memory of this process and its children on a timer."""

    def __init__(self, interval: float = 0.1):
        super().__init__(daemon=True)
        self.interval = interval
        self.peak = 0
        self._stop = threading.Event()

    def run(self) -> None:
        proc = psutil.Process()
        while not self._stop.wait(self.interval):
            total = proc.memory_info().rss
            for child in proc.children(recursive=True):
                try:
                    total += child.memory_info().rss
                except psutil.NoSuchProcess:
                    pass
            self.peak = max(self.peak, total)

    def stop(self) -> int:
        self._stop.set()
        self.join(timeout=2)
        return self.peak


@contextmanager
def measure() -> Result:
    result = Result()
    sampler = _PeakSampler()
    sampler.start()
    started = time.perf_counter()
    try:
        yield result
    finally:
        result.wall = time.perf_counter() - started
        result.peak_rss = sampler.stop()


@contextmanager
def unit(result: Result):
    started = time.perf_counter()
    try:
        yield
    finally:
        result.durations.append(time.perf_counter() - started)

Using it around a real batch:

import json
from pathlib import Path

from bench.harness import measure, unit
from gistools.warp import warp_one


def run(paths: list[str], dst_dir: Path, dst_crs: str) -> dict:
    with measure() as result:
        for path in paths:
            with unit(result):
                warp_one(path, dst_dir / Path(path).name, dst_crs)
    summary = result.summary()
    print(json.dumps(summary, indent=2))
    return summary

Step Annotations

  1. time.perf_counter, never time.time. The wall clock can step backwards under NTP adjustment, which produces negative durations in a long batch and is impossible to explain later.
  2. The sampler includes child processes. A pool’s memory lives in the children, and sampling only the parent reports a flat line while the node runs out.
  3. _stop.wait(interval) rather than sleep. It makes the thread stop promptly at the end of the run instead of after one more full interval.
  4. Percentiles from a sorted list, not statistics.quantiles. The latter interpolates, which is correct statistically and confusing when someone tries to match a p95 against an actual unit in the log.
  5. The summary is JSON, not a printed table. A committed history of these objects is what makes when did this get slower answerable.

One Named Gotcha: The First Unit Is Always Slow

The first unit in any batch pays for imports, driver registration and the PROJ database open — often several hundred milliseconds that have nothing to do with the work. Included in the sample, it inflates the maximum and, in a small corpus, the p95 as well.

# Warm up before measuring, with one unit whose result is discarded.
warp_one(paths[0], dst_dir / "_warmup.tif", dst_crs)

with measure() as result:
    for path in paths:
        with unit(result):
            warp_one(path, dst_dir / Path(path).name, dst_crs)

Discarding the first unit rather than warming up separately also works, but is worse: it makes the unit count differ from the corpus size, and someone will eventually compare two runs whose counts differ without noticing.

The three numbers the harness reports move independently, and a change that improves one usually costs one of the others.

Every optimisation trades one number against another Adding workers raises throughput and peak memory together. Enlarging windows raises throughput and both peak memory and the tail. Adding a cache raises throughput on reuse-heavy work and always raises memory. more workers throughput up, peak memory up proportionally watch the limit larger windows throughput up, tail and memory both up watch the tail a larger block cache throughput up only if tiles are revisited measure the hit rate block-aligned reads throughput up, memory unchanged the free one

The last row is the only unambiguous win on the list, which is why access-pattern fixes are worth looking for before any of the others.

Verification

# Two runs of the same variant should agree closely.
python -m bench.run --variant baseline --json a.json
python -m bench.run --variant baseline --json b.json
python -c "
import json
a, b = json.load(open('a.json')), json.load(open('b.json'))
for k in ('throughput_per_min', 'p95_ms', 'peak_rss_mb'):
    d = abs(a[k] - b[k]) / max(a[k], 1) * 100
    print(f'{k}: {d:.1f}% apart')
"

Under five percent on all three means the harness is stable enough to compare variants. A large spread on peak_rss_mb specifically usually means the sampler interval is too coarse to catch the spike — halve it and re-measure.

Where a result is stored decides whether it can answer questions later.

Storing results so they stay useful Printed to a terminal, a result answers a question once. Committed per release, results build a history that shows when a regression entered. Pushed to a metrics system, results support alerting but lose the link to a specific commit. printed to a terminal answers one question, once gone when the shell scrolls not enough committed per release a history in git answers when did this get slower reviewable in a diff do this pushed to metrics supports alerting loses the link to a specific commit in addition

Committing per release rather than per commit keeps the history readable, because release-to-release changes are the ones anyone actually investigates.

Reporting Per Stage, Not Just Per Unit

A per-unit duration tells you a unit is slow; a per-stage breakdown tells you which part of it is. Adding stage timing costs one more context manager and changes what the result can answer.

@contextmanager
def stage(result: Result, name: str):
    started = time.perf_counter()
    try:
        yield
    finally:
        result.stages.setdefault(name, []).append(time.perf_counter() - started)

With read, warp and write timed separately, a summary can report the share of total time each takes. That share is what makes an optimisation decision obvious: a pipeline spending seventy percent in read is bound by I/O and will not improve from a faster transform, however satisfying that transform is to optimise.

The one caution is that stage times do not sum to the unit time in a concurrent pipeline, because stages overlap. Reporting them as shares of the sum rather than as shares of wall clock avoids a summary that appears to account for more than a hundred percent.

Comparing Two Variants Fairly

Once the harness produces stable numbers, comparing variants has three requirements that are easy to overlook.

Randomise the order. Running the baseline first and the candidate second, always, lets thermal throttling or a warming cache favour whichever ran second. Alternating between them across several repeats removes the bias.

Take the minimum, not the mean. For a benchmark, noise is almost always additive — something else on the machine took time away — so the fastest observed run is closest to the true cost. The mean measures the machine’s other tenants as much as your code.

Report the difference with its uncertainty. A candidate that is three percent faster with a harness that varies by five is not faster; saying so plainly is more useful than a number with a decimal point. The rule of thumb that works: only claim an improvement larger than twice the observed run-to-run spread.

Finally, keep the harness itself out of the measurement. Timing code that builds the corpus, resolves the work list or writes the result file inflates every number and does so unevenly, because those costs are fixed while the work is not. Starting the clock after setup and stopping it before teardown is obvious in principle and easy to get wrong when a fixture is created lazily inside the timed block — which is exactly what a generator-based work list does unless it is materialised first.

Where the harness fits

This measurement belongs beside the test suite rather than inside it. The fixtures are shared — both want a corpus with a realistic size distribution — but the invocation differs: tests run on every commit and must be fast, benchmarks run on demand and must be stable.

A bench/ package importing fixtures from tests/ gets the sharing without the coupling, and a make target keeps the invocation memorable. What does not work is a benchmark implemented as a pytest test with a timing assertion, because the threshold has to be loose enough to pass on a busy machine and is then too loose to catch anything.

Reporting to People Who Did Not Run It

A summary written for the person who ran the benchmark is usually unreadable to anyone else, and the audience that matters is the second one — a reviewer deciding whether to merge, or a colleague deciding whether to adopt a setting.

Three habits make the difference. State the question the measurement answers before the numbers, in one sentence, because a table without a question is just data. Give the comparison rather than the absolute: twelve percent faster than the previous release, on the same corpus and machine is actionable in a way that 4,200 tiles per minute is not. And name what got worse, if anything — peak memory, tail latency, request count — because a report with no trade-off in it invites the reader to look for the one you omitted.