How to Use the Spatial Index Directly in GeoPandas (sindex)
Problem statement
sjoin covers the common cases. Then you need something it does not do:
- All features within 500 m of each point — not the nearest, all of them.
- Pairs of features that overlap each other within one layer.
- The candidates for a custom test that is not one of the standard predicates.
- A repeated filter inside a loop where rebuilding the join each time is wasteful.
The obvious code is quadratic:
for point in points.geometry:
nearby = parcels[parcels.geometry.distance(point) < 500]That computes 4 million distances per point. On 50,000 points it is 200 billion distance calculations, which is a job that finishes some time next year.
The index is right there and unused:
tree = parcels.sindex # already built, or built once on first accesssindex is the R-tree that sjoin uses internally, and it is a public API. Using it directly is a handful of lines and turns those quadratic loops into near-linear ones.
Quick answer
import geopandas as gpd
import numpy as np
tree = parcels.sindex
# candidates whose bounding box intersects a geometry
idx = tree.query(point)
# candidates that also pass an exact predicate — no second filter needed
idx = tree.query(point, predicate="intersects")
# every pair between two layers, vectorised
left_idx, right_idx = parcels.sindex.query(points.geometry, predicate="contains")
# k nearest, with distances
nearest_idx, dists = tree.nearest(points.geometry, return_distance=True,
max_distance=2_000)
| Call | Returns | Use for |
|---|---|---|
tree.query(geom) |
integer positions of bbox candidates | you will run your own test |
tree.query(geom, predicate="intersects") |
positions passing the exact predicate | a standard predicate |
tree.query(geoseries, predicate=...) |
a (2, n) array of index pairs |
a whole-layer join |
tree.nearest(geoseries) |
nearest positions, optionally distances | nearest-neighbour work |
query returns positional indices, not labels. Use .iloc, or map through gdf.index — this is the single most common mistake.
Step-by-step solution
1. Understand what the index gives you
The R-tree stores bounding boxes, so a query returns candidates whose boxes match — a superset of the true answer:
from shapely.geometry import Point
tree = parcels.sindex
p = Point(383_618, 398_050)
candidates = tree.query(p) # bbox test only
print(f"{len(candidates)} candidates")
exact = [i for i in candidates if parcels.geometry.iloc[i].contains(p)]
print(f"{len(exact)} actually contain the point")7 candidates
1 actually contain the point
Seven boxes overlap the point's box; one polygon actually contains it. That is the two-phase filter described in spatial indexes explained, and the index does the cheap phase.
Passing predicate= makes GeoPandas run the exact phase for you:
exact = tree.query(p, predicate="contains") # already filtered
print(len(exact)) # 1Use the predicate form for standard tests, and the bare form when your test is custom — "within 500 m and of the same class", say.
2. Get the index type right
print(type(parcels.sindex))<class 'geopandas.sindex.SpatialIndex'>
GeoPandas 1.0 uses Shapely 2.0's STRtree underneath. The important properties:
- It is built lazily on first access and cached on the GeoDataFrame.
- It is immutable. Editing geometries after building means the index is stale — GeoPandas invalidates it on assignment, but not if you mutate a geometry in place.
- It is per-object. Any operation returning a new frame produces a frame with no index yet.
import time
t0 = time.perf_counter(); parcels.sindex; build = time.perf_counter() - t0
t0 = time.perf_counter(); parcels.sindex; cached = time.perf_counter() - t0
subset = parcels[parcels["class"] == "residential"]
t0 = time.perf_counter(); subset.sindex; rebuild = time.perf_counter() - t0
print(f"build {build:.3f} s for {len(parcels):,} rows")
print(f"cached {cached:.6f} s")
print(f"rebuild {rebuild:.3f} s for the {len(subset):,}-row subset")build 4.812 s for 4,012,884 rows
cached 0.000002 s
rebuild 3.104 s for the 2,841,002-row subset
Nearly five seconds to build over 4 million rows. Inside a loop that filters and queries, that cost repeats every iteration — which is why building once outside the loop is often the whole optimisation.
3. Query positionally, and convert carefully
query returns positional indices into the GeoDataFrame, which are not the same as label indices unless the index happens to be a default RangeIndex:
parcels = parcels.set_index("parcel_id") # now labels are strings
idx = tree.query(p, predicate="contains")
print(idx) # array([412881]) ← a position
# ❌ this looks up a label, and will raise or return the wrong row
parcels.loc[idx]
# ✅ positional
parcels.iloc[idx]
# ✅ or convert to labels explicitly
labels = parcels.index[idx]
parcels.loc[labels]This is the mistake that produces silently wrong answers rather than an error, because on a default RangeIndex .loc and .iloc agree — so code works during development and breaks the moment someone sets a meaningful index.
4. Use the vectorised form for whole-layer work
Passing a GeoSeries queries every geometry at once, in C:
left_idx, right_idx = parcels.sindex.query(points.geometry, predicate="contains")
print(left_idx[:5], right_idx[:5])[0 0 1 2 2] [412881 412902 88104 4 918]
The result is two aligned arrays of positions: left_idx[k] is a position in the queried series (points) and right_idx[k] is a position in the indexed frame (parcels). A point matching two parcels appears twice.
Turn that into a joined frame:
import pandas as pd
pairs = pd.DataFrame({
"point_id": points.index[left_idx],
"parcel_id": parcels.index[right_idx],
})Note the argument order, which catches everyone once: the indexed frame is the one whose sindex you called, and the queried geometries are the argument. parcels.sindex.query(points.geometry, predicate="contains") asks "which parcels contain each point", so predicate describes the relationship of the indexed geometry to the queried one.
5. Handle distance queries
The R-tree indexes boxes, so a distance query needs the box expanded first:
def within_distance(tree, gdf, geom, distance):
"""All features within `distance` of geom, using the index."""
from shapely.geometry import box
minx, miny, maxx, maxy = geom.bounds
search = box(minx - distance, miny - distance,
maxx + distance, maxy + distance)
candidates = tree.query(search) # cheap: bbox only
sub = gdf.iloc[candidates]
return sub[sub.geometry.distance(geom) <= distance] # exact: on candidates only%timeit parcels[parcels.geometry.distance(p) < 500] # 3.44 s
%timeit within_distance(tree, parcels, p, 500) # 0.0021 s1,600× faster, and the answers are identical. The pattern is the general one: expand the box, query the index, run the exact test on the candidates.
GeoPandas also exposes it directly:
idx = tree.query(p.buffer(500), predicate="intersects")Buffering the query geometry is simpler and slightly more expensive, since it builds a real circular polygon rather than a rectangle. For a single query the difference is irrelevant; for millions it is not.
6. Use nearest for nearest-neighbour work
nearest_idx, dists = parcels.sindex.nearest(
points.geometry, return_distance=True, max_distance=2_000, return_all=False)
points["nearest_parcel"] = parcels.index[nearest_idx[1]]
points["distance_m"] = distsnearest returns a (2, n) array: row 0 is positions in the queried series, row 1 is positions in the indexed frame. return_all=False keeps one result per input even when distances tie exactly — otherwise ties produce extra rows, which is the same behaviour that makes sjoin_nearest return more rows than it received.
max_distance matters more than it looks. Without it, a query point far from everything walks a long way through the tree; with it, the search is bounded and points with no neighbour within the limit are simply omitted.
Code examples
Example 1: finding overlapping pairs within one layer
A self-join is where the index earns the most, because the naive form is O(n²):
import numpy as np
import pandas as pd
import geopandas as gpd
def find_overlaps(gdf, *, min_area=0.0):
"""Every pair of features in `gdf` whose geometries overlap."""
tree = gdf.sindex
left, right = tree.query(gdf.geometry, predicate="intersects")
# drop self-matches and keep each unordered pair once
keep = left < right
left, right = left[keep], right[keep]
a = gdf.geometry.values[left]
b = gdf.geometry.values[right]
import shapely
inter = shapely.intersection(a, b)
areas = shapely.area(inter)
# touching neighbours intersect along a line, with zero area
real = areas > min_area
return pd.DataFrame({
"a": gdf.index[left[real]],
"b": gdf.index[right[real]],
"overlap_m2": areas[real],
}).sort_values("overlap_m2", ascending=False)
overlaps = find_overlaps(gpd.read_file("parcels.gpkg"), min_area=0.01)
print(f"{len(overlaps):,} overlapping pairs, "
f"{overlaps['overlap_m2'].sum():,.0f} m² total")
print(overlaps.head())1,204 overlapping pairs, 88,412 m² total
a b overlap_m2
882 41882004 41882091 8,204.1
4 10004112 10004118 4,118.9
left < right does two things at once: it removes self-matches (where left == right) and it keeps each pair in one direction only, so a set of 1,204 real overlaps does not report as 2,408.
The min_area filter is essential for a polygon coverage. Adjacent parcels legitimately share a boundary, so they intersect along a line whose area is zero. Without the filter, every neighbouring pair in the layer is reported as an overlap. Using this for data quality is covered in how to fix gaps and overlaps in a polygon coverage.
shapely.intersection on two arrays is the vectorised pairwise form, so the exact phase runs in C rather than in a Python loop over candidate pairs.
Example 2: a reusable neighbourhood query
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
class Neighbourhood:
"""Repeated radius queries against one layer, with the index built once."""
def __init__(self, gdf, geometry_col=None):
self.gdf = gdf
self.geoms = gdf.geometry.values
self.tree = gdf.sindex # built here, once
self.crs = gdf.crs
def within(self, geom, distance, *, columns=None):
"""Features within `distance` of one geometry."""
minx, miny, maxx, maxy = geom.bounds
search = shapely.box(minx - distance, miny - distance,
maxx + distance, maxy + distance)
cand = self.tree.query(search)
if len(cand) == 0:
return self.gdf.iloc[[]]
d = shapely.distance(self.geoms[cand], geom)
hit = cand[d <= distance]
out = self.gdf.iloc[hit]
out = out[columns] if columns else out
return out.assign(distance_m=d[d <= distance])
def within_many(self, geoms, distance):
"""Vectorised: every (query, match) pair within `distance`."""
buffered = shapely.buffer(np.asarray(geoms), distance, quad_segs=2)
left, right = self.tree.query(buffered, predicate="intersects")
exact = shapely.distance(np.asarray(geoms)[left], self.geoms[right])
keep = exact <= distance
return pd.DataFrame({
"query_pos": left[keep],
"match_index": self.gdf.index[right[keep]],
"distance_m": exact[keep],
})
def counts_within(self, geoms, distance):
"""How many features are within `distance` of each query geometry."""
pairs = self.within_many(geoms, distance)
counts = pairs.groupby("query_pos").size()
return counts.reindex(range(len(geoms)), fill_value=0).to_numpy()
hood = Neighbourhood(gpd.read_file("amenities.gpkg").to_crs(27700))
homes = gpd.read_file("homes.gpkg").to_crs(27700)
homes["amenities_500m"] = hood.counts_within(homes.geometry, 500)
print(homes["amenities_500m"].describe())count 184204.000000
mean 12.408000
std 9.882000
min 0.000000
max 88.000000
Building the index once in __init__ is the point of the class. A loop that calls within 200,000 times pays the build cost once rather than 200,000 times.
quad_segs=2 on the buffer is a deliberate approximation: the buffered geometry is only used to query the index, and the exact distance test runs afterwards. A coarser circle means fewer vertices and a faster bbox computation, with no effect on the answer.
reindex(range(len(geoms)), fill_value=0) restores queries that matched nothing. Without it, groupby silently omits them and the result array is shorter than the input — a misalignment that assigns counts to the wrong rows.
Example 3: a custom join sjoin cannot express
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
def join_with_rule(left, right, *, max_distance, rule, right_cols):
"""Nearest match within max_distance that also satisfies `rule(l_row, r_row)`."""
tree = right.sindex
buffered = shapely.buffer(left.geometry.values, max_distance, quad_segs=2)
l_idx, r_idx = tree.query(buffered, predicate="intersects")
dist = shapely.distance(left.geometry.values[l_idx], right.geometry.values[r_idx])
order = np.lexsort((dist, l_idx)) # by left, then by distance
l_idx, r_idx, dist = l_idx[order], r_idx[order], dist[order]
chosen = {}
for li, ri, d in zip(l_idx, r_idx, dist):
if li in chosen:
continue # already found a nearer valid match
if d > max_distance:
continue
if rule(left.iloc[li], right.iloc[ri]):
chosen[li] = (ri, d)
out = left.copy()
for col in right_cols:
out[col] = [right.iloc[chosen[i]][col] if i in chosen else None
for i in range(len(left))]
out["match_distance_m"] = [chosen[i][1] if i in chosen else np.nan
for i in range(len(left))]
print(f"{len(chosen):,} of {len(left):,} matched")
return out
# "the nearest school of the right phase, within 2 km"
result = join_with_rule(
homes, schools,
max_distance=2_000,
rule=lambda home, school: school["phase"] == home["required_phase"],
right_cols=["school_name", "phase"],
)178,204 of 184,204 matched
The Python loop here is not a mistake — it exists because rule is arbitrary user code that cannot be vectorised. What matters is how much it runs: the index has already reduced the candidate set from 184,204 × 4,102 pairs to about 2.1 million, and lexsort puts the nearest candidate first for each left row, so the loop breaks on the first acceptable match and touches only a fraction of even that.
np.lexsort((dist, l_idx)) sorts by the last key first, which is the opposite of what most people expect — it gives groups ordered by l_idx with each group sorted by dist.
This is the general shape of a custom spatial join: use the index to get from quadratic to a manageable candidate set, then apply whatever logic the problem actually requires.
Explanation
sindex is a public interface to the same R-tree that sjoin uses internally, and reaching for it directly is worth it whenever the built-in joins do not express what you need.
The R-tree is a hierarchy of bounding boxes. Leaf nodes hold one box per geometry; internal nodes hold boxes enclosing their children. A query descends the tree, discarding any subtree whose box does not match — so a query touches a number of nodes proportional to the logarithm of the collection size rather than to its size. Since Shapely 2.0 this is an STRtree, packed with the Sort-Tile-Recursive algorithm, which builds a well-balanced tree in one pass and is why building over 4 million geometries takes seconds rather than minutes.
The index only answers box questions, which is why every use follows the same two-phase shape: query the tree for candidates, then run the exact test on them. predicate= performs the second phase for you; omitting it hands you the candidates so you can apply your own test. Neither is more correct — the choice depends on whether your test is one GEOS provides.
Positional versus label indexing is the trap, and it is worth being deliberate about because it fails silently. query returns positions, because the tree knows nothing about pandas labels. On a default RangeIndex positions and labels coincide, so .loc[idx] works — and then someone calls set_index("parcel_id") and the same line starts returning wrong rows or raising KeyError. Always .iloc, or convert explicitly with gdf.index[idx].
Distance queries need the box expanded, not the tree changed. An R-tree cannot answer "within 500 m" directly, because it indexes boxes and 500 m is not a box. Expanding the query geometry's box by the distance produces a region guaranteed to contain every feature within that distance — a superset, since a corner of the expanded box is more than 500 m from the original. The exact test then removes the extras. That is why the pattern is always expand-query-filter, and why the buffer used for the query can be a crude approximation.
The index's per-object caching has a practical consequence that catches people in loops. Any operation returning a new GeoDataFrame — a filter, a copy, a to_crs — produces an object with no index built. Code that filters and then queries inside a loop rebuilds the tree every iteration, which can cost more than the loop saves. Build the index on a stable object outside the loop and query it, as the class in Example 2 does.
Finally, the honest boundary: use sjoin when it fits. It is well tested, handles the join semantics and index alignment correctly, and returns a properly formed frame. Reach for sindex when you need something it does not express — a custom predicate, a radius query, a self-join, or a candidate set you will process yourself. The techniques here are for the gap, not a replacement for the tool.
Edge cases or notes
queryreturns positional indices. Use.iloc, orgdf.index[idx]to get labels.- The index is built lazily and cached per object. A filter or
copyproduces a frame with no index yet. - Editing a geometry in place does not invalidate the index. Assign a new GeoSeries instead.
predicate=runs the exact test; omitting it returns bbox candidates only.- The vectorised
query(geoseries)returns(2, n)positions, indexed frame second. nearest(..., return_all=False)keeps one result per input when distances tie.max_distanceonnearestbounds the search and omits inputs with no neighbour inside it.sindex.queryon an empty frame returns empty arrays rather than raising, so guard withlen(...).- Predicate cost scales with vertex count. A few enormous geometries dominate the exact phase; consider simplifying.
shapely.STRtreeis available directly if you want a tree over an array of geometries with no GeoDataFrame.
Internal links
- Spatial indexes explained: R-trees and why spatial joins are fast — the mechanism
- How to perform a spatial join in Python (GeoPandas) — use this when it fits
- Why GeoPandas is slow: the four real bottlenecks — the missing index as one of four
- How to replace a row loop with vectorised GeoPandas operations — the other half of most speedups
- How to find the nearest point in GeoPandas — the common nearest-neighbour case
- Nearest-neighbour joins explained — ties, radius and cardinality
- How to fix gaps and overlaps in a polygon coverage — the self-join in Example 1, applied
- PostGIS spatial indexes explained — the same idea, on disk
FAQ
What does sindex.query return?
Positional indices into the GeoDataFrame, not label indices. Use .iloc[idx], or convert with gdf.index[idx] — on a non-default index, .loc[idx] gives wrong rows or raises.
When should I use sindex instead of sjoin?
When you need a custom predicate, a radius query, a self-join, or the candidate set itself. sjoin handles the standard cases better than hand-written code.
How do I find everything within a distance?
Expand the query geometry's bounding box by the distance, query the index for candidates, then run an exact distance test on those. Roughly 1,600× faster than testing every feature.
Why does my query return more features than expected?
The index tests bounding boxes, which is a superset of the true answer. Pass predicate= to run the exact test, or filter the candidates yourself.
Is the index rebuilt every time I use it?
No — it is built on first access and cached on that GeoDataFrame. But any operation returning a new frame produces one with no index, so building inside a loop is a real cost.
How do I find overlapping features within one layer?
Query the layer's index with its own geometries, keep pairs where left < right to remove self-matches and duplicates, then filter by intersection area to exclude neighbours that merely touch.
Does sindex work on points, lines and polygons alike?
Yes — it indexes bounding boxes, which every geometry type has. A point's box is degenerate, and the tree handles that fine.