Benchmarking H3 Lookup Against an R-tree Query

A cell lookup is a dictionary access and a tree query is a traversal, so the benchmark appears settled before it starts. It is not, because the two indexes answer different questions and the work they push downstream differs by more than the lookup times differ.

Prerequisites

Conceptual minimum

The two indexes are not interchangeable, and a benchmark that ignores that measures the wrong thing.

An R-tree answers “which geometries’ bounding boxes intersect this region”. The answer is a candidate set that must then be refined by an exact geometric test, because a bounding box overlap is not containment. The refinement is usually the larger cost.

An H3 index answers “which features are in these cells”. The answer is exact for the cells, and approximate for the region, because a set of hexagons only approximates any shape that is not made of hexagons. Refinement is needed at the boundary, and how much depends on the resolution chosen.

A fair comparison therefore measures total time to a correct answer, including the refinement each approach requires, and reports candidates examined alongside, because that number explains the timing rather than merely restating it.

The second fairness issue is the query shape. H3 is at its best when the query is itself a cell or a set of cells — an aggregation, a join on a grid. An R-tree is at its best for an arbitrary polygon or a nearest-neighbour search. Benchmarking only one shape answers only for that shape.

What each index actually answers, and what has to happen afterwards A grid of four properties against the two indexes. An R-tree returns geometries whose bounding boxes intersect the query region, requires an exact geometric refinement afterwards, handles arbitrary query shapes, and supports nearest-neighbour search directly. A cell index returns features assigned to the queried cells, requires refinement only near the region boundary, handles cell-shaped queries natively and arbitrary shapes only through approximation, and supports nearest-neighbour search poorly. Two indexes, two different questions answered R-tree Cell index Returns box intersections cell members Refinement always needed boundary only Arbitrary shapes native approximated Nearest neighbour direct poorly A benchmark that runs only cell-shaped queries measures the second column at its best and the first at its worst.
Choosing the query shape is choosing the winner, which is why the workload has to come from the real application.

Runnable solution

python 
from __future__ import annotations

import logging
import time
from collections import defaultdict
from dataclasses import dataclass

import h3
import numpy as np
from rtree import index as rtree_index
from shapely.geometry import Point, Polygon, shape
from shapely.prepared import prep

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger("osm.index.benchmark")


@dataclass(frozen=True)
class Timing:
    label: str
    build_seconds: float
    query_seconds: float
    refine_seconds: float
    candidates: int
    results: int

    @property
    def total_seconds(self) -> float:
        return self.query_seconds + self.refine_seconds


def build_rtree(points: np.ndarray) -> tuple[rtree_index.Index, float]:
    started = time.perf_counter()
    def stream():
        for i, (x, y) in enumerate(points):
            yield (i, (x, y, x, y), None)
    tree = rtree_index.Index(stream())
    return tree, time.perf_counter() - started


def build_h3(lonlat: np.ndarray, resolution: int) -> tuple[dict, float]:
    started = time.perf_counter()
    cells: dict[str, list[int]] = defaultdict(list)
    for i, (lon, lat) in enumerate(lonlat):
        cells[h3.latlng_to_cell(lat, lon, resolution)].append(i)
    return dict(cells), time.perf_counter() - started


def query_rtree(tree, points: np.ndarray, polygons: list[Polygon]) -> Timing:
    query_seconds = refine_seconds = 0.0
    candidates = results = 0
    for polygon in polygons:
        started = time.perf_counter()
        hits = list(tree.intersection(polygon.bounds))
        query_seconds += time.perf_counter() - started
        candidates += len(hits)

        # Refinement is NOT optional: a bounding box hit is not containment.
        started = time.perf_counter()
        ready = prep(polygon)
        results += sum(1 for i in hits
                       if ready.contains(Point(points[i][0], points[i][1])))
        refine_seconds += time.perf_counter() - started
    return Timing("r-tree", 0.0, query_seconds, refine_seconds, candidates, results)


def query_h3(cells: dict, lonlat: np.ndarray, polygons: list[Polygon],
             resolution: int) -> Timing:
    query_seconds = refine_seconds = 0.0
    candidates = results = 0
    for polygon in polygons:
        started = time.perf_counter()
        # Cover the polygon with cells; this is the approximation step.
        covering = h3.geo_to_cells(polygon.__geo_interface__, resolution)
        hits: list[int] = []
        for cell in covering:
            hits.extend(cells.get(cell, ()))
        query_seconds += time.perf_counter() - started
        candidates += len(hits)

        # Cells that lie wholly inside need no refinement; boundary cells do.
        started = time.perf_counter()
        ready = prep(polygon)
        results += sum(1 for i in hits
                       if ready.contains(Point(lonlat[i][0], lonlat[i][1])))
        refine_seconds += time.perf_counter() - started
    return Timing(f"h3 r{resolution}", 0.0, query_seconds, refine_seconds,
                  candidates, results)


def report(timings: list[Timing]) -> None:
    for t in timings:
        logger.info("%-10s query %6.3fs + refine %6.3fs = %6.3fs  "
                    "%8d candidate(s) -> %7d result(s)  (%.1fx selectivity)",
                    t.label, t.query_seconds, t.refine_seconds, t.total_seconds,
                    t.candidates, t.results,
                    t.candidates / max(1, t.results))
    # Identical results are the precondition for comparing anything else.
    distinct = {t.results for t in timings}
    if len(distinct) > 1:
        logger.error("indexes returned different result counts %s — the "
                     "comparison is meaningless until they agree", distinct)


if __name__ == "__main__":
    logger.info("run with a real workload; random boxes favour whichever "
                "index the box shape happens to suit")

Step-by-step walkthrough

  1. Time the refinement separately. An index query that returns quickly and hands back ten times as many candidates has not won; separating the two numbers is what makes that visible.
  2. Use a prepared geometry for the refinement. Preparing the polygon once per query rather than per candidate removes a cost that would otherwise dominate and distort the comparison.
  3. Count candidates alongside time. The candidate count is the structural measurement and it explains the timing; without it, a result is a number with no mechanism behind it.
  4. Report selectivity. Candidates divided by results says how much work each index wasted, which is the number that transfers to a different machine.
  5. Assert the results agree. Two indexes returning different counts are answering different questions, and every other comparison is meaningless until that is resolved.
  6. Cover the polygon at query time for the cell index. The covering step is part of the query cost and omitting it from the timing flatters the cell index substantially.
  7. Use a real workload. Randomly generated query boxes favour whichever index the box shape suits, which is a property of the benchmark rather than of the application.
Where the time goes for each index on a polygon containment workload Four measurements over the same set of polygon containment queries against a million OSM points. The R-tree's index query is fast but returns many candidates because bounding boxes overlap the polygon generously. Its refinement dominates its total. The cell index's query is slower because covering the polygon with cells is real work, but it returns far fewer candidates for a compact polygon. Its refinement is correspondingly smaller, and the totals end up close. Query and refinement, for a polygon containment workload R-tree: index query fast R-tree: refinement dominates Cell index: covering real work Cell index: refinement fewer candidates Totals are close here and diverge sharply with query shape: a long thin polygon reverses the candidate advantage entirely.
Comparing only the first and third bars — the index query alone — would report a result the totals contradict.
Three ways a spatial index benchmark reaches a conclusion that does not transfer Three panels. Omitting the refinement measures only the index lookup, which flatters whichever index returns candidates fastest regardless of how many it returns. Using randomly generated query boxes measures a workload nobody has, and the box shape decides the winner before any code runs. Testing one data distribution hides that a clustered dataset and a dispersed one favour different structures, so the conclusion reverses on the next region. Three benchmark flaws, three wrong conclusions Refinement omitted Only the lookup timed Candidate count ignored Fastest lookup wins Total time disagrees Random query boxes A workload nobody has Shape decides the winner Before any code runs Does not transfer One distribution Clustered or dispersed Not both Conclusion reverses On the next region All three produce a confident number, which is why a benchmark result should always be reported with the workload it used.
The fix for all three is the same: measure the application's real queries on the application's real data.

Verification

  • Result counts match exactly. Any difference means the two are not answering the same question.
  • Selectivity is reported. A ratio near one means the index is doing nearly all the work; a large ratio means the refinement is.
  • Query shape is varied. Run compact and elongated polygons; the candidate advantage should move between the indexes.
  • The cell resolution is swept. Too coarse and the covering is inaccurate; too fine and the covering itself becomes expensive.
  • Build time is reported but separated. It matters for a one-shot job and not at all for a long-lived index.

Common errors and fixes

Symptom Root cause One-line fix
Cell index looks unbeatable Covering excluded from the query timing Time the covering as part of the query
R-tree looks unbeatable Refinement excluded from the timing Time the exact test as part of the answer
Results differ between indexes Different predicates applied Assert equal result counts before comparing
Refinement dominates everything Polygon prepared per candidate Prepare the geometry once per query
Benchmark does not transfer Random query boxes used Benchmark the application’s real query shapes
Cell index slow at fine resolution Covering cost grows with cell count Sweep resolutions and report the curve
Conclusions reverse on other data Single distribution tested Test both clustered and dispersed data

Specification reference

H3 assigns each location to a hexagonal cell at a chosen resolution, so a containment query is a set membership test over the cells covering the query region, with approximation error at the boundary. An R-tree returns entries whose bounding rectangles intersect the query rectangle, which is a necessary but not sufficient condition for geometric containment, so an exact test must follow. See the H3 documentation for the covering operations and Spatial Index Selection: R-tree vs H3 vs Quadkey for the wider comparison.

Frequently Asked Questions

Why is comparing index lookup times alone misleading?

Because neither index returns a final answer. An R-tree returns bounding-box candidates that must be tested exactly, and a cell index returns cell members that must be tested near the region boundary. The refinement is frequently the larger cost, so a benchmark that stops at the index lookup can report a winner that the total time contradicts. Measuring time to a correct answer is the only comparison that transfers.

Does the query shape really change the outcome?

Substantially. A compact polygon is covered efficiently by hexagons and matches an R-tree’s bounding box loosely, favouring the cell index. A long thin diagonal polygon is covered by many hexagons of which most are nearly empty, and its bounding box is enormous relative to its area, which hurts both — but differently. Benchmarking one shape and generalising is the most common way these comparisons mislead.

What cell resolution should the benchmark use?

Several, reported as a curve. Too coarse and the covering approximates the query region badly, inflating candidates. Too fine and the covering itself becomes expensive, because the number of cells grows quickly with resolution. The best resolution depends on the size of the query regions relative to the cells, so a single value answers only for the region size it was chosen against.

Should build time be included in the comparison?

Reported separately, and weighted by how the index is used. For a long-lived index queried millions of times, build cost is irrelevant and query cost is everything. For a one-shot join where the index is built and discarded, build cost can dominate and a cell assignment — a single arithmetic operation per feature — is hard to beat. Stating which regime applies is part of the result.

Up one level: Spatial Index Selection: R-tree vs H3 vs Quadkey.