R-tree vs Grid Index for Point-in-Polygon Queries
For point-in-polygon joins the GeoPandas R-tree (.sindex behind gpd.sjoin) is the right default because it adapts to skewed density and variable polygon sizes; a hand-built uniform grid only wins when points are near-uniform, polygons are similar in size, and you can tune the cell size precisely.
Jump to heading Why this matters
Assigning millions of points to the polygon that contains them — GPS pings to districts, addresses to sales territories, sensor readings to catchment zones — is the workhorse operation of spatial dashboards, and it runs on every filter change. The exact contains test is expensive, so both contenders exist to shrink the candidate set before that test runs. They shrink it differently. An R-tree groups polygon bounding boxes into a balanced tree and descends it per query; a uniform grid slices space into fixed cells and hashes each point straight to its bucket. Which structure is faster is entirely a function of how your data is shaped, and picking wrong can cost an order of magnitude. This comparison sits under the Spatial Index Types Reference, which catalogues both structures in the broader Spatial Data Reference.
Jump to heading Prerequisites
- Python 3.10+,
geopandas>=1.0,shapely>=2.0, andnumpy>=1.24. - Both layers reprojected to a shared projected CRS (metres). Building either index in geographic degrees distorts cell sizes and bounding boxes — see the CRS & Coordinate Systems Reference.
- A points GeoDataFrame and a polygons GeoDataFrame that already share the same CRS.
Jump to heading Step-by-step comparison
Jump to heading Step 1 — Point-in-polygon with the R-tree
The R-tree path is one call. gpd.sjoin consults the .sindex of the right-hand frame, prunes candidate polygons per point, then runs the exact within test on survivors.
import geopandas as gpd
import numpy as np
# Polygons: sales territories in WGS 84, reprojected to UTM 33N metres
territories = gpd.read_file("territories.gpkg").to_crs("EPSG:32633")
# Points: 200k customer locations around Berlin (lon 13.40, lat 52.52)
rng = np.random.default_rng(7)
pts = gpd.GeoDataFrame(
{"cust_id": np.arange(200_000)},
geometry=gpd.points_from_xy(
rng.normal(13.404, 0.09, 200_000),
rng.normal(52.520, 0.06, 200_000),
crs="EPSG:4326",
),
).to_crs("EPSG:32633")
rtree_join = gpd.sjoin(
pts, territories[["terr_id", "geometry"]],
how="left", predicate="within",
)
print(rtree_join["terr_id"].notna().sum(), "points assigned")
Jump to heading Step 2 — Point-in-polygon with a manual uniform grid
A grid index hashes each polygon into every cell its bounding box touches, then hashes each point to its single cell to fetch candidates. Cell size is the one tuning knob that decides everything.
from collections import defaultdict
from shapely.geometry import Point
def build_grid_index(polys: gpd.GeoDataFrame, cell: float) -> dict:
"""Map each grid cell (col, row) -> list of polygon positions."""
grid = defaultdict(list)
for pos, geom in enumerate(polys.geometry.values):
minx, miny, maxx, maxy = geom.bounds
for col in range(int(minx // cell), int(maxx // cell) + 1):
for row in range(int(miny // cell), int(maxy // cell) + 1):
grid[(col, row)].append(pos)
return grid
def grid_point_in_polygon(pts, polys, grid, cell):
poly_geoms = polys.geometry.values
poly_ids = polys["terr_id"].values
out = np.full(len(pts), None, dtype=object)
xs, ys = pts.geometry.x.values, pts.geometry.y.values
for i in range(len(pts)):
key = (int(xs[i] // cell), int(ys[i] // cell))
p = Point(xs[i], ys[i])
for pos in grid.get(key, ()): # candidates in this cell only
if poly_geoms[pos].contains(p): # exact test
out[i] = poly_ids[pos]
break
return out
# Cell sized near the median territory width (~1.5 km in EPSG:32633 metres)
CELL = 1500.0
grid = build_grid_index(territories, CELL)
grid_ids = grid_point_in_polygon(pts, territories, grid, CELL)
The grid replaces tree descent with a dictionary lookup, but it only pays off if each cell holds a handful of polygons. Get the cell size wrong and each lookup returns dozens of candidates, dragging the exact test back toward a linear scan.
Jump to heading Step 3 — Benchmark both on identical data
import time
t0 = time.perf_counter()
gpd.sjoin(pts, territories, how="left", predicate="within")
rtree_ms = (time.perf_counter() - t0) * 1000
t0 = time.perf_counter()
grid = build_grid_index(territories, CELL)
build_ms = (time.perf_counter() - t0) * 1000
t0 = time.perf_counter()
grid_point_in_polygon(pts, territories, grid, CELL)
grid_query_ms = (time.perf_counter() - t0) * 1000
print(f"R-tree total: {rtree_ms:8.0f} ms")
print(f"Grid build: {build_ms:8.0f} ms")
print(f"Grid query (loop): {grid_query_ms:8.0f} ms")
# Vectorised C sjoin typically beats the pure-Python grid loop here;
# a grid written in NumPy/Numba closes the gap on uniform data.
On skewed, real-world data the vectorised C R-tree join usually wins outright, partly because the grid query above is a Python loop. On genuinely uniform points with a well-sized cell, a vectorised grid can edge ahead — the point of the benchmark is to measure your distribution rather than trust a rule of thumb.
Jump to heading Step 4 — Choose by the shape of your data
| Factor | R-tree (.sindex) | Uniform grid |
|---|---|---|
| Point density | Adapts to skew automatically | Degrades toward O(n) in hot cells |
| Polygon size variance | Handles mixed sizes well | Large polygons flood many cells |
| Data mutability | Rebuilt on geometry reassignment | Cheap incremental inserts |
| Build cost | O(n log n), lazy on first use | O(n), but re-tune cell on new data |
| Query cost | O(log n) tree descent | O(1) hash, if cells stay small |
| Tuning required | None | Cell size must match polygon scale |
| Best fit | Default for real-world spatial joins | Uniform density, streaming inserts |
Choose the R-tree unless you can honestly answer “yes” to uniform density, similar polygon sizes, and a stable cell size. Those conditions hold for regular sensor grids and tiled raster footprints; they rarely hold for administrative or commercial polygons.
Jump to heading Verification
Confirm both indexes produce the same assignment before trusting either one’s timing:
rtree_ids = (
gpd.sjoin(pts, territories, how="left", predicate="within")
.sort_index()["terr_id"].values
)
agreement = np.mean(
[a == b for a, b in zip(rtree_ids, grid_ids)]
)
print(f"assignment agreement: {agreement:.4%}")
assert agreement > 0.999, "index disagreement — check cell size and boundary handling"
Agreement below 100% almost always points to boundary points or a grid cell that missed a polygon whose bounding box straddled a cell edge — both covered next.
Jump to heading Edge cases and gotchas
- Points exactly on a boundary: the index only selects candidates; the
withinpredicate decides membership. A point on the shared edge of two territories may match neither underwithinor both underintersects. Snap coordinates to a tolerance (shapely.set_precision(geom, 0.001)) before the join for deterministic results. - Grid cell sizing: a cell far smaller than the median polygon makes large polygons register in hundreds of buckets, inflating build time and memory; a cell far larger collapses candidates back to a linear scan. Target a cell near the median polygon bounding-box width and measure average candidates per cell.
- Memory of dense grids: hashing every polygon into every overlapping cell can explode memory when a few huge polygons cover most of the extent. Cap it by storing only polygon positions (ints), not geometries, in each bucket, and consider a quadtree if polygon sizes span several orders of magnitude.
Jump to heading FAQ
Is a hand-built grid index ever faster than the R-tree in GeoPandas?
Yes, in a narrow case: when points are close to uniformly distributed, polygons are similar in size, and you size the grid cell so each holds only a few polygons. A single hash lookup then beats tree descent. As soon as density is skewed or polygon sizes vary widely, the R-tree wins because it adapts to the data while a fixed grid does not. Vectorising the grid query in NumPy or Numba also matters — the pure-Python loop shown here hands the C-level sjoin an unfair advantage.
How do I choose a grid cell size?
Aim for a cell roughly the size of the median polygon bounding box. Too large and each cell holds many candidates, degrading toward a linear scan; too small and large polygons register in hundreds of cells, inflating memory and build time. Compute territories.geometry.bounds and take the median of (maxx - minx) as a starting cell width, then measure the average candidates returned per lookup and tune toward single digits.
Do points exactly on a polygon boundary behave the same in both indexes?
The index only selects candidates; the exact predicate decides membership, so boundary behaviour depends on the predicate, not the index. With within, a point on a shared edge may match neither polygon; with intersects it may match both. Snap coordinates to a tolerance before the join if you need deterministic assignment, and document which predicate the dashboard uses so downstream counts are reproducible.
Back to Spatial Index Types Reference
Related
- Spatial Index Types Reference — the full catalogue of index structures these two are drawn from
- Speeding Up GeoPandas Spatial Joins with .sindex — tune the R-tree path benchmarked here
- CRS & Coordinate Systems Reference — align both layers before any index lookup
- Query Result Caching — cache the join result so repeated point-in-polygon queries never recompute