Spatial Indexes Explained: R-trees and Why Spatial Joins Are Fast
On this page
Problem statement
Two ways to find which of Leicestershire's 114,187 mapped buildings fall inside each of its wards. Measured on 2,000 of them:
sample = buildings.sample(2_000, random_state=0)
t0 = time.perf_counter()
matches = [(i, j) for i, b in zip(sample.index, sample.geometry) # the obvious way
for j, w in zip(wards.index, wards.geometry) if b.within(w)]
loop_s = time.perf_counter() - t0
t0 = time.perf_counter()
joined = gpd.sjoin(sample, wards, predicate="within") # the same answer, via the index
join_s = time.perf_counter() - t0
print(f"loop : {loop_s:6.2f} s, {len(sample) * len(wards):,} comparisons, {len(matches):,} matches")
print(f"sjoin: {join_s:6.3f} s, {len(joined):,} matches — {loop_s / join_s:,.0f}× faster")loop : 0.73 s, 322,000 comparisons, 1,990 matches
sjoin: 0.017 s, 1,990 matches — 44× faster
The difference is not that GeoPandas is written in C — it is that sjoin builds a spatial index first, and the index turns a quadratic problem into a nearly linear one. At all 114,187 buildings the loop would take about fifty-seven times as long again; the join barely notices. The same structure is behind sindex.query, PostGIS's GIST indexes, GeoPackage's R-tree tables, and every fast clip, overlay and nearest you have ever run.
Quick answer
A spatial index is a tree of bounding boxes that lets you skip almost every candidate:
- Every geometry has a bounding box — a cheap, four-number approximation
- Boxes are grouped into a tree (R-tree), so whole branches can be ruled out at once
- Query the index to get candidates, then run the exact predicate on those few
- GeoPandas exposes it as
gdf.sindex, built lazily and cached sjoin,clip,overlayandsjoin_nearestall use it automatically
print(type(buildings.sindex).__module__, type(buildings.sindex).__name__)
ward = wards.geometry.iloc[0]
candidates = buildings.sindex.query(ward) # bounding boxes only
wrong = buildings.sindex.query(ward, predicate="within") # ward within a building?
exact = buildings.sindex.query(ward, predicate="contains") # ward contains the building
print(f"{wards['WD24NM'].iloc[0]}: {len(candidates):,} bounding-box candidates of {len(buildings):,}")
print(f"predicate='within': {len(wrong):,} | predicate='contains': {len(exact):,} buildings inside the ward")geopandas.sindex SpatialIndex
Ambien: 8,474 bounding-box candidates of 114,187
predicate='within': 0 | predicate='contains': 1,379 buildings inside the ward
That two-stage pattern — cheap filter, then exact test — is the whole idea, and it is what every spatial database does too. Note the direction of the predicate: sindex.query(geometry, predicate=p) tests p(geometry, tree_geometry) — the ward within a building, which is never true. To find the buildings within the ward, ask whether the ward contains them.
Filter, then refine
From bounding boxes to R-trees
Bounding boxes: the cheap approximation
Testing whether two complex polygons intersect is expensive: GEOS has to compare their edges. Testing whether two rectangles overlap is four comparisons. Two neighbouring real wards:
a, b = wards.geometry.iloc[0], wards.geometry.iloc[1]
print("vertices:", shapely.get_num_coordinates(a), "and", shapely.get_num_coordinates(b))
t0 = time.perf_counter()
for _ in range(2_000):
a.intersects(b)
exact_s = time.perf_counter() - t0
box_a, box_b = a.bounds, b.bounds
t0 = time.perf_counter()
for _ in range(2_000):
box_a[0] <= box_b[2] and box_a[2] >= box_b[0] and box_a[1] <= box_b[3] and box_a[3] >= box_b[1]
box_s = time.perf_counter() - t0
print(f"exact: {exact_s * 1000:.1f} ms bbox: {box_s * 1000:.2f} ms ratio {exact_s / box_s:,.0f}×")vertices: 3411 and 1386
exact: 38.4 ms bbox: 0.20 ms ratio 190×
A bounding-box test never gives a false negative — if the boxes do not overlap, the geometries cannot — but it does give false positives, which is why the exact test still has to run on the survivors.
The R-tree: boxes of boxes
An R-tree groups nearby bounding boxes into a parent box, groups those parents into grandparents, and so on up to a single root. A query walks down from the root, and any node whose box does not intersect the query is discarded together with everything beneath it.
GeoPandas uses Shapely 2's STRtree: an R-tree bulk-loaded with the Sort-Tile-Recursive algorithm, which sorts the boxes into vertical slices and then into tiles so that each full node holds boxes that are near each other, built bottom-up. Construction is fast and queries are efficient, at the cost of being read-only — which suits GeoPandas, where a frame's geometry does not change under the index. (FlatGeobuf uses a close relative, a packed Hilbert R-tree, which orders the boxes along a space-filling curve instead.)
Querying the index directly
area = box(456000, 303000, 460000, 307000) # a 4 km square of Leicester
print(len(buildings.sindex.query(area)), "bounding-box candidates")
print(len(buildings.sindex.query(area, predicate="intersects")), "exact intersects")
pairs = buildings.sindex.query(wards.geometry, predicate="intersects") # many query geometries at once
print("pairs array:", pairs.shape, "(ward index, building index)")
nearest, distance = wards.sindex.nearest(buildings.geometry.iloc[:5], return_distance=True)
print("nearest ward to each of five buildings:", nearest[1].tolist(), distance.round(1).tolist())5670 bounding-box candidates
5669 exact intersects
pairs array: (2, 114091) (ward index, building index)
nearest ward to each of five buildings: [122, 126, 130, 130, 85] [0.0, 0.0, 0.0, 0.0, 0.0]
query(..., predicate=...) is the one to reach for: it does the filter and the refine inside Shapely, without a Python loop over candidates.
Where the index is used for you
sjoin, clip, overlay and sjoin_nearest all build and use an index internally. You rarely need to touch sindex yourself; knowing it is there matters because it explains the performance cliff. Any operation that cannot use it — a Python loop, an apply over rows — falls back to the quadratic behaviour of the loop at the top of the page.
When the index does not help
The filter's usefulness depends on boxes being tight. A long, diagonal geometry has a box far larger than itself — here, a 100 m corridor along the River Soar:
soar = waterways[waterways["name"] == "River Soar"].union_all().buffer(100)
print(f"corridor: {soar.area / 1e6:.1f} km², its bounding box: {shapely.box(*soar.bounds).area / 1e6:.0f} km²")
t0 = time.perf_counter()
candidates = buildings.sindex.query(soar)
hits = buildings.sindex.query(soar, predicate="intersects")
whole_s = time.perf_counter() - t0
print(f"one geometry : {len(candidates):,} candidates for {len(hits):,} hits ({whole_s:.2f} s)")
pieces = gpd.GeoSeries(shapely.intersection(grid, soar), crs=27700) # cut on a 2 km grid
pieces = pieces[~pieces.is_empty]
t0 = time.perf_counter()
left, right = buildings.sindex.query(pieces, predicate=None)
hits_split = np.unique(buildings.sindex.query(pieces, predicate="intersects")[1])
split_s = time.perf_counter() - t0
print(f"{len(pieces)} pieces : {len(np.unique(right)):,} candidates for {len(hits_split):,} hits ({split_s:.2f} s)")corridor: 12.8 km², its bounding box: 549 km²
one geometry : 50,047 candidates for 860 hits (0.03 s)
39 pieces : 4,607 candidates for 860 hits (0.00 s)
The whole corridor's box passed tens of thousands of buildings to the exact test; cut into 2 km pieces, each with a tight box, the same hits came from a fraction of the candidates. The same applies to a multipart boundary — explode() it first — and to any long linear feature.
Indexes in files and databases
A GeoPackage carries an R-tree in the file, so a bbox read touches only the rows it needs; FlatGeobuf puts its index at the front of the file, which also works over HTTP through GDAL's /vsicurl/:
for path in ("data/buildings.gpkg", "data/buildings.fgb"):
t0 = time.perf_counter()
everything = gpd.read_file(path)
full_s = time.perf_counter() - t0
t0 = time.perf_counter()
subset = gpd.read_file(path, bbox=tuple(area.bounds))
bbox_s = time.perf_counter() - t0
print(f"{Path(path).suffix:5} full read {len(everything):,} rows in {full_s:.2f} s; bbox read {len(subset):,} rows in {bbox_s:.3f} s").gpkg full read 114,187 rows in 0.16 s; bbox read 5,669 rows in 0.012 s
.fgb full read 114,187 rows in 0.16 s; bbox read 5,669 rows in 0.010 s
In PostGIS the index is a GIST index on the geometry column, and the planner uses it for the && bounding-box test that ST_Within performs before its exact test. With the buildings and wards loaded into a local PostGIS:
query = "SELECT count(*) FROM buildings b JOIN wards w ON ST_Within(b.geom, w.geom)"
with engine.begin() as conn:
print("indexes after to_postgis:", [r[0] for r in conn.exec_driver_sql("SELECT indexname FROM pg_indexes WHERE tablename = 'buildings'")])
for label, statement in (("no index", "DROP INDEX IF EXISTS idx_buildings_geom; ANALYZE buildings"),
("GIST index", "CREATE INDEX buildings_geom_idx ON buildings USING GIST (geom); ANALYZE buildings")):
if statement:
conn.exec_driver_sql(statement)
plan = [row[0] for row in conn.exec_driver_sql("EXPLAIN ANALYZE " + query)]
used = next((line.strip() for line in plan if "Index" in line), "no index scan")
total = next(line for line in plan if line.startswith("Execution Time"))
print(f"{label:10} {total:28} | {used[:70]}")indexes after to_postgis: ['idx_buildings_geom']
no index Execution Time: 10708.467 ms | -> Index Scan using idx_wards_geom on wards w (cost=0.15..12.67 rows
GIST index Execution Time: 371.292 ms | -> Index Scan using buildings_geom_idx on buildings b (cost=0.28..14
to_postgis had already created a GIST index — GeoAlchemy2 adds one by default — so the first run drops it to measure the difference. Even with the wards still indexed, the join took about ten seconds without the buildings' index and well under half a second with it. A spatial query that is slow in PostGIS is very often a missing GIST index, which is worth checking even when the tool that loaded the table might have made one.
Two-stage joins on Leicestershire's buildings
The data used on this page
Leicestershire's buildings and waterways from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL), and the wards of Leicester and the county's districts from the ONS December 2024 boundaries (OGL), all in British National Grid. The PostGIS examples use a throwaway container, as other guides on this site do.
import os
import time
from pathlib import Path
import geopandas as gpd
import numpy as np
import requests
import shapely
from shapely.geometry import box
from sqlalchemy import create_engine
Path("data").mkdir(exist_ok=True)
if not Path("data/leicestershire.shp.zip").exists():
Path("data/leicestershire.shp.zip").write_bytes(requests.get(
"https://download.geofabrik.de/europe/united-kingdom/england/leicestershire-260101-free.shp.zip", timeout=600).content)
LADS = ["E06000016", "E07000129", "E07000130", "E07000131", "E07000132", "E07000133", "E07000134", "E07000135"]
if not Path("data/wards.geojson").exists():
Path("data/wards.geojson").write_bytes(requests.get(
"https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
"Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
params={"where": "LAD24CD IN ('" + "','".join(LADS) + "')", "outFields": "WD24NM,LAD24NM",
"outSR": 27700, "f": "geojson"}, timeout=180).content)
OSM = "/vsizip/data/leicestershire.shp.zip/"
buildings = gpd.read_file(OSM + "gis_osm_buildings_a_free_1.shp", columns=["osm_id"]).to_crs(27700)
waterways = gpd.read_file(OSM + "gis_osm_waterways_free_1.shp", columns=["name"]).to_crs(27700)
wards = gpd.read_file("data/wards.geojson")
for suffix in ("gpkg", "fgb"):
if not Path(f"data/buildings.{suffix}").exists():
buildings.to_file(f"data/buildings.{suffix}")
minx, miny, maxx, maxy = buildings.total_bounds
xs, ys = np.meshgrid(np.arange(minx, maxx, 2000), np.arange(miny, maxy, 2000))
grid = shapely.box(xs.ravel(), ys.ravel(), xs.ravel() + 2000, ys.ravel() + 2000) # 2 km squares
os.environ.update(PGPORT="55433", PGPASSWORD="guide")
if os.system("docker inspect guide-postgis >/dev/null 2>&1") != 0:
os.system("docker run -d --name guide-postgis -e POSTGRES_USER=user -e POSTGRES_PASSWORD=guide "
"-e POSTGRES_DB=gis -p 127.0.0.1:55433:5432 postgis/postgis:16-3.4 >/dev/null")
time.sleep(10)
engine = create_engine("postgresql+psycopg://user:guide@localhost:55433/gis")
buildings.rename_geometry("geom").to_postgis("buildings", engine, if_exists="replace", index=False)
wards.rename_geometry("geom").to_postgis("wards", engine, if_exists="replace", index=False)
print(f"{len(buildings):,} buildings, {len(wards)} wards in {wards['LAD24NM'].nunique()} districts, "
f"{len(waterways):,} waterways")114,187 buildings, 161 wards in 8 districts, 4,701 waterways
Example 1: a manual two-stage join, to see the mechanism
def manual_sjoin(left, right, predicate="within"):
"""What sjoin does, written out."""
left_idx, right_idx = right.sindex.query(left.geometry, predicate=None) # stage 1: boxes
print(f"stage 1: {len(left_idx):,} candidate pairs (brute force would be {len(left) * len(right):,})")
keep = getattr(shapely, predicate)(left.geometry.values[left_idx], right.geometry.values[right_idx]) # stage 2
print(f"stage 2: {keep.sum():,} exact matches")
return left.iloc[left_idx[keep]].assign(index_right=right.index[right_idx[keep]])
result = manual_sjoin(buildings.sample(5_000, random_state=0), wards)
print("same as sjoin:", len(result) == len(gpd.sjoin(buildings.sample(5_000, random_state=0), wards, predicate="within")))stage 1: 10,562 candidate pairs (brute force would be 805,000)
stage 2: 4,975 exact matches
same as sjoin: True
Example 2: the cost of building the index, and when it is kept
fresh = buildings.copy()
t0 = time.perf_counter()
_ = fresh.sindex
build_s = time.perf_counter() - t0
t0 = time.perf_counter()
_ = fresh.sindex
cached_s = time.perf_counter() - t0
print(f"build: {build_s:.3f} s for {len(fresh):,} features | cached: {cached_s * 1e6:.0f} µs")
subset = fresh[fresh.index % 2 == 0]
print("index kept on a filtered frame:", subset.has_sindex, "| on a .copy():", fresh.copy().has_sindex)build: 0.023 s for 114,187 features | cached: 160 µs
index kept on a filtered frame: False | on a .copy(): False
Because the index is cached on the frame, reusing the same GeoDataFrame across many queries is far cheaper than re-reading it — and filtering or copying the frame starts without one, so build it once on the frame you will query repeatedly rather than filtering inside a loop.
Example 3: nearest-neighbour queries
sample = buildings.sample(1_000, random_state=1)
rivers = waterways[waterways["name"].notna()]
joined = gpd.sjoin_nearest(sample, rivers, how="left", distance_col="dist_m", max_distance=1_000)
print(f"{joined['dist_m'].notna().sum()} of {len(sample)} buildings have a named waterway within 1 km; "
f"median {joined['dist_m'].median():.0f} m")699 of 1000 buildings have a named waterway within 1 km; median 415 m
Nearest queries use the same tree, walking it in order of box distance so that the true nearest is found without testing every feature, and max_distance lets them stop early.
Why grouping boxes makes queries fast
The R-tree was described by Antonin Guttman in 1984, and the idea has not needed changing since: group objects by their bounding rectangles into a balanced tree, so that a query can discard an entire subtree with one rectangle comparison.
The reason it works so well is the asymmetry between the two tests. A bounding-box overlap is four floating-point comparisons. An exact polygon intersection may require comparing thousands of edges. If the filter stage discards almost every candidate, the expensive test runs on a tiny share of the pairs, and the total cost is dominated by the cheap operation.
Crucially, the filter is conservative: a bounding box always contains its geometry, so a box that does not overlap the query guarantees the geometry does not either. There are no false negatives, which is what makes the two-stage approach exact rather than approximate. The false positives — boxes that overlap where the geometries do not — are what the refinement stage exists to remove.
This also explains the failure mode. The filter's usefulness depends on boxes being tight: a compact building has a box close to its own area, while a river corridor has a box many times larger than itself. For those, the filter passes far too many candidates through, and performance approaches brute force. Exploding multipart geometries, or splitting long features into pieces, restores tight boxes and the speed that goes with them.
The same structure appears everywhere in the stack because the same reasoning applies at every level. GeoPandas builds an STRtree in memory. GeoPackage stores an R-tree in an SQLite virtual table. FlatGeobuf packs one at the front of the file so a client can range-request the right bytes. PostGIS builds a GIST index on the geometry column and its planner uses && before any exact predicate. Learning the pattern once explains performance in all four.
What the index cannot do
- The index is on bounding boxes, not geometry: It can never be more selective than the boxes allow.
sindexis built lazily and cached: The first spatial operation on a frame pays for construction; subsequent ones do not.- Slicing a frame drops the cached index: Build the index on the frame you will query repeatedly, not inside a loop.
- Shapely 2's STRtree is immutable: Adding geometries means rebuilding. That is fine for analysis, wrong for a live editing session.
queryreturns positional indices: Use.iloc, not.loc, unless you map them back to labels.- Long diagonal geometries defeat it: Explode multiparts and consider segmentising long lines.
- Run
ANALYZEafter loading PostGIS data: the planner chooses between the index and a scan from table statistics, and fresh tables have none.
Internal links
- How to Perform a Spatial Join in Python (GeoPandas)
- How to Speed Up GeoPandas: Tips for Large Datasets
- How to Find the Nearest Point in GeoPandas
- How to Select Features by Location in GeoPandas
- Spatial Join Returns Empty Results in GeoPandas: How to Fix It
- GIS Vector File Formats Compared: Shapefile, GeoPackage, GeoJSON, Parquet
- Discrete global grids explained — a cell id as an index you can store in a column
FAQ
What is a spatial index?
A tree of bounding boxes that lets a query discard most features without testing them. GeoPandas exposes it as gdf.sindex; it is the reason sjoin is fast.
Do I need to build the index myself?
No. It is created on first use and cached on the GeoDataFrame, and sjoin, clip, overlay and sjoin_nearest all use it automatically.
Why is my spatial join still slow?
Usually one of three things: a geometry whose bounding box covers most of the layer, a new frame in each loop iteration forcing a rebuild, or a Python loop that never reaches the index at all.
What is an R-tree exactly?
A balanced tree whose nodes are bounding rectangles containing their children. Shapely 2 uses a packed Hilbert variant, which sorts geometries along a space-filling curve before packing them.
Does the index give exact results?
The index stage is approximate — it returns candidates — but because a bounding box always contains its geometry there are no false negatives, and the exact predicate on the candidates makes the overall result exact.
Do file formats have indexes too?
GeoPackage stores an R-tree in the file, and FlatGeobuf puts one at the front so bbox reads work over HTTP. Shapefile has an optional .qix; GeoJSON has none, so bbox reads still parse the whole file.
How do I speed up queries against one huge geometry?
Explode it into parts so each has a tight bounding box, or clip it into tiles. A single geometry whose box spans the dataset makes the filter stage useless.