🧭 Concept Performance & scale Joins & relationships Batch Processing · step 5 of 26 Updated

Spatial Indexes Explained: R-trees and Why Spatial Joins Are Fast

Problem statement

Two ways to find which of 500,000 buildings fall inside 400 wards. One takes a week, the other takes a second.

# the obvious way: 200 million comparisons
matches = []
for i, building in buildings.iterrows():
    for j, ward in wards.iterrows():
        if building.geometry.within(ward.geometry):
            matches.append((i, j))
# the same answer, using the index
matches = gpd.sjoin(buildings, wards, predicate="within")

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. 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.

Understanding it explains when spatial operations are fast, when they are not, and what you can do about it.

Quick answer

A spatial index is a tree of bounding boxes that lets you skip almost every candidate:

  1. Every geometry has a bounding box — a cheap, four-number approximation
  2. Boxes are grouped into a tree (R-tree), so whole branches can be ruled out at once
  3. Query the index to get candidates, then run the exact predicate on those few
  4. GeoPandas exposes it as gdf.sindex, built lazily and cached
  5. sjoin, clip, overlay and sjoin_nearest all use it automatically
import geopandas as gpd

buildings = gpd.read_file("data/raw/buildings.gpkg")
wards = gpd.read_file("data/ref/wards.gpkg").to_crs(buildings.crs)

# the index is built on first access and cached on the GeoDataFrame
print(type(buildings.sindex))

ward = wards.geometry.iloc[0]
candidates = buildings.sindex.query(ward, predicate="intersects")
print(f"{len(candidates)} candidates out of {len(buildings):,} features")

hits = buildings.iloc[candidates]
exact = hits[hits.geometry.within(ward)]
print(f"{len(exact)} exact matches")

That two-stage pattern — cheap filter, then exact test — is the whole idea, and it is what every spatial database does too.

Filter, then refine

A query rectangle against bounding boxes, showing candidates from the index and exact matches after refinement.
The index returns boxes that might match; the exact predicate decides which actually do.

Step-by-step solution

Layered R-tree structure from root node through internal nodes to leaf boxes and geometries.
An R-tree groups boxes into boxes — so one comparison can discard thousands of features.

Bounding boxes: the cheap approximation

Testing whether two complex polygons intersect is expensive: GEOS has to node their edges and compare segments. Testing whether two rectangles overlap is four comparisons.

from shapely.geometry import Polygon
import time

a = Polygon([(0, 0), (100, 0), (100, 100), (0, 100)]).buffer(0)
complex_a = a.buffer(10, quad_segs=64)         # many vertices
complex_b = complex_a.buffer(5)

t0 = time.perf_counter()
for _ in range(10_000):
    complex_a.intersects(complex_b)
exact_s = time.perf_counter() - t0

box_a, box_b = complex_a.bounds, complex_b.bounds
t0 = time.perf_counter()
for _ in range(10_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:.1f} ms   ratio {exact_s/box_s:.0f}×")

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.

import geopandas as gpd

gdf = gpd.read_file("data/raw/buildings.gpkg")
index = gdf.sindex

print("index type   :", type(index).__name__)
print("features     :", len(gdf))
print("index bounds :", [round(v, 1) for v in gdf.total_bounds])

GeoPandas 0.13+ uses Shapely 2's STRtree, a packed Hilbert R-tree: it sorts geometries along a space-filling curve so that neighbours in the tree are neighbours in space, then packs them into full nodes bottom-up. That makes construction fast and queries efficient, at the cost of being read-only — which suits GeoPandas, where a frame's geometry does not change under the index.

Querying the index directly 

import geopandas as gpd
from shapely.geometry import box

buildings = gpd.read_file("data/raw/buildings.gpkg")
area = box(325000, 673000, 326000, 674000)

# bounding-box candidates only
candidates = buildings.sindex.query(area)
print(f"{len(candidates)} candidates")

# ask the index to apply the exact predicate too (Shapely 2)
exact = buildings.sindex.query(area, predicate="intersects")
print(f"{len(exact)} exact intersects")

# many query geometries at once — returns (input index, tree index) pairs
pairs = buildings.sindex.query(wards.geometry, predicate="intersects")
print(pairs.shape)      # (2, n_pairs)

# nearest neighbours
nearest = buildings.sindex.nearest(wards.geometry, return_distance=True)

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 

import geopandas as gpd

# all of these build and use an index internally
joined = gpd.sjoin(buildings, wards, predicate="within")
clipped = gpd.clip(buildings, wards.geometry.iloc[0])
overlaid = gpd.overlay(parcels, zones, how="intersection")
near = gpd.sjoin_nearest(schools, roads, distance_col="dist_m", max_distance=500)

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, a predicate the index does not support — falls back to the quadratic behaviour.

Measure the difference 

import time
import geopandas as gpd

buildings = gpd.read_file("data/raw/buildings.gpkg")
wards = gpd.read_file("data/ref/wards.gpkg").to_crs(buildings.crs)
sample = buildings.sample(2_000, random_state=0)

t0 = time.perf_counter()
brute = [
    (i, j)
    for i, b in zip(sample.index, sample.geometry)
    for j, w in zip(wards.index, wards.geometry)
    if b.within(w)
]
brute_s = time.perf_counter() - t0

t0 = time.perf_counter()
indexed = gpd.sjoin(sample, wards, predicate="within")
index_s = time.perf_counter() - t0

print(f"brute force : {brute_s:6.2f} s  ({len(sample) * len(wards):,} comparisons)")
print(f"sjoin       : {index_s:6.2f} s  ({len(indexed)} matches)")
print(f"speed-up    : {brute_s / index_s:.0f}×")

On 2,000 buildings and 400 wards the ratio is already large; at 500,000 buildings the brute-force version is not something you would wait for.

When the index does not help 

import geopandas as gpd
from shapely.geometry import box

# 1. everything overlaps the query — the filter discards nothing
huge = box(*buildings.total_bounds)
print(len(buildings.sindex.query(huge)), "candidates of", len(buildings))

# 2. one enormous geometry whose bbox covers everything
national_boundary = gpd.read_file("data/ref/country.gpkg").geometry.iloc[0]
print(len(buildings.sindex.query(national_boundary)), "candidates — the bbox is the country")

# 3. building the index costs more than the query saves
tiny = buildings.head(50)

The pathological case is a geometry whose bounding box is much larger than the geometry itself — a long diagonal river, a coastline, a multipart national boundary. Its box overlaps almost everything, so the filter stage passes almost everything through to the expensive exact test.

The fix is to break such geometries up:

import geopandas as gpd

# split a sprawling multipart boundary into its parts, so each has a tight bbox
parts = gpd.GeoDataFrame(geometry=[national_boundary], crs=buildings.crs)
parts = parts.explode(index_parts=False, ignore_index=True)
print(f"1 feature → {len(parts)} parts with much tighter boxes")

hits = gpd.sjoin(buildings, parts, predicate="within")

Indexes in files and databases 

import geopandas as gpd

# GeoPackage carries an R-tree in the file, so bbox reads are cheap
subset = gpd.read_file("data/raw/buildings.gpkg", bbox=(325000, 673000, 326000, 674000))

# FlatGeobuf puts a packed Hilbert R-tree at the front of the file — works over HTTP
remote = gpd.read_file("/vsicurl/https://example.org/buildings.fgb",
                       bbox=(325000, 673000, 326000, 674000))
-- PostGIS: a GIST index is a spatial index on the geometry column
CREATE INDEX buildings_geom_idx ON buildings USING GIST (geom);
ANALYZE buildings;

-- the planner then uses it for && (bbox overlap) before the exact predicate
EXPLAIN ANALYZE
SELECT b.id FROM buildings b JOIN wards w ON ST_Within(b.geom, w.geom);

Note that ST_Within in PostGIS is internally "bbox overlap AND exact within" — the same two-stage pattern, with the index doing the first half. A spatial query that is slow in PostGIS is very often a missing GIST index.

Code examples

Example 1: a manual two-stage join, to see the mechanism 

import geopandas as gpd
import numpy as np

def manual_sjoin(left: gpd.GeoDataFrame, right: gpd.GeoDataFrame, predicate="within"):
    """What sjoin does, written out."""
    tree_index = right.sindex

    # stage 1 — bounding-box filter, vectorised over all left geometries
    left_idx, right_idx = tree_index.query(left.geometry, predicate=None)
    print(f"stage 1: {len(left_idx):,} candidate pairs "
          f"(brute force would be {len(left) * len(right):,})")

    # stage 2 — exact predicate on the survivors only
    keep = np.array([
        getattr(left.geometry.iloc[i], predicate)(right.geometry.iloc[j])
        for i, j in zip(left_idx, right_idx)
    ])
    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)

Example 2: the cost of building the index 

import time
import geopandas as gpd

gdf = gpd.read_file("data/raw/buildings.gpkg")

t0 = time.perf_counter()
_ = gdf.sindex                       # first access builds it
build_s = time.perf_counter() - t0

t0 = time.perf_counter()
_ = gdf.sindex                       # cached
cached_s = time.perf_counter() - t0

print(f"build : {build_s:.3f} s for {len(gdf):,} features")
print(f"cached: {cached_s*1e6:.0f} µs")

Because the index is cached on the frame, reusing the same GeoDataFrame across many queries is far cheaper than re-reading it — and copying or filtering the frame discards the cache.

Example 3: keep the index by filtering carefully 

import geopandas as gpd

buildings = gpd.read_file("data/raw/buildings.gpkg")
_ = buildings.sindex                                  # built once

# a slice creates a new frame, which needs its own index
residential = buildings[buildings["class"] == "residential"]

# for many queries against the same subset, build once and reuse
residential = residential.copy()
_ = residential.sindex

for ward in wards.geometry:
    hits = residential.sindex.query(ward, predicate="within")

The mistake to avoid is filtering inside a loop: each new frame rebuilds the tree, which can dominate the runtime.

Example 4: nearest-neighbour queries 

import geopandas as gpd

schools = gpd.read_file("data/raw/schools.gpkg").to_crs(27700)
roads = gpd.read_file("data/raw/roads.gpkg").to_crs(27700)

# high level
joined = gpd.sjoin_nearest(schools, roads[["road_id", "geometry"]],
                           how="left", distance_col="dist_m", max_distance=1000)
print(joined[["name", "road_id", "dist_m"]].head())

# low level, when you want the distances only
tree_idx, distances = roads.sindex.nearest(
    schools.geometry, return_distance=True, max_distance=1000)
print(distances[:5])

Nearest queries use the same tree, walking it in order of box distance so that the true nearest is found without testing every feature.

Explanation

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.

Bar chart of query time with and without an index at several dataset sizes.
Without an index the work grows with the product of the two layers; with one it barely moves.

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 noding hundreds of edges. If the filter stage discards 99.9% of candidates, the expensive test runs on a thousandth 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 diagonal river or a multipart archipelago has a box many times larger than the geometry. For those, the filter passes far too many candidates through, and performance approaches brute force. Exploding multipart geometries, or splitting long linear features into segments, 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 && (bbox overlap) before any exact predicate. Learning the pattern once explains performance in all four.

Edge cases or notes

  • The index is on bounding boxes, not geometry: It can never be more selective than the boxes allow.
  • sindex is 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.
  • query returns 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.
  • PostGIS needs ANALYZE: An index the planner does not know about will not be used.

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.