Geospatial Orchestration Hub

Instrumenting gdalwarp with Prometheus Counters

gdalwarp is a subprocess, so it cannot emit metrics — but everything worth measuring about it is observable from the outside. Wrap the invocation, record what it was asked to do, how long it took, whether it succeeded, how many output pixels it produced and what its peak memory was, and you have throughput, latency, error rate and a memory profile without touching GDAL at all. The only trap is trying to parse its progress output, which is unstable across versions and tells you less than the boundary already does.

When to Use This Pattern

  • GDAL utilities run as subprocesses, which is the usual arrangement for anything that must be cancellable or memory-bounded.
  • Throughput needs to be attributable to a step, so “the pipeline is slow” can become “the warp is slow” or “the fetch is slow”.
  • Memory is the binding constraint, and you need evidence for the pool sizing rather than an estimate.
  • A GDAL upgrade is coming, and you want a before-and-after comparison that is not anecdotal.

Complete Working Example

The wrapper measures at the boundary and returns the metrics alongside the result, so nothing has to be re-derived later.

from __future__ import annotations

import os
import resource
import subprocess
import time
from dataclasses import dataclass
from typing import Sequence

from prometheus_client import Counter, Histogram

GDAL_RUNS = Counter(
    "gdal_invocations_total",
    "GDAL utility invocations",
    ("utility", "layer", "resampling", "outcome"),
)
GDAL_PIXELS = Counter(
    "gdal_output_pixels_total",
    "Output pixels written by GDAL",
    ("utility", "layer", "resampling"),
)
GDAL_SECONDS = Histogram(
    "gdal_duration_seconds",
    "Wall-clock of one GDAL invocation",
    ("utility", "layer", "resampling"),
    buckets=(1, 2, 5, 10, 20, 40, 80, 160, 320, float("inf")),
)
GDAL_PEAK_MB = Histogram(
    "gdal_peak_memory_mb",
    "Peak resident memory of the GDAL child process",
    ("utility", "layer"),
    buckets=(128, 256, 512, 1024, 2048, 4096, 8192, 16384, float("inf")),
)


@dataclass(frozen=True)
class GdalRun:
    outcome: str
    seconds: float
    peak_mb: float
    output_pixels: int


def run_and_measure(
    argv: Sequence[str], utility: str, layer: str, resampling: str, dest: str
) -> GdalRun:
    """Invoke a GDAL utility and record everything observable from outside it."""
    labels = {"utility": utility, "layer": layer, "resampling": resampling}

    # RUSAGE_CHILDREN accumulates, so take a before-reading and difference it.
    before = resource.getrusage(resource.RUSAGE_CHILDREN).ru_maxrss
    started = time.monotonic()

    proc = subprocess.run(argv, capture_output=True, start_new_session=True)

    seconds = time.monotonic() - started
    after = resource.getrusage(resource.RUSAGE_CHILDREN).ru_maxrss
    # ru_maxrss is the HIGH WATER MARK across all children, so this is a lower
    # bound when several run concurrently. Good enough to size a pool from.
    peak_mb = max(0.0, (after - before)) / 1024

    outcome = "ok" if proc.returncode == 0 else "failed"
    pixels = 0
    if outcome == "ok":
        pixels = output_pixel_count(dest)          # one header read, not a scan

    GDAL_RUNS.labels(outcome=outcome, **labels).inc()
    GDAL_SECONDS.labels(**labels).observe(seconds)
    GDAL_PEAK_MB.labels(utility=utility, layer=layer).observe(peak_mb)
    if pixels:
        GDAL_PIXELS.labels(**labels).inc(pixels)

    if outcome == "failed":
        raise RuntimeError(proc.stderr.decode("utf-8", "replace")[-500:])

    return GdalRun(outcome, seconds, peak_mb, pixels)


def output_pixel_count(path: str) -> int:
    """Pixels actually written — read from the output's header, not predicted."""
    import rasterio

    with rasterio.open(path) as src:
        return src.width * src.height * src.count

Calling it is unremarkable, which is the point — the instrumentation lives in one place and every GDAL step inherits it:

def warp(source: str, dest: str, layer: str, dst_crs: str = "EPSG:3857") -> GdalRun:
    return run_and_measure(
        [
            "gdalwarp", "-t_srs", dst_crs, "-r", "bilinear",
            "-dstnodata", "-9999",
            "-co", "COMPRESS=DEFLATE", "-co", "TILED=YES",
            "-wm", "512",                # bounded, so peak memory is meaningful
            source, dest,
        ],
        utility="gdalwarp", layer=layer, resampling="bilinear", dest=dest,
    )
Everything measurable at the boundaryBefore the call: the arguments, which give layer and resampling. During: wall-clock and peak resident memory. After: exit code and the output file's header, which gives the pixel count.beforeargv → labelsgdalwarp — opaquemeasurable from outside:wall-clockpeak RSS of the childexit codeafteroutput header → pixelsNothing here parses GDAL’s progress output, which changes between versions and adds no information the boundary lacks.
Four numbers from outside a process that offers no instrumentation of its own — enough for throughput, latency, error rate and pool sizing.

Reading the output’s header rather than predicting the pixel count is worth defending, because prediction is available and cheaper. The difference is what the metric means. A predicted count reports what the pipeline intended; a measured one reports what exists. Those diverge in exactly the cases worth knowing about — a warp clipped by a -te that was narrower than expected, a band dropped by a creation option, an output whose resolution was snapped to a grid. If the two are recorded side by side, their ratio becomes a cheap correctness check that costs one extra header read.

Predicted and actual, side by sideThree warps produce exactly the predicted pixel count. A fourth produces sixty per cent of it because the target extent clipped the output, which no other check would have noticed.warppredictedactualratio12/2100/13001.0012/2101/13001.0012/2102/13000.60 — clipped by -te12/2103/13001.00
The third row exited zero and produced a valid GeoTIFF. Only the ratio says anything is wrong, and it costs one header read to have.

Parameter & Option Reference

Item Source Spatial notes
utility label argv[0] gdalwarp, gdal_translate, gdaladdo. Bounded and useful for splitting a pipeline’s time.
resampling label -r value Six values. Correlates strongly with duration, so it belongs in the labels.
duration buckets measured 1 s … 320 s. The client library’s defaults stop at 10 s and would collapse every warp into +Inf.
memory buckets measured 128 MB … 16 GB. The distribution this reveals is what the pool size should be derived from.
ru_maxrss RUSAGE_CHILDREN A high-water mark across children — a lower bound under concurrency, and adequate for sizing.
output pixels output header Read after the fact rather than predicted, so the metric reports what was produced.
-wm 512 argv Bounds GDAL’s own memory so the peak-RSS distribution is stable rather than machine-dependent.

Verification & Testing

Assert that the wrapper records on both paths and that the pixel count matches the file.

def test_failure_is_recorded_then_raised(tmp_path) -> None:
    before = sample_counter(GDAL_RUNS, outcome="failed")
    with pytest.raises(RuntimeError):
        run_and_measure(["gdalwarp", "/nonexistent.tif", str(tmp_path / "o.tif")],
                        utility="gdalwarp", layer="ortho", resampling="bilinear",
                        dest=str(tmp_path / "o.tif"))
    assert sample_counter(GDAL_RUNS, outcome="failed") == before + 1


def test_pixels_match_the_output(tmp_path, sample_cog) -> None:
    dest = tmp_path / "warped.tif"
    run = warp(sample_cog, str(dest), layer="ortho")
    with rasterio.open(dest) as src:
        assert run.output_pixels == src.width * src.height * src.count


def test_peak_memory_is_plausible(tmp_path, large_cog) -> None:
    run = warp(large_cog, str(tmp_path / "big.tif"), layer="ortho")
    # A warp with -wm 512 should never report a peak below its own cache budget.
    assert 400 < run.peak_mb < 8000

Once the series exist, three PromQL expressions cover the questions people actually ask during an incident. The third is the one worth having pinned, because it separates a pipeline that is slow from one that is failing:

# Throughput, in megapixels per second, by layer.
sum by (layer) (rate(gdal_output_pixels_total[5m])) / 1e6

# p95 warp duration, by resampling — where the tail actually lives.
histogram_quantile(0.95, sum by (le, resampling) (rate(gdal_duration_seconds_bucket[10m])))

# Failure ratio. Rising while throughput holds means retries are absorbing it.
sum(rate(gdal_invocations_total{outcome="failed"}[10m]))
  / sum(rate(gdal_invocations_total[10m]))
Peak memory is a distribution, not a numberMost warps peak below one gigabyte, but a long tail reaches six. A pool sized from the mean would be overcommitted by a factor of four against the tail.02 GB4 GB6 GBpeak resident memory of one gdalwarpmean — misleadingp99 — size the pool here
Without this histogram, pool sizing is a guess. With it, the number is a lookup: usable memory divided by the p99, which is exactly the calculation the process-pool recipe uses.

Common Pitfalls

  • Parsing GDAL’s progress output. The percentage line is meant for humans, changes format between versions, and interleaves unhelpfully when several processes share a terminal. Everything it would tell you is available from the boundary.
  • Recording only successes. A counter incremented after check=True misses every failure, so the error rate reads as zero exactly when it is not. Record in a finally, or before raising.
  • Predicting the pixel count instead of reading it. The prediction is what the pipeline asked for; the header is what GDAL wrote. When they differ — a clipped extent, a dropped band — only the second one is a metric.
  • Using ru_maxrss without differencing. It accumulates across all children of the process, so an absolute reading grows monotonically and reports the largest warp of the day for every tile. Difference it around the call.
  • Omitting -wm. Without a warp-memory bound, GDAL sizes its buffers from available memory, so the peak-RSS distribution reflects how busy the machine was rather than how large the tile was. The metric becomes unusable for sizing.

Frequently Asked Questions

Is `ru_maxrss` accurate enough under concurrency?

It is a lower bound, because it reports the high-water mark across all children rather than per child. With a process pool running four warps, the difference around one call may attribute another warp’s peak to this one, or miss its own. For pool sizing that is acceptable — you want the upper tail, and the bias is in the safe direction. For per-tile attribution, sample /proc/<pid>/status from a monitoring thread instead.

Should the wrapper also record the source's size?

As a label, no — sizes are unbounded. As a separate observation, yes: a histogram of input megapixels alongside the output one reveals whether a slow run was a big job or a slow machine. Two histograms are cheap and they answer a question that neither answers alone.

How do I compare before and after a GDAL upgrade?

Keep a gdal_version label for the duration of the rollout only, then remove it. It is bounded (two values during a migration), it makes the comparison a single query, and leaving it in permanently accumulates a series per version forever. This pairs with the rollout process in rolling out GDAL upgrades without breaking flows.

What about `gdal_translate` and `gdaladdo`?

Same wrapper, different utility label. Their duration distributions are quite different — gdaladdo is dominated by reading the finished file — so give them their own buckets if the shared ones put most observations in one bucket. The rest of the instrumentation is identical, which is the argument for one wrapper rather than three.

Does this add measurable overhead?

No. Four metric observations and one header read per invocation, against a subprocess that runs for seconds to minutes. The header read is the largest of them and is typically under a millisecond on a local file. Instrumentation cost is not a reason to skip this.

← Prometheus Metrics for Raster Throughput

Related in this section