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.

Jump to heading# 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.

Jump to heading# 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.

Jump to heading# Step-by-step solution

Jump to heading# 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.

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"

Jump to heading# 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.

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:

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

Jump to heading# 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.

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

Jump to heading# 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:

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.

Jump to heading# 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.

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.

Jump to heading# Verification

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

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.

Jump to heading# 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.

Jump to heading# FAQ

Back to Spatial Index Types Reference

Related

• Spatial Index Types Reference — how the R-tree behind .sindex compares to grid, H3, and GiST indexes
• R-tree vs Grid Index for Point-in-Polygon Queries — when a hand-built grid beats the .sindex R-tree
• CRS & Coordinate Systems Reference — align both layers before joining to keep the index valid
• Query Result Caching — cache the join output so repeated joins never recompute