Benchmarking Spatial Joins with pytest-benchmark

8 min read

A spatial join is the operation whose cost most often surprises a team, because it is the one whose complexity depends on something invisible: whether the index was used. This guide sits beneath performance benchmarking spatial suites and shows how to build a join benchmark with pytest-benchmark and GeoPandas 0.14+ that measures the right things — candidate counts as well as timings, several input sizes rather than one, and a growth exponent that fails when the join quietly becomes quadratic.

The measurement that matters is not how long the join took. It is how many candidate pairs the filter stage produced, because that number is exact, machine-independent, and it is the quantity that changes when an index stops being used.

Root cause: the timing and the complexity are different questions

A join that takes 40 ms today and 46 ms next month has got 15 per cent slower, which is a nuisance. A join that has gone from filtering candidates through an index to comparing every pair has changed complexity, and at fixture scale those two look identical.

The reason they are indistinguishable is arithmetic. At a thousand features the difference between O(nlog⁡n)O(n \log n) and O(n2)O(n^2) is a factor of about a hundred in the pair count but only a small factor in wall-clock time, because both are fast enough that fixed costs dominate. At a million features the same change is the difference between a query and an outage.

Timing hides at fixture scale; candidate count does not A comparison of two measurements taken at fixture scale for the same join run with and without a usable index. The wall-clock durations are close together, differing by only a small factor because fixed costs such as fixture loading and interpreter overhead dominate both. The candidate pair counts differ by two orders of magnitude, because the indexed join visits only spatial neighbours while the unindexed one visits every pair. A closing note records that the candidate count therefore detects the change at fixture scale, while the duration does not, and that the count is also machine-independent so it can be asserted rather than merely watched. Wall-clock at fixture scale — nearly identical indexed · 41 ms no index · 52 ms indexed no index a 1.3× difference — inside ordinary run-to-run variance Candidate pairs at the same scale — two orders apart 8 200 pairs 1 000 000 pairs indexed no index

Parameter reference

Setting Value that discriminates Why
Fixture Fixed seed, recorded config hash A changed fixture is a changed benchmark
Warm-up Discarded explicitly A cold import dominates the first round
Rounds Enough for a stable median The median is far more robust than the mean
Reported statistic Median and interquartile range The mean follows the tail
Companion metric Candidate pair count Machine-independent; the real signal
Sizes Four or five over an order of magnitude One size cannot show a slope
Worker count Fixed, usually one Concurrency is a variable, not a constant

Step-by-step implementation

The benchmark targets pytest-benchmark, GeoPandas 0.14+ and Shapely 2.x.

Step 1 — Build the fixtures once, at several sizes

The sizes are the whole point, and they must come from one generator so that only the count varies.

import pytest, geopandas as gpd
from shapely import STRtree

SIZES = [1_000, 4_000, 16_000, 64_000]

@pytest.fixture(scope="session", params=SIZES, ids=[f"n{n}" for n in SIZES])
def joined_pair(request):
    n = request.param
    left = build_points(n=n, seed=20260811)      # seeded — identical every run
    right = build_zones(n=max(64, n // 100), seed=20260811)
    return n, left, right

Step 2 — Separate the index build from the query

Benchmarking both together measures construction, which is a different question and usually the larger number.

def test_join_query_only(benchmark, joined_pair):
    n, left, right = joined_pair
    tree = STRtree(right.geometry.values)        # built outside the measurement
    def query():
        return [tree.query(g, predicate="within") for g in left.geometry.values]
    result = benchmark(query)
    assert sum(len(r) for r in result) > 0

Step 3 — Record the candidate count alongside the timing

This is the measurement that discriminates, and pytest-benchmark carries arbitrary extra info per benchmark.

def test_join_candidates(benchmark, joined_pair):
    n, left, right = joined_pair
    tree = STRtree(right.geometry.values)
    candidates = 0
    def query():
        nonlocal candidates
        candidates = sum(len(tree.query(g)) for g in left.geometry.values)
        return candidates
    benchmark(query)
    benchmark.extra_info["n"] = n
    benchmark.extra_info["candidates"] = candidates
    benchmark.extra_info["candidates_per_feature"] = candidates / n

candidates_per_feature is the number to watch. For a planar dataset with a working index it is roughly constant as nn grows; if it starts rising with nn, the index is not filtering.

Step 4 — Assert the exponent, not the duration

Fit the slope across the sizes and fail when it exceeds the bound. This is the one performance assertion worth gating on, because it is machine-independent.

import math

def test_join_complexity_is_subquadratic(benchmark_results):
    """benchmark_results: [(n, candidates), ...] collected across the sizes."""
    xs = [math.log(n) for n, _ in benchmark_results]
    ys = [math.log(c) for _, c in benchmark_results]
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    slope = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    assert slope <= 1.4, (
        f"candidate count grows as n^{slope:.2f} — the index is not filtering; "
        f"expected <= 1.4"
    )
Candidates per feature: flat means the index works A chart of candidate pairs per feature against input size, with four measured sizes. The working-index series is essentially flat across the whole range, because a feature in a planar dataset has approximately the same number of spatial neighbours whatever the total size of the dataset. The broken-index series rises linearly with the input size, because each feature is compared against every other one. Annotations identify the flat line as the signature of a functioning filter and the rising line as the signature of its absence, and note that the distinction is already unmistakable at fixture scale where a duration measurement would show almost nothing. candidates per feature input size n 1k4k16k64k index works no usable index Flat means the filter is doing its job; rising means every pair is being compared. Both are already unmistakable at the smallest fixture size.

Verify the fix

Break the index deliberately and confirm the exponent assertion fires:

pytest -q tests/bench/ --benchmark-only --benchmark-columns=median,iqr

Replace the tree query with a brute-force comparison and the slope assertion should fail with the measured exponent in its message. If instead only the timing changes and the slope assertion passes, the sizes span too narrow a range to fit anything.

Reading the output without being misled

Three habits keep a benchmark report honest, and all three are about which statistic is being read.

Read the median, not the mean. A single slow round caused by unrelated work on the runner moves the mean substantially and the median barely at all. Most harnesses report both; the median is the one that describes typical behaviour.

Read the interquartile range before the central value. A wide spread means the measurement is not stable enough for the central value to mean anything, and the right response is to fix the measurement rather than to record the number. On a shared runner a spread wider than the regression you are trying to detect is normal, which is exactly why the candidate count matters more than the duration.

Compare like against like. A comparison across different fixture sizes, different worker counts or different runner classes is measuring the difference in those rather than in the code. Record all three alongside every result, and refuse to compare when they differ.

There is a fourth habit worth adopting for spatial work specifically: check the result, not only the timing. A join that got dramatically faster has often started returning fewer rows — a predicate that no longer matches, a filter applied earlier than intended — and a benchmark that asserts only on duration reports that as an improvement. Asserting the result count alongside the timing turns a suspicious speed-up into a failure.

Benchmarking the database side

When the join runs in PostGIS rather than in memory, the same principles apply and the instruments change. The candidate count becomes the number of rows the planner examined, which the database reports directly, and the exponent is fitted over that instead of over an in-process counter.

Two database-specific measurements are worth adding. The plan shape is exact and machine-independent: capturing it and asserting that the spatial index still appears is a stronger check than any timing, because a plan change is the regression rather than a symptom of one. And the buffer statistics distinguish a query that is slow because it read a great deal from one that is slow because it computed a great deal, which points at completely different fixes.

The practical arrangement is a benchmark that runs the query with statistics collection enabled, records the plan alongside the rows examined, and asserts on both the exponent and the presence of the expected index scan. Timing is recorded and watched rather than asserted, exactly as it is in process.

Four fields that make a benchmark number comparable Four properties that must be recorded alongside every benchmark measurement. The fixture seed together with the configuration hash fixes the input, so a regenerated fixture is not mistaken for a code change. The engine versions and the container image digest fix the environment, so a runtime upgrade is attributable. The worker count fixes the concurrency, since a measurement taken with several workers and one taken serially measure different things. The warm-up policy fixes what the measurement includes, because a mean that contains a cold import is dominated by it. A closing note states that a number recorded without these four cannot be compared with any other number. fixture seed + config hash fixes the input — a regenerated fixture is not a code change engine versions + image digest fixes the environment — a runtime upgrade becomes attributable worker count fixes the concurrency — serial and parallel are different measurements warm-up policy fixes what is included — a cold import dominates a first round A number recorded without these four cannot be compared with any other number, because each of them moves the measurement independently of the code. Recording them costs four lines and is the difference between a series and a collection of unrelated figures. It is the same provenance discipline that makes a correctness failure attributable, applied to a different quantity.

Failure modes and edge cases

  1. Benchmarking construction and query together. Index construction usually dominates, so the measurement moves when construction changes and not when the query does. Build outside the measured callable.
  2. A fixture regenerated per run. Different data is a different benchmark. Seed it, record the seed and the config hash, and treat a change in either as a break in the series.
  3. Too few sizes, or too narrow a span. A slope fitted over three sizes within a factor of two is noise. Four or five over an order of magnitude is the minimum that says anything.
  4. Comparing medians across worker counts. Parallelism changes cache behaviour and memory pressure. Fix the worker count for the benchmark lane.
  5. A speed-up that is actually a correctness regression. Assert the row count alongside the timing, or a broken predicate reads as an optimisation.
  6. Gating the merge on a duration. It fails on a busy runner and teaches everyone to re-run. Gate on the exponent, alert on the duration.

Conclusion

A spatial join benchmark that measures only elapsed time cannot distinguish a modest slowdown from a complexity change, and the second is the one that matters. Fixing the fixture with a recorded seed, separating construction from query, recording candidates per feature alongside the timing, and asserting the fitted exponent gives a benchmark that catches a lost index at fixture scale — which is the whole purpose of performance benchmarking spatial suites.