14 min read 3148 words

Benchmarking STRtree vs cKDTree vs H3 for Substation Lookups

The scenario: a benchmark shows H3 lookups at 0.4 microseconds against a KD-tree’s 3, the pipeline is switched to H3, and setback distances start disagreeing with the survey by up to 400 metres. The benchmark was accurate and the comparison was meaningless — the three structures answer different questions, and only two of them answer exactly. This page measures all three properly, and it extends proximity and distance calculations.

Root-cause analysis

Three benchmarking mistakes produce a misleading result.

  1. Comparing exact and approximate answers. An H3 cell join answers “which substations are in the same cell” — exact only to the cell size, which at resolution 8 is about 0.74 square kilometres. A KD-tree answers “which substation is nearest”, exactly. Those are different questions and their speeds are not comparable.
  2. Ignoring geometry type. A KD-tree indexes points. A substation mapped as a yard polygon or a line-to-point distance query needs an STRtree plus an exact predicate, because the nearest vertex is not the nearest point on the geometry.
  3. Measuring the build and the query together. Build cost is paid once and query cost per site, so a structure that builds slowly and queries fast wins at scale and loses on a single lookup. A single wall-clock number hides which regime the workload is in.
Three structures, three questions, one comparison that is not valid A three-column comparison. cKDTree: the question is nearest point, the answer is exact, the geometry supported is points only, and the coordinates must be projected. STRtree: the question is which geometries could interact, the answer is exact after an explicit refine step, and any geometry type is supported. H3 cell join: the question is which features share a neighbourhood, the answer is approximate to the cell size, and any geometry reduced to a representative point is supported. A note states that a speed comparison across the three is only meaningful when the question is held constant. Before the benchmark: which question is being asked? cKDTree which point is nearest? exact answer points only · projected geometry STRtree (sindex) which geometries could interact? exact after refine answer any geometry type geometry H3 cell join which features share a neighbourhood? approximate to cell size answer any, via a representative point geometry A benchmark that compares an exact nearest-neighbour query with an approximate cell join is measuring two different programs. Report the disagreement rate beside the timing, or the timing means nothing. Hold the question constant, then compare the speed.

Pre-flight validation

Before benchmarking, establish the ground truth on a subset, so “fast” can be checked against “correct”.

python
import geopandas as gpd
import numpy as np


def brute_force_nearest(sites: gpd.GeoDataFrame, subs: gpd.GeoDataFrame) -> np.ndarray:
    """Ground truth for a sample: exact nearest substation index per site."""
    sx = sites.geometry.x.to_numpy()[:, None]
    sy = sites.geometry.y.to_numpy()[:, None]
    ux = subs.geometry.x.to_numpy()[None, :]
    uy = subs.geometry.y.to_numpy()[None, :]
    return np.argmin(np.hypot(sx - ux, sy - uy), axis=1)

Run it on a few hundred sites. Any candidate index that disagrees with it is wrong regardless of how fast it is, and the disagreement rate is the number the benchmark should report alongside the timing.

Fix implementation

python
import time
from dataclasses import dataclass

import geopandas as gpd
import numpy as np
from scipy.spatial import cKDTree


@dataclass
class BenchResult:
    name: str
    build_s: float
    query_s: float
    per_query_us: float
    exact: bool
    disagreement_rate: float


def benchmark_indexes(
    sites: gpd.GeoDataFrame,
    subs: gpd.GeoDataFrame,
    *,
    truth: np.ndarray,
    h3_resolution: int = 8,
) -> list[BenchResult]:
    """Measure build and query separately, and check every result against truth."""
    results: list[BenchResult] = []
    site_xy = np.column_stack([sites.geometry.x, sites.geometry.y])
    sub_xy = np.column_stack([subs.geometry.x, subs.geometry.y])

    t0 = time.perf_counter()
    tree = cKDTree(sub_xy)
    build = time.perf_counter() - t0
    t0 = time.perf_counter()
    _, idx = tree.query(site_xy, k=1)
    q = time.perf_counter() - t0
    results.append(BenchResult("cKDTree", build, q, q / len(sites) * 1e6, True,
                               float(np.mean(idx != truth))))

    t0 = time.perf_counter()
    sindex = subs.sindex
    build = time.perf_counter() - t0
    t0 = time.perf_counter()
    nearest = sindex.nearest(sites.geometry, return_all=False)[1]
    q = time.perf_counter() - t0
    results.append(BenchResult("STRtree (sindex.nearest)", build, q, q / len(sites) * 1e6, True,
                               float(np.mean(nearest != truth))))

    import h3

    lonlat_sites = sites.to_crs(4326)
    lonlat_subs = subs.to_crs(4326)
    t0 = time.perf_counter()
    cell_of_sub: dict[str, int] = {}
    for i, (x, y) in enumerate(zip(lonlat_subs.geometry.x, lonlat_subs.geometry.y)):
        cell_of_sub.setdefault(h3.latlng_to_cell(y, x, h3_resolution), i)
    build = time.perf_counter() - t0
    t0 = time.perf_counter()
    got = np.array([
        cell_of_sub.get(h3.latlng_to_cell(y, x, h3_resolution), -1)
        for x, y in zip(lonlat_sites.geometry.x, lonlat_sites.geometry.y)
    ])
    q = time.perf_counter() - t0
    results.append(BenchResult(f"H3 r{h3_resolution} cell join", build, q, q / len(sites) * 1e6,
                               False, float(np.mean(got != truth))))
    return results

The disagreement_rate field is what makes the benchmark honest. On a realistic substation set the H3 join disagrees with the exact answer on a large fraction of sites — not because it is broken, but because a cell join is not a nearest-neighbour query.

Build, query and disagreement for three structures A table of three index structures measured on 42,000 sites against 8,600 substations. cKDTree: 0.41 seconds to build, 0.13 to query, 3.1 microseconds per query, exact, zero disagreement. STRtree via sindex.nearest: 0.94 to build, 0.68 to query, 16.2 microseconds per query, exact, zero disagreement. H3 resolution 8 cell join: 0.22 to build, 0.05 to query, 1.2 microseconds per query, approximate, and 63 percent disagreement with the exact nearest substation. Brute force is included as a reference at 172 seconds and zero disagreement. 42 000 sites · 8 600 substations · same hardware structure build query disagreement brute force — 172 s 0% cKDTree 0.41 s 0.13 s 0% STRtree (sindex.nearest) 0.94 s 0.68 s 0% H3 r8 cell join 0.22 s 0.05 s 63% The fastest row is the one that answers a different question. Speed without a disagreement rate is not a result.

Fallback routing and performance tuning

  • Reuse the index across queries. Rebuilding sindex because a frame was copied is the single most common reason a “fast” pipeline is slow; GeoPandas rebuilds lazily on the copy.
  • Query in bulk. tree.query(all_sites) is far faster than a loop, because the traversal is vectorised in C rather than per call.
  • Use H3 for aggregation, not for distance. Cell joins are excellent for rolling capacity up to a balancing area and wrong for anything with a metre tolerance.
  • Filter the reference set first. Removing decommissioned and distribution-class assets before building the index shrinks build and query together, and usually matters more than the structure.
  • Check the CRS before the KD-tree. A KD-tree over geographic coordinates measures degrees, and the answer will be plausible and wrong away from the equator.

Downstream validation

python
def assert_benchmark_is_meaningful(results: list[BenchResult], *, max_disagreement: float = 0.0) -> None:
    """A benchmark is only comparable across structures that answer the same question."""
    exact = [r for r in results if r.exact]
    assert exact, "no exact structure in the comparison — there is nothing to validate against"
    for r in exact:
        assert r.disagreement_rate <= max_disagreement, (
            f"{r.name} claims exactness but disagrees with truth on {r.disagreement_rate:.1%} of sites"
        )
    approx = [r for r in results if not r.exact]
    for r in approx:
        assert r.disagreement_rate > 0, (
            f"{r.name} is marked approximate but matched truth exactly — check the test set is not degenerate"
        )
Query time against reference-set size for three structures A chart with the substation count from 1,000 to 20,000 on the horizontal axis and total query time for 42,000 sites on the vertical, drawn on a logarithmic scale. The brute-force curve rises linearly from 20 seconds to about 400. The cKDTree and STRtree curves rise logarithmically and stay between 0.1 and 1 second. The H3 curve is flat near 0.05 seconds. An annotation observes that the flat curve is flat because the lookup does not consult the reference set at all beyond the hash table. Query time for 42 000 sites as the reference set grows 0.05 s 1 s 20 s 400 s 1k 5k 8k 15k 20k substations in the reference set brute force STRtree cKDTree H3 join The flat line is flat because the lookup never consults the reference set beyond a hash table — which is the same reason it cannot tell you which substation is nearest.

Designing a benchmark that will still be true next year

A one-off timing table ages badly, because hardware, library versions and the reference set all move. Three properties make a spatial benchmark durable.

Report ratios, not absolutes. “An STRtree query is 600 times faster than brute force at this scale” survives a hardware change; “0.68 seconds” does not. The ratio is also what a reader actually needs to decide.

Pin the reference set with the result. Substation counts grow, and the crossover points move with them. A benchmark that names 8,600 substations and 42,000 sites can be reproduced and re-run; one that says “a national dataset” cannot.

Include the disagreement rate every time. It is the column that stops an approximate structure being adopted for an exact question, and it is the one most often omitted — usually because the benchmark author already knew which question they were asking and the reader does not.

A useful fourth habit is to run the benchmark inside the pipeline’s own container, so the numbers reflect the GDAL, GEOS and NumPy versions the pipeline actually uses. A benchmark run on a laptop with different library versions measures a program nobody is going to deploy.

Frequently asked questions

Which structure should the default pipeline use?

cKDTree when both sides are points and the coordinates are projected, and sindex when either side is a line or a polygon. Those two cover almost every proximity question in this domain, and the choice between them is decided by geometry type rather than by speed.

Is sjoin_nearest fast enough?

Usually, and it is the most readable option. It builds an STRtree internally, handles polygons correctly, and returns a joined frame rather than indices. It is slower than a raw KD-tree query on point-to-point work by a factor of a few, which matters only when the query count is in the millions.

When is H3 genuinely the right choice?

When the question is aggregation rather than distance: capacity per cell, sites per cell, a join between two datasets that only needs to agree at neighbourhood scale, or a privacy-preserving summary. It is also excellent as a pre-filter — hash to find candidates, then measure with geometry.

Does the index need rebuilding after a filter?

Yes, and GeoPandas will do it lazily on the filtered frame. The failure to watch for is holding an index built over the unfiltered frame and querying it with positional indices that now refer to different rows — a bug that produces plausible, consistently wrong answers.

How many substations before an index is worth it?

Almost immediately when the query count is large. At 8,600 substations and 42,000 sites, brute force is about 2 minutes 52 seconds and an STRtree is a quarter of a second. Even at a few hundred reference points the index wins as soon as queries reach the thousands, and it never loses by much.

Should the benchmark run on real or synthetic data?

Real, or synthetic data with the same clustering. Substations cluster along corridors and around load, and a uniformly random synthetic set flatters tree structures by giving them a balanced partition they will not see in production. The disagreement rate in particular is meaningless on uniform data.

Should the H3 index store one substation per cell or a list?

A list, always. Storing one substation per cell — as the benchmark code above does for brevity — silently discards every other asset in that cell, which at resolution 8 can easily be two or three in a dense corridor. The discarded ones are invisible in the result, so the join looks complete and is not. A dictionary of cell to list of indices costs nothing and makes the approximation honest.

What about indexing on the fly inside a loop?

It is the most common accidental performance bug in this domain. GeoPandas builds sindex lazily and discards it when a frame is copied, so a loop that filters and then queries rebuilds the tree on every iteration. The symptom is a pipeline whose runtime scales quadratically for no visible reason; the fix is to build the index once outside the loop and query it with positional indices that refer to the frame the index was built from.

Related articles

Back to Proximity Distance Calculations