R-tree vs H3 vs Quadkey: Spatial Index Selection Jump to heading#

The spatial index you pick for an OpenStreetMap workload is a decision made once and paid on every query for the life of the pipeline, and the wrong choice fails quietly rather than loudly. A team indexing a national point-of-interest layer into H3 resolution 8 to answer “which shops fall inside this administrative boundary” ships a plausible-looking service that is subtly wrong at every border: hexagons approximate the polygon, straddle its edge, and count cafés on the far side of a river as inside the district. Another team reaches for an R-tree to drive a country-wide density heatmap and discovers the tree gives them fast candidate lookups but no natural cell to aggregate into, so they bolt on an ad-hoc grid and re-derive it inconsistently across reports. Neither failure is a bug in the library — both are a mismatch between the index’s geometry and the question being asked. This page is the selection guide for that decision, sitting inside the broader OSM Data Fundamentals & Architecture layer, and it is the companion to the how-to-build reference: where Spatial Indexing for OSM Extracts shows how to construct a disk-backed R-tree, this page is the decision matrix it points to for which index to construct in the first place.

Start here: the three questions that decide it Jump to heading#

Before comparing families, resolve three questions about the workload, because they eliminate options faster than any benchmark. First, is the query exact or approximate? A topology join, a nearest-road snap, or a point-in-polygon test that must be correct at the boundary needs the true geometry — only an R-tree gives you a candidate set you can refine to an exact answer. A “roughly how many features per area” question tolerates cell approximation and opens the door to grids. Second, does the output have to align to something external? If results feed a slippy-map renderer, a tile server, or a Web Mercator pyramid, the cell scheme is already chosen for you — it is the quadkey. Third, must every cell cover comparable ground area? Density statistics, choropleths, and machine-learning features over a wide latitude span want near-uniform cells, which is exactly what planar grids fail to give and hexagonal H3 is built to provide.

These questions map onto three families that are complementary rather than rival: bounding-box hierarchies (R-tree), planar tile grids (Quadkey), and global discrete grids (H3 and S2). Most mature OSM stacks end up running two of them for different questions, so the goal is not to crown a winner but to route each query to the structure whose geometry matches it.

Prerequisites Jump to heading#

This guide assumes you already understand what you are indexing. Read the Node-Way-Relation Data Model first, because whether a primitive is a point, a line, or a polygon changes which index even applies — you can drop a node straight into an H3 cell, but a way is only indexable once its node references resolve into a geometry. Second, the Coordinate Reference Systems in OSM reference matters here more than anywhere: OSM stores unprojected WGS 84, quadkeys live in Web Mercator, and H3 works on the sphere, so every family answers the latitude-distortion question differently. Third, keep the Tag Taxonomy & Key-Value Standards conventions in view, since the smart move is often per-feature-class indexing — a separate structure for highway=* versus amenity=* rather than one monolithic index.

The decision matrix Jump to heading#

The table below is the fast path. Read down the column for your dominant query and the winning family is usually obvious; the sections after it explain the trade-offs when two columns tie.

Two rows deserve emphasis. Area uniformity is the single reason to leave a quadkey grid behind: because a Web Mercator cell’s ground area scales with $\cos(\varphi)$, a fixed-zoom grid over a country spanning many degrees of latitude gives Oslo cells a fraction of the ground area of Naples cells, silently biasing any per-cell statistic. Point-in-polygon prefilter is the reason to keep an R-tree even in a grid-heavy stack: only bounding-box candidate resolution followed by an exact predicate returns geometrically correct membership at a boundary.

R-tree: when exact geometry wins Jump to heading#

Choose an R-tree whenever the answer must be geometrically correct — spatial joins between two OSM layers, snapping GPS traces to roads, clipping features to an administrative polygon, or nearest-feature search. The tree stores minimum bounding rectangles and adapts its node shapes to the data, so it handles the wild density gradient between a dense city core and empty countryside without the wasted cells a fixed grid would spend on the empty half. The catch is that it returns candidates, never final answers, and it offers no free notion of a neighbour cell to aggregate into — so it is a poor fit for “count features per area” questions.

import rtree
from shapely.geometry import shape

idx = rtree.index.Index()
geoms: dict[int, object] = {}
for fid, feature in enumerate(features):
    geom = shape(feature["geometry"])
    geoms[fid] = geom
    idx.insert(fid, geom.bounds)  # bounds == (minx, miny, maxx, maxy)

# Coarse candidates by MBR overlap, then refine to an exact answer.
window = (13.30, 52.45, 13.50, 52.55)  # (min_lon, min_lat, max_lon, max_lat)
candidates = [geoms[i] for i in idx.intersection(window)]

The construction details — streaming a PBF into a disk-backed index, the flex_mem versus sparse_file_array location store, and STR-packing for read-heavy trees — belong to the build guide and are not repeated here.

Quadkey and grids: when tiles rule Jump to heading#

Choose a quadkey or an equivalent fixed grid when the output must line up with map tiles, or when you need a dirt-cheap, embarrassingly parallel bucketing key and can tolerate area distortion. A quadkey is a pure function of longitude, latitude, and zoom — no tree to build, no state to persist, and the key is a short string you can group by in any database. It is the right prefilter in front of an exact test: bucket features by tile, and a spatial query only has to open the handful of tiles its window touches. The weaknesses are the $\cos(\varphi)$ area skew and boundary fragmentation — a feature straddling a tile edge belongs to every tile it touches and must be de-duplicated on read.

import mercantile

def quadkey_at(lon: float, lat: float, zoom: int = 14) -> str:
    """Map a coordinate to its slippy-map tile and quadkey string."""
    tile = mercantile.tile(lng=lon, lat=lat, zoom=zoom)
    return mercantile.quadkey(tile)  # e.g. '12020210233' — z digits long

# Bucket POIs into tiles for a coarse, tile-aligned prefilter.
from collections import defaultdict
buckets: dict[str, list[int]] = defaultdict(list)
for fid, (lon, lat) in enumerate(points):
    buckets[quadkey_at(lon, lat)].append(fid)

Because the quadkey shares the exact tiling scheme of the slippy-map pyramid, this is also the family that hands off cleanly to a tile renderer without a reprojection step.

H3 and S2: when uniform area matters Jump to heading#

Choose H3 when the question is aggregation — counts, densities, rates, or gridded ML features — and every cell should cover comparable ground so the numbers are comparable across latitudes. H3’s hexagons have a single neighbour distance in all six directions (no diagonal-versus-orthogonal ambiguity), so ring traversal with grid_disk is clean, and the fixed set of resolutions gives you a hierarchy to roll up or drill down. Choose S2 instead when the workload is planet-scale and you want cells whose 64-bit ids fall on a space-filling curve, so a geographic region maps to a small set of contiguous integer ranges — ideal for range scans in an ordered key-value store. Both approximate polygon boundaries, so neither is a substitute for an R-tree when membership must be exact.

import h3
from collections import Counter

def aggregate_h3(points: list[tuple[float, float]], res: int = 8) -> Counter:
    """Count POIs per H3 cell at a chosen resolution."""
    counts: Counter = Counter()
    for lon, lat in points:
        cell = h3.latlng_to_cell(lat, lon, res)  # note: (lat, lon) order
        counts[cell] += 1
    return counts

cell = h3.latlng_to_cell(52.520, 13.404, 8)
ring = h3.grid_disk(cell, 1)  # the cell plus its six neighbours

Picking the resolution number is its own decision — too coarse blurs the signal, too fine leaves most cells empty — and that trade-off is worked through in Choosing an H3 Resolution for OSM Point Aggregation.

Combine them: an R-tree join plus an H3 column Jump to heading#

The most common production answer is not one index but two, because exact joins and coarse rollups are genuinely different questions. Carry an R-tree for geometry-correct spatial joins and, on the very same feature rows, materialize an H3 cell id as an ordinary column so aggregation is a plain GROUP BY. The R-tree lives in memory or on disk as a tree; the H3 id is just a 64-bit integer that travels with the row into Parquet, PostGIS, or a warehouse, where it costs nothing to filter or roll up.

import h3
import rtree

idx = rtree.index.Index()
rows: list[dict] = []
for fid, feat in enumerate(features):
    lon, lat = feat["lon"], feat["lat"]
    idx.insert(fid, (lon, lat, lon, lat))          # exact spatial join path
    rows.append({
        "fid": fid,
        "h3_r8": h3.latlng_to_cell(lat, lon, 8),    # GROUP BY aggregation path
        "tags": feat["tags"],
    })
# idx answers "which features intersect this polygon?"
# rows answer "how many features per hexagon?" — no re-derivation, one pass.

The single-pass build keeps the two indexes consistent by construction: every feature that enters the tree also gets its hex id, so the join layer and the analytics layer can never drift out of agreement about which features exist.

Failure modes and gotchas Jump to heading#

• Approximate membership sold as exact. Snapping a boundary polygon to H3 or S2 cells and treating cell membership as polygon membership double-counts and mis-counts along every edge. Use a discrete grid to bucket, then refine with an exact predicate when correctness at the boundary matters.
• Latitude area skew in grids. A fixed-zoom quadkey grid gives high-latitude cells a small fraction of the equatorial ground area, biasing any per-cell density. Reach for H3 the moment you compute a statistic per cell across a wide latitude span.
• Coordinate order flips. R-tree bounds are (minx, miny, maxx, maxy) — that is (lon, lat) — but h3.latlng_to_cell takes (lat, lon). Mixing the conventions is the classic silent relocation bug; pin the order per call site.
• Resolution or zoom chosen by taste. Both H3 resolution and quadkey zoom trade cell count against statistical stability; choose them from cell area versus expected feature count, not from a round number.
• Pentagon and antimeridian edge cases. H3 has twelve pentagon cells and S2 wraps the antimeridian; distance and neighbour logic must tolerate both rather than assuming a perfect hex or a flat plane.
• One monolithic index for every feature class. Indexing highway=* and amenity=* into the same structure forces every query to wade through irrelevant features. Filter by feature class at insertion and index each class separately.

Choosing a resolution, in depth Jump to heading#

This section has one child guide, which turns the “which H3 resolution?” question into a runnable, reproducible procedure rather than a guess:

• Choosing an H3 Resolution for OSM Point Aggregation — how to pick a resolution 0–15 for aggregating OSM points, trading cell area against cell count and statistical stability, with the H3 area formula and a per-cell counting script.

Frequently Asked Questions Jump to heading#

Related Jump to heading#

• Spatial Indexing for OSM Extracts — the how-to-build companion that constructs the disk-backed R-tree this page tells you when to reach for.
• Choosing an H3 Resolution for OSM Point Aggregation — the resolution decision once you have chosen H3.
• Coordinate Reference Systems in OSM — why WGS 84, Web Mercator, and the sphere make each family answer latitude distortion differently.
• Node-Way-Relation Data Model — the primitives that decide whether a feature is point-, line-, or polygon-indexable.
• Tag Taxonomy & Key-Value Standards — filtering feature classes so each class gets its own index rather than one monolith.

This guide is part of the OSM Data Fundamentals & Architecture section — return there for the full map of the OSM data model, serialization formats, and indexing foundations.