Spatial Dashboard Streamlit · Panel · Production
Data Reference

Speeding Up GeoPandas Spatial Joins with the .sindex R-tree

gpd.sjoin already uses the .sindex R-tree for you — the real speed wins come from aligning the CRS first, choosing the right predicate, and shrinking both the polygons and the candidate set before the exact geometric test runs.

Why this matters

A spatial join is the most expensive routine operation in most dashboards: it pairs every point with every polygon it relates to, and it reruns whenever a filter, slider, or map bounds change. GeoPandas is not naive about this — under the hood gpd.sjoin builds an R-tree spatial index and uses it to discard the vast majority of point-polygon pairs before evaluating the exact predicate. The reason joins still feel slow is almost never the index itself; it is the things around it: a CRS mismatch that defeats the bounding-box pre-filter, polygons carrying tens of thousands of vertices, or the large layer placed where the tree gets built over it. This guide fixes each of those, working within the Spatial Data Reference.

How gpd.sjoin uses the .sindex R-treePoints and polygons enter the sjoin pipeline. The right-hand polygon layer's .sindex R-tree is queried with each point's bounding box, producing a small set of candidate pairs. The exact predicate, such as within, then tests only those candidates, yielding the final matched join rows. A note shows that skipping candidate reduction sends every pair to the exact test instead.pointsleft framepolygonsright frame.sindexR-tree bbox prunecandidate pairsfew per pointexactpredicatewithinSkip the R-tree (mismatched CRS, manual loop) → every pair reaches the exact predicate

Prerequisites

  • Python 3.10+, geopandas>=1.0, shapely>=2.0, numpy>=1.24, and the rtree package (or a Shapely 2 build that supplies the tree).
  • A points layer and a polygons layer. This guide reprojects both from EPSG:4326 to EPSG:32633 (UTM 33N) so the index measures metres.
  • Comfort with the difference between the bounding-box pre-filter and the exact predicate, introduced in the Spatial Index Types Reference.

Step-by-step solution

Step 1 — Align the CRS before anything else

The single biggest sjoin performance and correctness bug is a CRS mismatch. If the two layers use different coordinate systems, their bounding boxes live in different number ranges and the R-tree pre-filter either matches everything or nothing.

python
import geopandas as gpd
import numpy as np

# Polygons: catchment zones, WGS 84 on disk → projected metres
zones = gpd.read_file("catchments.gpkg").to_crs("EPSG:32633")

# Points: 500k sensor readings around Berlin (lon 13.40, lat 52.52)
rng = np.random.default_rng(11)
readings = gpd.GeoDataFrame(
    {"sensor_id": np.arange(500_000)},
    geometry=gpd.points_from_xy(
        rng.normal(13.404, 0.10, 500_000),
        rng.normal(52.520, 0.07, 500_000),
        crs="EPSG:4326",
    ),
).to_crs("EPSG:32633")

assert readings.crs == zones.crs, "align CRS or the R-tree prunes nothing"

Step 2 — Let sjoin use the index

gpd.sjoin builds the .sindex of the right-hand frame and consults it automatically. Put the smaller layer on the right so the tree is smaller and cheaper to descend.

python
joined = gpd.sjoin(
    readings,                              # large: left
    zones[["zone_id", "geometry"]],        # small: right → tree built here
    how="left",
    predicate="within",
)
print(joined["zone_id"].notna().sum(), "readings matched to a zone")

To avoid the first user interaction paying the build cost, warm the tree once at startup:

python
_ = zones.sindex          # build the R-tree eagerly, then cache the frame
assert zones.has_sindex

Step 3 — Use query and query_bulk for custom joins

When the join logic is bespoke — nearest zone within a distance, first match only, a custom scoring test — call the index directly. Pass an array of geometries to sindex.query; it returns a two-row array of input positions and tree positions.

python
# Candidate (point_pos, zone_pos) pairs straight from the R-tree
input_pos, tree_pos = zones.sindex.query(
    readings.geometry.values, predicate="intersects"
)

# Run any custom exact test on the candidate pairs only
pt = readings.geometry.values
zn = zones.geometry.values
matches = [
    (i, j) for i, j in zip(input_pos, tree_pos)
    if zn[j].contains(pt[i])
]
print(len(matches), "confirmed containments from",
      len(input_pos), "candidates")

query_bulk was the older name for querying with many geometries and is deprecated in GeoPandas 1.0 — the array form of query above replaces it and returns the same candidate structure.

An index built before a reprojection is not stale — it is wrongTwo orderings of the same three operations are traced. In the correct ordering both frames are reprojected to a common system, the index is built afterwards, and the join's candidate boxes describe exactly the coordinates the exact predicate will then test — so the result is complete. In the incorrect ordering the index is built first and the frame is reprojected afterwards. Reprojection replaces the coordinate arrays but the index's stored bounding boxes are not recomputed, so the candidate search runs against boxes that describe where the geometry used to be. The failure is not an exception: the join returns rows, just not the right ones, and it returns them quickly. What makes this specific bug persistent is that both a row count and a spot check of one feature can look entirely reasonable — the only reliable test is comparing the indexed result against an unindexed run on a small sample.REPROJECT, THEN INDEX, THEN JOINto_crs on both frames.sindex built heresjoin — completeboxes describe thecoordinates being testedINDEX, THEN REPROJECT, THEN JOIN.sindex built hereto_crs replaces thecoordinates, not the boxessjoin — fast and wrongrows are returned,just not the right onesNo exception is raised, and both a row count and a spot check of one feature can look entirely reasonable — which iswhy this survives review so often.The one reliable test is comparing the indexed result against an unindexed run on a small sample. ### Step 4 — Pick the right predicate

The predicate decides both correctness and how many candidates survive the exact test. Match it to the geometric question:

PredicateGeometric meaningTypical join
intersectsGeometries share any point, boundaries includedPoints possibly on zone edges; permissive matching
withinLeft geometry lies fully inside the right, boundary excludedAssign each point to its containing polygon
containsRight geometry lies fully inside the leftFind polygons enclosed by a drawn selection
touchesGeometries share a boundary but no interiorAdjacency between neighbouring polygons

For point-to-polygon assignment, within is the standard choice. Using intersects where within is meant inflates the candidate survivors and can double-count points that land on a shared boundary.

Why the index returns more than the answerA long diagonal river polygon is drawn with its axis-aligned bounding box around it. The box is what the index actually stores, and for a diagonal or elongated shape it covers a great deal of ground the geometry does not: points scattered across the box are all returned as candidates, and only the few that fall inside the polygon itself survive the exact predicate. The ratio between the two counts is the selectivity of the index for this shape, and it is poor precisely for the geometries dashboards care about — rivers, roads, coastlines, administrative strips. This is why the two-pass structure is not an optimisation but the design: the index is a cheap way to eliminate what cannot match, never a way to determine what does, and any code that treats a candidate list as an answer is wrong for every shape that is not roughly square.WHAT THE INDEX STORES, AND WHAT YOU ASKED FORdashed: the bounding box the index stores · solid: the geometry you asked about15 candidates returnedeverything inside the dashed box3 survive the exact testthe other twelve are nowhere near the riverSelectivity is worst for exactly theshapes dashboards care about —rivers, roads, coastlines, strips.The two passes are the design, not an optimisation: the index cheaply eliminates what cannot match and never determines what does. ### Step 5 — Reduce candidates before the exact test

The R-tree gets you to a small candidate set; the exact predicate then walks polygon vertices. Fewer vertices means a faster exact test. Simplify display-only polygons and clip to the working extent before joining.

python
from shapely.geometry import box

# 1. Clip to the current viewport so the tree is smaller
viewport = box(380000, 5810000, 400000, 5830000)
zones_v = zones.clip(viewport)

# 2. Simplify complex boundaries for display-tier joins (metres tolerance)
zones_v = zones_v.copy()
zones_v["geometry"] = zones_v.geometry.simplify(10, preserve_topology=True)

fast_join = gpd.sjoin(readings, zones_v, how="inner", predicate="within")

Only simplify geometry used for display or approximate binning — never for legally precise boundaries such as parcel or tax lines.

Verification

Prove the tuned join is both correct and faster than the un-indexed baseline:

python
import time

# Baseline: force a per-point Python loop (no R-tree)
sample = readings.iloc[:3000]
t0 = time.perf_counter()
brute = [
    next((zid for zid, g in zip(zones.zone_id, zones.geometry)
          if g.contains(p)), None)
    for p in sample.geometry
]
brute_ms = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
idx = gpd.sjoin(sample, zones, how="left", predicate="within")
idx_ms = (time.perf_counter() - t0) * 1000

agree = np.mean(
    [a == b for a, b in zip(brute, idx.sort_index()["zone_id"].values)]
)
print(f"brute={brute_ms:.0f}ms  indexed={idx_ms:.0f}ms  agreement={agree:.2%}")
assert agree > 0.999, "indexed join disagrees with brute force"

The indexed join should be several times faster on this sample and agree with the brute-force result to within boundary-point tolerance. The speed gap grows super-linearly as the layers scale.

Edge cases and gotchas

  • CRS must match before the join: a geographic-vs-projected mismatch does not raise an error — it silently returns empty or all-matching results. Assert left.crs == right.crs and reproject the smaller frame first.
  • Predicate correctness: within excludes boundary points, intersects includes them. Choosing the wrong one changes counts without any error. Decide the convention up front and encode it in a test that snapshots the matched row count.
  • Index invalidation after geometry edits: the cached .sindex describes the geometry as it was when the tree was built. Editing vertices in place leaves the tree stale, so queries return wrong candidates. Reassign the geometry column (gdf = gdf.set_geometry("geometry")) to force a rebuild, then cache the result with Query Result Caching so the rebuilt join is not recomputed on every rerun.

FAQ

Do I need to call .sindex before gpd.sjoin?

No. gpd.sjoin builds and uses the .sindex R-tree of the right-hand GeoDataFrame automatically on first use. You only touch .sindex directly when writing a custom join with sindex.query, or when you want to build the tree ahead of time — for example warming it at dashboard startup with _ = zones.sindex so the first user interaction does not pay the build cost.

Why is my sjoin still slow even though .sindex exists?

The usual causes are a mismatched CRS that makes every bounding box overlap, oversized polygons with tens of thousands of vertices that make the exact predicate expensive, or joining the large frame as the right-hand side so the tree is built over the wrong layer. Align the CRS, simplify display geometries with simplify(tolerance, preserve_topology=True), and put the smaller layer on the right.

What is the difference between sindex.query and sindex.query_bulk?

query_bulk was the older API for querying the tree with many geometries at once and is deprecated in GeoPandas 1.0. Modern code calls sindex.query with an array of geometries, which returns a two-row array of input positions and tree positions. Both prune candidates via the R-tree; query is the current, vectorised form and the one to use in new code.

Back to Spatial Index Types Reference

Related