Geospatial Orchestration Hub

Measuring Tile Generation Latency Percentiles

A latency percentile is only as good as the buckets underneath it. Prometheus histograms interpolate within a bucket, so a p95 computed from a histogram whose top bucket is ten seconds — while half the tiles take ninety — is not an estimate, it is arithmetic on a number that does not exist. Measure a real run first, choose bucket boundaries that bracket the observed distribution, and then read the result as a distribution rather than as three numbers, because tile latency is almost always bimodal and a percentile is exactly the wrong shape for that.

When to Use This Pattern

  • A latency dashboard exists and nobody believes it, usually because every observation is in the +Inf bucket.
  • Tile runtimes span orders of magnitude, which they do in any job covering both sea and mountains.
  • An alert on p95 fires constantly or never, both of which are symptoms of buckets that do not bracket the data.
  • You are about to set an SLO on tile latency and need the measurement to mean something first.

Complete Working Example

Start by measuring, without Prometheus in the loop at all. Twenty minutes of this decides every bucket boundary.

from __future__ import annotations

import statistics
import time
from collections import Counter as Tally
from pathlib import Path
from typing import Iterable


def measure_run(items: Iterable, render) -> list[float]:
    """Time every unit of one real run, and return the raw durations."""
    durations: list[float] = []
    for item in items:
        started = time.monotonic()
        render(item)
        durations.append(time.monotonic() - started)
    return durations


def summarise(durations: list[float]) -> dict[str, float]:
    ordered = sorted(durations)
    def q(p: float) -> float:
        return ordered[min(len(ordered) - 1, int(p * len(ordered)))]
    return {
        "n": len(ordered), "min": ordered[0], "p50": q(0.50), "p90": q(0.90),
        "p95": q(0.95), "p99": q(0.99), "max": ordered[-1],
        "mean": statistics.fmean(ordered),
    }


def suggest_buckets(durations: list[float]) -> tuple[float, ...]:
    """Roughly logarithmic boundaries bracketing the observed range.

    Two rules: the lowest bucket must sit below the observed minimum, and the
    highest FINITE bucket must sit above the observed p99. Between them,
    doubling gives enough resolution without producing thirty series per label.
    """
    lo = max(0.05, 2 ** (len(bin(int(min(durations)))) - 4))
    hi = max(durations)
    buckets, edge = [], lo
    while edge < hi * 1.5:
        buckets.append(round(edge, 2))
        edge *= 2
    return tuple(buckets) + (float("inf"),)


# A real run of 4 218 tiles:
#   {'n': 4218, 'min': 0.4, 'p50': 8.9, 'p90': 41.0, 'p95': 63.0,
#    'p99': 141.0, 'max': 218.0, 'mean': 19.4}
# suggest_buckets → (0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, inf)

With the boundaries chosen, the histogram declaration is unremarkable — and the labels stay bounded, because a per-tile label would defeat the whole exercise:

from prometheus_client import Histogram

TILE_SECONDS = Histogram(
    "tile_generation_seconds",
    "Wall-clock to generate one tile, from task start to output written",
    ("layer", "zoom", "resampling"),
    # Measured, not defaulted. The client library's defaults stop at 10 s, which
    # would place 52% of these tiles in +Inf and make every percentile fiction.
    buckets=(0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, float("inf")),
)
Default buckets against measured onesWith the library's default buckets, most tile durations fall into the overflow bucket and the distribution is invisible. With measured buckets, the distribution's two modes are clearly separated.DEFAULT BUCKETS — 0.005 s … 10 s+Inf — 52% of observationsAny percentile above p48 is interpolated inside an unbounded bucket. It is a number, not a measurement.MEASURED BUCKETS — 0.5 s … 256 s+Inf — 0.2%, as it should be
The lower row is the same data. The second mode around 64 seconds is entirely invisible in the upper one, and it is the mode that generates the incidents.

The bucket count deserves a moment’s thought, because more resolution is not free. Every boundary is a separate time series per label combination: with three labels producing 4 320 combinations and twelve buckets, the histogram alone is over fifty thousand series. Ten to twelve boundaries spanning three orders of magnitude gives roughly ±30% resolution near any quantile, which is enough to answer every question a latency dashboard is asked. Twenty boundaries doubles the storage and improves no decision anyone actually makes.

Boundaries against series count and resolutionSix boundaries give coarse resolution at low series count. Eleven give useful resolution at moderate cost. Twenty-two double the series for a resolution improvement nobody acts on.boundariesseriesquantile error626 000±70% — too coarse to act on1148 0002295 000±14%, and no decision changes22 boundaries
The middle row is the useful compromise, and it is where doubling boundaries across a measured range naturally lands.

Parameter & Option Reference

Choice Rule Spatial notes
lowest bucket below observed min Tiles over open sea finish in under a second; a lowest bucket of 1 s hides them all.
highest finite bucket above observed p99 With p99 at 141 s, 256 s is right and 64 s makes every high percentile fiction.
spacing doubling Ten to twelve boundaries. Finer spacing multiplies series without improving decisions.
labels 3, bounded layer, zoom, resampling. Each bucket is a series per label combination, so cardinality multiplies.
recalibration per worker class A change of instance type shifts the whole distribution; the buckets have to move with it.
native histograms consider Prometheus 2.40+ removes the bucket decision entirely, at the cost of a newer stack and less tooling support.

Verification & Testing

Two checks: that the buckets bracket the data, and that a query on them returns something close to the truth.

def test_buckets_bracket_the_observed_range(measured_durations) -> None:
    buckets = suggest_buckets(measured_durations)
    finite = [b for b in buckets if b != float("inf")]
    assert finite[0] < min(measured_durations), "the lowest bucket hides the fast tiles"
    p99 = sorted(measured_durations)[int(0.99 * len(measured_durations))]
    assert finite[-1] > p99, "the top finite bucket is below p99 — percentiles will be fiction"
    assert 8 <= len(finite) <= 14, "between eight and fourteen boundaries is the useful range"


def test_histogram_percentile_is_close_to_the_truth(measured_durations) -> None:
    for d in measured_durations:
        TILE_SECONDS.labels(layer="ortho", zoom="12", resampling="bilinear").observe(d)
    estimated = query_histogram_quantile(0.95, TILE_SECONDS)
    actual = sorted(measured_durations)[int(0.95 * len(measured_durations))]
    # Within a bucket width is the best a histogram can do, and it is enough.
    assert abs(estimated - actual) / actual < 0.35

Reading the result matters as much as recording it. A heatmap shows the whole distribution and makes bimodality obvious; three percentile lines average the two modes into a number that describes no tile that ever ran:

# For a heatmap panel: the bucket rates, which Grafana renders as a distribution.
sum by (le) (rate(tile_generation_seconds_bucket[5m]))

# For a single number, be explicit that it is per zoom — aggregating across
# zooms mixes populations whose latencies differ by an order of magnitude.
histogram_quantile(
  0.95,
  sum by (le, zoom) (rate(tile_generation_seconds_bucket{layer="ortho"}[10m]))
)
The percentile that describes nothingTile latencies form two modes, one around six seconds and one around seventy. The p95 falls between them, in a region where almost no tile actually lands.030 s80 s160 stile generation timep95 — in the empty valleyBlue is coastal and sea tiles; amber is dense terrain. They are two workloads, and one number describes neither.
The right response to this shape is usually to split the metric by a label that separates the modes — zoom, region, or a data-density band — rather than to keep quoting a percentile between them.

Common Pitfalls

  • Library default buckets. They are tuned for HTTP request latencies and top out around ten seconds. On raster work this puts most observations in +Inf, and every percentile above the overflow point is invented.
  • Aggregating percentiles across zooms. A zoom-8 tile and a zoom-16 tile are different populations. histogram_quantile over their union produces a number that belongs to neither, and it moves whenever the mix of zooms changes.
  • Averaging percentiles. avg(p95) across workers is not a percentile of anything. Aggregate the bucket rates with sum by (le, …) and compute the quantile once, at the end.
  • Setting an SLO before looking at the shape. A bimodal distribution cannot be summarised by a threshold on a percentile; the SLO will either exclude the entire second mode or permit the entire first one to degrade.
  • Never recalibrating. Buckets sized on a 4-vCPU worker are wrong on a 16-vCPU one. Recalibration is a five-minute job and belongs in the same change that resizes the workers.

Frequently Asked Questions

How accurate is `histogram_quantile` really?

Accurate to within a bucket width, because it interpolates linearly inside the bucket containing the quantile. With doubling boundaries that means roughly ±30% near the top of the range, which is fine for “is this getting worse?” and inadequate for “is this exactly 63 seconds?”. If you need the second, record the raw durations somewhere queryable and stop asking a histogram to be a database.

Should I use native histograms instead?

If your Prometheus and Grafana are recent enough, yes — they remove the bucket decision entirely and store an exponential representation with far better resolution. The caveats are ecosystem support, which still lags, and remote-write compatibility. For a new deployment they are worth evaluating; for an existing one, measured classic buckets are perfectly adequate.

What about summaries rather than histograms?

Summaries compute quantiles client-side, which means they cannot be aggregated across workers — the p95 of ten workers is not the average of their p95s, and a summary gives you no way to combine them. For a distributed pipeline that alone rules them out. Histograms aggregate correctly because the buckets are counters.

How do I find the label that separates a bimodal distribution?

Try the ones you already have: zoom, layer, resampling, source scheme. If none of them split it, the discriminator is probably data density, which is not a label you have — but tile-level pixel counts are, so bucketing tiles by megapixels and looking at the latency per bucket usually identifies it. That analysis belongs in a notebook rather than in a dashboard, and it only has to be done once.

Does the histogram include queue time?

Only if you start the timer when the task is scheduled rather than when it starts executing. Both are useful and they answer different questions: execution time tells you whether the work is slow, and queue-plus-execution tells you whether the pipeline is slow. Record both if the fan-out is large enough for queueing to matter, and label them clearly, because confusing the two produces alarming graphs during a backlog.

← Prometheus Metrics for Raster Throughput

Related in this section