Spatial Index Types Reference for Python Geospatial Dashboards

Every interactive spatial dashboard eventually hits the same wall: a filter changes, and the app has to decide which of a million geometries fall inside a viewport, a drawn polygon, or a buffer. Done naively, that decision is an O(n) scan that tests the query shape against every feature in the layer — and it runs again on every widget interaction. A spatial index replaces that linear scan with a tree or bucket lookup that discards almost all of the layer before a single expensive geometric predicate executes. Choosing the right index type is the difference between a dashboard that answers in 20 milliseconds and one that freezes for several seconds on each pan.

This reference catalogues the spatial index types you will encounter across the Python geospatial stack — the R-tree behind GeoPandas .sindex, the packed STRtree in Shapely, uniform grid and geohash bucket indexes, H3 hexagonal indexing, quadtrees, the KD-tree used for nearest-neighbour search, and the PostGIS GiST index that pushes filtering down into the database tier. It sits under Spatial Data Reference alongside the CRS & Coordinate Systems Reference, because an index is only correct when both layers share a coordinate system. A comparison matrix, a runnable .sindex workflow, and a set of production troubleshooting notes follow.

Jump to heading# Why a spatial index matters

A point-in-polygon test is not free. Shapely evaluates the ray-casting or winding-number algorithm across every vertex of a polygon ring, so testing one point against one 500-vertex census tract is hundreds of coordinate comparisons. A brute-force join of 100,000 GPS pings against 2,000 districts is 200 million such tests — the classic O(n·m) blow-up.

A spatial index short-circuits almost all of that work with a cheap pre-filter. Every geometry has a minimum bounding rectangle (MBR): the axis-aligned box that fully contains it. Comparing two boxes is four float comparisons. An index organises those boxes so that, given a query box, it returns only the handful of geometries whose MBRs could overlap — the candidate set. The exact, expensive predicate then runs against candidates only. For a query that touches ten districts out of two thousand, the index discards 99.5% of the polygons before Shapely ever runs.

Not every index solves the same problem. R-trees excel at range and overlap queries against irregular polygons. KD-trees answer nearest-neighbour questions on point clouds. H3 turns aggregation into a group-by. The matrix below is the decision surface.

Jump to heading# Prerequisites

• Python 3.10+ with geopandas>=0.14 (the .sindex API and predicate= keyword are stable from 0.10; the array-input form of sindex.query shown here needs geopandas>=1.0).
• shapely>=2.0 — provides the vectorised STRtree and releases the GIL during predicate evaluation.
• rtree>=1.0 (bundling libspatialindex) as the default .sindex backend, or pygeos-era builds where Shapely 2 supplies the tree directly.
• Optional: h3>=4.0 for hexagonal indexing, scipy>=1.10 for cKDTree, and a PostGIS-enabled PostgreSQL for the database-tier examples.
• Familiarity with EPSG codes and reprojection — an index built in geographic degrees (EPSG:4326) measures distance incorrectly, which is covered in the CRS & Coordinate Systems Reference.

Install the core stack:

pip install "geopandas>=1.0" "shapely>=2.0" rtree h3 scipy

Jump to heading# Spatial index comparison matrix

The table below is the reference dataset for this page. “Query complexity” is the average cost of a single lookup once the index is built; skewed data degrades grid and geohash indexes toward the linear worst case.

Two rules of thumb fall out of this table. First, for exact geometric predicates against irregular polygons, an R-tree family index (.sindex, STRtree, or GiST) is almost always the right default. Second, a grid, geohash, or H3 index wins when you can reframe the question as bucketing rather than exact overlap — and that reframing usually simplifies the downstream aggregation as much as it speeds the lookup.

Jump to heading# Index families beyond the R-tree

The R-tree is the default, but three other families solve problems it handles poorly, and knowing when to reach for them is what separates a fast dashboard from a slow one.

Geohash encodes a longitude/latitude pair as a base-32 string by recursively bisecting the world into quadrants — each added character narrows the cell roughly eightfold. Its defining property is that nearby locations usually share a string prefix, so a prefix range scan in a plain key-value store (Redis, a sorted database column) retrieves a neighbourhood without any geometry library at all. That makes geohash ideal for coarse proximity lookups and for sharding spatial data across a distributed store. Its weakness is the edge case: two points on opposite sides of a cell boundary can be metres apart yet share no prefix, so geohash proximity is an approximation, not a guarantee.

import pygeohash as pgh

# Central Berlin at precision 7 ≈ 150 m cells
gh = pgh.encode(52.520, 13.404, precision=7)     # 'u33dc0d'
neighbours = pgh.geohash_approximate_distance(gh, pgh.encode(52.521, 13.405, 7))

Quadtree subdivides space into four quadrants recursively, but only where data actually exists — dense regions split deeper than empty ones. This adaptivity is its advantage over a uniform grid: a metropolitan hotspot gets fine cells while surrounding countryside stays coarse, so lookup cost stays roughly balanced regardless of skew. Quadtrees underpin level-of-detail rendering, where you draw coarse cells when zoomed out and descend the tree as the user zooms in. They are a natural fit when point density varies by orders of magnitude across the extent and a single grid cell size cannot serve both regimes.

KD-tree answers a different question entirely: not “which polygon contains this point?” but “what are the k nearest points to this one?” It partitions a point set with alternating axis-aligned splits, giving logarithmic nearest-neighbour queries. For a dashboard that snaps a click to the closest sensor, finds the five nearest facilities, or builds a distance-weighted interpolation, scipy.spatial.cKDTree is the right tool — an R-tree can approximate this with expanding bounding-box queries, but a KD-tree does it directly and faster.

import numpy as np
import geopandas as gpd
from scipy.spatial import cKDTree

# Sensor GeoDataFrame projected to EPSG:32633 metres
sensors_gdf = gpd.read_file("sensors.gpkg").to_crs("EPSG:32633")
coords = np.column_stack(
    [sensors_gdf.geometry.x.values, sensors_gdf.geometry.y.values]
)                                            # shape (n, 2)
tree = cKDTree(coords)

# Five nearest sensors to a clicked location, with distances in metres
click = np.array([[391000.0, 5819000.0]])
dist, pos = tree.query(click, k=5)
print("nearest sensor rows:", pos[0], "distances (m):", dist[0].round(1))

Use a KD-tree for point-to-point proximity, a quadtree for adaptive density and level-of-detail, geohash for prefix-based sharding, and the R-tree family for exact polygon predicates. Most production dashboards end up combining two — an R-tree for the join and a KD-tree for the “nearest” widget.

Jump to heading# Core implementation workflow

This walkthrough builds the GeoPandas R-tree once and reuses it for both a bounding-box query and a spatial join. The scenario: 120,000 GPS pings joined to Berlin administrative districts. Both layers are reprojected to EPSG:32633 (UTM zone 33N) so the index operates on metres, not degrees.

Jump to heading# Step 1 — Load geometries and align the CRS

import geopandas as gpd
import numpy as np
from shapely.geometry import Point

# Districts (polygons) — loaded from a local GeoPackage in WGS 84
districts = gpd.read_file("berlin_districts.gpkg").to_crs("EPSG:32633")

# Synthetic GPS pings around central Berlin (lon 13.40, lat 52.52)
rng = np.random.default_rng(42)
lon = rng.normal(13.404, 0.06, size=120_000)
lat = rng.normal(52.520, 0.04, size=120_000)
pings = gpd.GeoDataFrame(
    {"ping_id": np.arange(120_000)},
    geometry=gpd.points_from_xy(lon, lat, crs="EPSG:4326"),
).to_crs("EPSG:32633")

assert pings.crs == districts.crs, "CRS mismatch will corrupt every index lookup"

Aligning the CRS before indexing is not optional. An R-tree built over one layer’s bounding boxes is meaningless against another layer expressed in a different coordinate system.

Jump to heading# Step 2 — Trigger the R-tree index

GeoPandas builds .sindex lazily. Touch it once and the libspatialindex R-tree is constructed over the geometry column’s bounding boxes and cached on the frame.

sindex = districts.sindex          # builds the R-tree over district MBRs
print(type(sindex))                # <class 'geopandas.sindex.SpatialIndex'>
print("indexed geometries:", len(districts))

The tree now lives on districts until its geometry column is reassigned. Every later sjoin, clip, or .cx slice against this frame reuses it.

Jump to heading# Step 3 — Run a bounding-box query

To find every district that intersects a viewport rectangle, pass a single query geometry to sindex.query. It returns integer positions into districts, not the geometries themselves.

from shapely.geometry import box

# Viewport in EPSG:32633 metres (roughly central Berlin)
viewport = box(389000, 5817000, 393000, 5821000)

candidate_pos = districts.sindex.query(viewport, predicate="intersects")
visible = districts.iloc[candidate_pos]
print(f"{len(visible)} of {len(districts)} districts intersect the viewport")

The predicate="intersects" argument makes the index run the exact predicate on candidates internally, so visible is already precise — no second filtering pass is required. Omit the predicate to get the looser bounding-box-only candidate set when you want to run a custom test yourself.

Jump to heading# Step 4 — Run an indexed spatial join

gpd.sjoin uses the .sindex R-tree of the right-hand frame internally. You do not call the index directly; you pick the predicate that matches your question. For “which district is each ping in?”, the predicate is within.

joined = gpd.sjoin(
    pings,
    districts[["district_name", "geometry"]],
    how="left",
    predicate="within",
)
counts = joined.groupby("district_name")["ping_id"].count()
print(counts.sort_values(ascending=False).head())

Internally, sjoin queries the R-tree with every ping’s bounding box (a degenerate point box), collects candidate district positions, then runs the exact within test only on those pairs. On this dataset it evaluates a few hundred thousand exact predicates instead of 120,000 × N districts.

Jump to heading# Step 5 — Verify candidate reduction and correctness

Prove the index actually pruned the search by comparing against a brute-force scan on a small sample and by measuring the candidate count.

import time

sample = pings.iloc[:2000]

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

t0 = time.perf_counter()
brute = sample.assign(
    district_name=sample.geometry.apply(
        lambda p: next((n for n, g in zip(districts.district_name, districts.geometry)
                        if g.contains(p)), None)
    )
)
brute_ms = (time.perf_counter() - t0) * 1000

matched = (idx_join["district_name"].values == brute["district_name"].values).mean()
print(f"indexed={idx_ms:.0f}ms  brute={brute_ms:.0f}ms  agreement={matched:.1%}")
# Expected: indexed far faster than brute; agreement ~100%

Agreement should be effectively 100% (modulo points exactly on a shared boundary — see Troubleshooting), and the indexed path should be several times faster even on this small sample. The gap widens dramatically as the layer grows.

Jump to heading# Advanced patterns

Jump to heading# STRtree for static reference layers

When a polygon layer never changes during the session — administrative boundaries, service areas, election precincts — a Shapely STRtree is faster to query than the general-purpose .sindex because it is bulk-loaded and immutable. Build it once at startup and cache it as a Streamlit resource so every session shares one tree.

import streamlit as st
from shapely import STRtree

@st.cache_resource
def build_district_tree(gpkg_path: str):
    gdf = gpd.read_file(gpkg_path).to_crs("EPSG:32633")
    tree = STRtree(gdf.geometry.values)          # packed, read-only
    return tree, gdf

tree, districts = build_district_tree("berlin_districts.gpkg")

# Query returns positions into the input geometry array
hit_positions = tree.query(viewport, predicate="intersects")
visible = districts.iloc[hit_positions]

Because the tree is immutable, caching it with @st.cache_resource is safe across the query result caching layer — there is no per-user mutation to corrupt.

Jump to heading# H3 hexagonal indexing for aggregation

When the goal is how many pings per area, not which exact polygon, skip geometry entirely. Assign each point an H3 cell id at a chosen resolution and group by the string. This converts a spatial join into a hash aggregation.

import h3
import pandas as pd

# Resolution 8 ≈ 0.7 km² hexagons; use lon/lat in EPSG:4326
pings_wgs = pings.to_crs("EPSG:4326")
cells = [
    h3.latlng_to_cell(y, x, 8)
    for x, y in zip(pings_wgs.geometry.x, pings_wgs.geometry.y)
]
density = pd.Series(cells).value_counts()
print(density.head())

# Neighbouring cells for a smoothing kernel are a cheap grid walk
ring = h3.grid_disk(density.index[0], 1)   # centre cell + 6 neighbours

H3 shines for choropleth binning and hotspot detection because equal-area hexagons remove the visual bias of irregular administrative polygons, and grid_disk gives neighbour traversal without any geometry math.

Jump to heading# PostGIS GiST pushed down via SQL

For layers too large to hold in the dashboard process, keep the geometry in PostGIS and let the GiST index do the filtering server-side. The && bounding-box operator hits the index; ST_Within refines the candidates — all before a single row crosses the wire.

CREATE INDEX IF NOT EXISTS pings_geom_gist
    ON pings USING GIST (geom);

-- The && operator uses the GiST index; ST_Within refines candidates
SELECT d.district_name, count(*) AS ping_count
FROM pings p
JOIN districts d
  ON p.geom && d.geom            -- index-accelerated bbox pre-filter
 AND ST_Within(p.geom, d.geom)   -- exact predicate on survivors
GROUP BY d.district_name
ORDER BY ping_count DESC;

Run EXPLAIN ANALYZE on this query and confirm the plan shows an Index Scan on pings_geom_gist rather than a Seq Scan. Pushing the filter into the database keeps the dashboard from ever materialising the full ping table in Python.

Jump to heading# Verification and testing

Beyond correctness, assert that the index is doing its job. A silent fall-back to a linear scan looks correct but performs terribly under load.

# 1. Confirm the index exists and covers every row
assert districts.has_sindex, "sindex was never built"
assert districts.sindex.size == len(districts)

# 2. Confirm the bbox query returns a strict subset, not the whole layer
cand = districts.sindex.query(viewport)          # no predicate = bbox candidates
assert len(cand) < len(districts), "index returned every row — check CRS/extent"

# 3. Confirm an indexed join and a manual bbox+contains filter agree in count
manual = pings[pings.geometry.within(districts.geometry.union_all())]
indexed = gpd.sjoin(pings, districts, predicate="within", how="inner")
assert abs(len(manual) - len(indexed)) / len(indexed) < 0.001

For regression protection in CI, snapshot the join row count and the per-district totals. A schema change upstream — a district polygon re-digitised, a CRS drift on ingest — shows up immediately as a count that no longer matches the snapshot.

Jump to heading# Troubleshooting

sjoin is slow and EXPLAIN-style profiling shows a full scan : The index was never triggered, usually because you iterated geometries in a Python loop with .contains() instead of calling sjoin, .cx, or sindex.query. Manual loops bypass .sindex entirely. Route the filter through gpd.sjoin(points, polys, predicate="within") or polys.sindex.query(...) so the R-tree is actually consulted.

Query returns wrong candidates after editing geometries in place : The cached R-tree still describes the old bounding boxes. GeoPandas invalidates .sindex when you assign a new geometry column, but mutating existing geometries in place (for example via a Shapely operation that edits vertices) leaves the stale tree in place. Rebuild by reassigning the column: gdf = gdf.set_geometry(gdf.geometry.simplify(5)) forces a fresh index on next access.

TypeError: 'op' is an invalid keyword argument for sjoin : The op= parameter was renamed to predicate= in GeoPandas 0.10. Update calls to gpd.sjoin(left, right, predicate="intersects"). If you must support very old versions, gate on geopandas.__version__.

Every point matches every polygon, or none do : The two layers are in different coordinate systems, so their bounding boxes never (or always) overlap. Assert left.crs == right.crs before the join and reproject the smaller layer with .to_crs(...). This is the single most common .sindex bug — the index is correct, the coordinates are not.

Points on a shared district boundary appear in two joins or none : within is strict — a point exactly on the shared edge of two polygons may satisfy neither within or both intersects. For deterministic assignment, snap coordinates to a tolerance before the join, or choose the predicate deliberately: intersects includes boundaries, within excludes them. Document which convention the dashboard uses.

Jump to heading# Performance considerations

Build cost is amortised, not free. The first .sindex access on a large layer can take hundreds of milliseconds. In a dashboard, build reference-layer indexes once at startup behind @st.cache_resource rather than on every rerun. Point layers that change with user filters should be indexed after the filter, since re-indexing a 10,000-row subset is far cheaper than querying a stale million-row tree.

Match the index to data density. Uniform grids degrade to O(n) when points concentrate into a few cells — a national dataset with 90% of its rows in three metro areas will overload those buckets. R-trees and quadtrees adapt to density automatically, which is why .sindex is the safer default for real-world skewed geospatial data.

Reduce candidates before the exact test. The exact predicate is the expensive part. Simplify display-only polygons with gdf.geometry.simplify(tolerance, preserve_topology=True) before building the index so each MBR check and each contains call touches fewer vertices. Pair this with the Query Result Caching layer so repeated joins over the same viewport never recompute.

Prefer bulk queries. Calling sindex.query once per point in a Python loop reintroduces the overhead the index was meant to remove. Pass an array of geometries so the query vectorises in C, or use gpd.sjoin, which does this for you.

Jump to heading# Frequently asked questions

Back to Spatial Data Reference

Related

• R-tree vs Grid Index for Point-in-Polygon Queries — head-to-head benchmark of the two dominant point-in-polygon index strategies
• Speeding Up GeoPandas Spatial Joins with .sindex — practical tuning of the R-tree behind gpd.sjoin
• CRS & Coordinate Systems Reference — why an index is only valid when both layers share a coordinate system
• Query Result Caching — cache indexed join results so repeated viewport queries never recompute