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]

On the 293,729 Land Registry parcels of Leeds used on this page, that computes 293,729 distances per point. For the 9,227 crimes West Yorkshire Police recorded in Leeds in June 2026 it is 2.7 billion distance calculations — about half an hour at the 0.18 s per point measured below, for a question the index answers in a few seconds.

The index is right there and unused:

tree = parcels.sindex        # already built, or built once on first access

sindex 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="within")   # "each point within which parcel"

# k nearest, with distances
nearest_idx, dists = tree.nearest(points.geometry, return_distance=True,
                                  max_distance=2_000)

query returns positional indices, not labels. Use .iloc, or map through gdf.index — this is the single most common mistake.

2 candidates
1 actually contain the point

What the index returns, positional queries, bulk and nearest

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(429_700, 433_900)                      # Park Square, in central Leeds

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")

2 candidates
1 actually contain the point

Two parcels' boxes cover the point; only the polygon that actually contains it survives the exact test. 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="within")        # already filtered: p within the parcel
print(exact)

[48737]

Use 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

fresh = parcels.copy()                            # a new frame: no index built yet
t0 = time.perf_counter(); fresh.sindex; build = time.perf_counter() - t0
t0 = time.perf_counter(); fresh.sindex; cached = time.perf_counter() - t0
subset = fresh[fresh.geometry.area < 500]         # mostly house plots
t0 = time.perf_counter(); subset.sindex; rebuild = time.perf_counter() - t0

print(f"build   {build:.3f} s for {len(fresh):,} rows")
print(f"cached  {cached:.6f} s")
print(f"rebuild {rebuild:.3f} s for the {len(subset):,}-row subset")

build   0.087 s for 293,729 rows
cached  0.000165 s
rebuild 0.073 s for the 246,282-row subset

A quarter of a second over the whole city is cheap once. Inside a loop that filters and queries, that cost repeats every iteration — a few hundred iterations and the rebuilds cost more than everything else — 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("INSPIREID")          # now the labels are title-plan ids

idx = tree.query(p, predicate="within")
print("position:", idx)

# ❌ this looks up a label, and will raise or return the wrong row
try:
    parcels.loc[idx]
except KeyError as err:
    print("KeyError:", str(err)[:60])

# ✅ positional
print("iloc:", parcels.iloc[idx].index.tolist())

# ✅ or convert to labels explicitly
labels = parcels.index[idx]
print("loc: ", parcels.loc[labels].index.tolist())

position: [48737]
KeyError: "None of [Index([48737], dtype='int64', name='INSPIREID')] a
iloc: [26535072]
loc:  [26535072]

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="within")
print(left_idx[:5], right_idx[:5])

hits = np.bincount(left_idx, minlength=len(points))
print(f"{len(left_idx):,} pairs; {(hits == 0).sum():,} crimes in no parcel, "
      f"{(hits == 1).sum():,} in one, {(hits > 1).sum():,} in more than one")

[ 0  1  2  5 10] [241032  96039 139117 112597  96484]
6,329 pairs; 2,983 crimes in no parcel, 6,165 in one, 79 in more than one

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 — and on real data that happens, because INSPIRE polygons overlap. Most of the crimes in no parcel are on the street: police.uk snaps every crime to a nearby anonymised point, usually on a road, and most of the highway is not a registered title.

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. The predicate reads from the queried geometry to the indexed one: parcels.sindex.query(points.geometry, predicate="within") asks "is each point within a parcel", which is the question here. Shapely's own documentation puts it the same way — the predicate is evaluated as predicate(input_geometry, tree_geometry).

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

t0 = time.perf_counter(); brute = parcels[parcels.geometry.distance(p) < 500]; t_brute = time.perf_counter() - t0
t0 = time.perf_counter(); fast = within_distance(tree, parcels, p, 500); t_fast = time.perf_counter() - t0
print(f"every parcel {t_brute:.3f} s · index {t_fast:.4f} s · {t_brute / t_fast:.0f}x · "
      f"{len(fast):,} parcels, same answer: {brute.index.sort_values().equals(fast.index.sort_values())}")

every parcel 0.116 s · index 0.0012 s · 93x · 623 parcels, same answer: True

Much 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 does all three steps in one call:

idx = tree.query(p, predicate="dwithin", distance=500)
print(f"{len(idx):,} parcels, same as within_distance: {len(idx) == len(fast)}")

623 parcels, same as within_distance: True

Reach for dwithin rather than buffering the query geometry and testing intersects. A buffer is a polygon inside the true circle — its vertices are on the circle and its edges cut across — so a feature just short of the distance, between two vertices, can fall outside it. Example 2 measures how much that loses.

6. Use nearest for nearest-neighbour work

nearest_idx, dists = parcels.sindex.nearest(
    points.geometry, return_distance=True, max_distance=50, return_all=False)

# row 0 says which point each result belongs to — use it, because points with
# nothing inside max_distance are left out and the arrays are shorter than `points`
found = points.index[nearest_idx[0]]
points["nearest_parcel"] = pd.Series(parcels.index[nearest_idx[1]], index=found)
points["distance_m"] = pd.Series(dists, index=found)
print(f"{len(found):,} of {len(points):,} crimes have a parcel within 50 m; "
      f"median distance {points['distance_m'].median():.1f} m")

9,213 of 9,227 crimes have a parcel within 50 m; median distance 0.0 m

nearest returns a (2, n) array: row 0 is positions in the queried series, row 1 is positions in the indexed frame. Assigning row 1 straight into a column — the obvious line — raises a length mismatch as soon as max_distance drops a single point, and a later refactor that drops the same number from somewhere else would make it silently misalign instead. 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.

Overlap pairs and custom joins on Leeds parcels

The data used on this page

Everything here runs on Leeds. parcels is HM Land Registry's INSPIRE Index Polygons for the city — the extent of every registered freehold title, 293,729 polygons. points is every street-level crime West Yorkshire Police published for June 2026 inside the Leeds district boundary, from data.police.uk. homes is the 2021 census output areas of Leeds as centroids, from the ONS, and stops is every active bus stop in West Yorkshire from the Department for Transport's NaPTAN. The parcels contain HM Land Registry data © Crown copyright and database right 2026, licensed under the Open Government Licence v3.0.

import io
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import pyogrio
import requests
import shapely

# HM Land Registry INSPIRE Index Polygons: the extent of every registered freehold in Leeds.
# The download sets a cookie and redirects to itself, so it needs a session to keep the cookie.
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})    # otherwise GDAL spends minutes fetching the XSD
parcels = gpd.read_file("zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"])

# the Leeds district boundary, to keep only crimes inside it
ONS = "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services"
leeds = gpd.read_file(f"{ONS}/Local_Authority_Districts_DEC_2025_Boundaries_UK_BGC/FeatureServer/0/query"
                      "?where=LAD25NM%3D%27Leeds%27&outFields=LAD25CD&outSR=27700&f=geojson").set_crs(27700, allow_override=True)

# every street-level crime West Yorkshire Police published for June 2026, as points
west, south, east, north = leeds.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            r = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                              data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if r.status_code != 429:
                break
            time.sleep(2 ** attempt)
        r.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in r.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id").sort_values("id")   # the API's order varies
points = gpd.GeoDataFrame(crimes[["id", "category"]], crs=4326,
                          geometry=gpd.points_from_xy(crimes["lon"], crimes["lat"])).to_crs(27700)
points = points[points.within(leeds.geometry.iloc[0])].reset_index(drop=True)

# homes: the 2021 census output areas of Leeds, as their centroids
OA = f"{ONS}/Output_Areas_2021_EW_BGC_V2/FeatureServer/0/query"
rows, offset = [], 0
while True:                                        # the service pages at 2,000 records
    page = requests.get(OA, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "OA21CD,BNG_E,BNG_N",
                                    "returnGeometry": "false", "f": "json", "resultOffset": offset}, timeout=120).json()
    rows += [f["attributes"] for f in page["features"]]
    if not page.get("exceededTransferLimit"):
        break
    offset += len(page["features"])
oa = pd.DataFrame(rows).sort_values("OA21CD")
homes = gpd.GeoDataFrame({"oa": oa["OA21CD"]}, crs=27700, geometry=gpd.points_from_xy(oa["BNG_E"], oa["BNG_N"]))

# NaPTAN: the active bus stops of West Yorkshire, with the direction buses travel from each
# NaPTAN area 450 as saved on 9 October 2026 (the register changes weekly); live:
# naptan.api.dft.gov.uk/v1/access-nodes?dataFormat=csv&atcoAreaCodes=450
naptan = pd.read_csv(io.StringIO(requests.get("https://data.spatialworkflow.io/naptan-area-450-2026-10-09.csv", timeout=300).text),
                     low_memory=False)
naptan = naptan[(naptan["Status"] == "active") & naptan["StopType"].isin(["BCT", "BCS"])]
stops = gpd.GeoDataFrame(naptan[["ATCOCode", "CommonName", "Bearing"]].reset_index(drop=True), crs=27700,
                         geometry=gpd.points_from_xy(naptan["Easting"], naptan["Northing"]))

print(f"{len(parcels):,} parcels, {len(points):,} crimes, {len(homes):,} homes, {len(stops):,} bus stops")

293,729 parcels, 9,227 crimes, 2,607 homes, 16,450 bus stops

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)

t0 = time.perf_counter()
overlaps = find_overlaps(parcels, min_area=0.01)
print(f"{len(overlaps):,} overlapping pairs, {overlaps['overlap_m2'].sum():,.0f} m² total, "
      f"{time.perf_counter() - t0:.1f} s")
print(overlaps["overlap_m2"].describe().round(2).to_string())
print(overlaps.head(3).to_string(index=False))

47,911 overlapping pairs, 34,109,227 m² total, 15.3 s
count      47911.00
mean         711.93
std        17848.28
min            0.01
25%            0.03
50%            0.13
75%            1.47
max      1834536.62
       a        b   overlap_m2
58826449 58748224 1.834537e+06
28916412 55095526 1.507681e+06
62036776 28878818 1.232667e+06

Those are real overlaps between registered titles in Leeds. INSPIRE polygons are drawn to index the register, not surveyed as a seamless coverage, so neighbouring titles routinely overlap by slivers where their boundaries were captured separately. Most overlaps are tiny; a few are large, where one title sits over another — a lease-like arrangement, a flying freehold or an error.

left < right does two things at once: it removes self-matches (where left == right) and it keeps each pair in one direction only, so every real overlap is reported once rather than twice.

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`."""
        left, right = self.tree.query(np.asarray(geoms), predicate="dwithin", distance=distance)
        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(stops)

homes["stops_500m"] = hood.counts_within(homes.geometry, 500)
print(homes["stops_500m"].describe().round(1).to_string())
print(f"{(homes['stops_500m'] == 0).sum()} output areas have no bus stop within 500 m")

centre = hood.within(Point(429_932, 433_280), 250, columns=["CommonName"])   # Leeds station
print(f"{len(centre)} stops within 250 m of Leeds station")

count    2607.0
mean       14.4
std         8.9
min         0.0
25%         9.0
50%        14.0
75%        19.0
max       103.0
72 output areas have no bus stop within 500 m
14 stops within 250 m of Leeds station

Building the index once in __init__ is the point of the class. A loop that calls within for each of 2,607 homes pays the build cost once rather than 2,607 times.

predicate="dwithin" is the exact test, not an approximation. An earlier version of this class buffered each home with quad_segs=2 and queried with intersects, on the reasoning that the exact distance test afterwards would clean up. It cannot clean up what the query never returned: an eight-sided buffer reaches only 92% of the radius midway between its vertices, and on these homes and stops it found 34,101 of the 37,623 pairs within 500 m — 9% of the answer silently missing, as the run below shows. The later distance test removes false positives; nothing removes false negatives.

true = len(hood.within_many(homes.geometry, 500))
for segs in (2, 4, 8, 16):                          # a buffer polygon, then intersects: the earlier way
    rings = shapely.buffer(homes.geometry.values, 500, quad_segs=segs)
    left, right = hood.tree.query(rings, predicate="intersects")
    found = (shapely.distance(homes.geometry.values[left], hood.geoms[right]) <= 500).sum()
    print(f"quad_segs={segs:<2}  {found:,} of the {true:,} pairs within 500 m ({found / true:.1%})")

quad_segs=2   34,101 of the 37,623 pairs within 500 m (90.6%)
quad_segs=4   36,768 of the 37,623 pairs within 500 m (97.7%)
quad_segs=8   37,408 of the 37,623 pairs within 500 m (99.4%)
quad_segs=16  37,568 of the 37,623 pairs within 500 m (99.9%)

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
    l_idx, r_idx = tree.query(left.geometry.values, predicate="dwithin", distance=max_distance)

    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][0]][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 bus stop within 2 km where buses head towards the city centre"
COMPASS = {"N": 0, "NE": 45, "E": 90, "SE": 135, "S": 180, "SW": 225, "W": 270, "NW": 315}
CENTRE = Point(429_932, 433_280)                  # Leeds station


def heads_to_centre(home, stop):
    """Buses at `stop` travel within 45° of the direction from `home` to the centre."""
    if stop["Bearing"] not in COMPASS:            # 267 stops have no bearing recorded
        return False
    towards = np.degrees(np.arctan2(CENTRE.x - home.geometry.x, CENTRE.y - home.geometry.y)) % 360
    return abs((COMPASS[stop["Bearing"]] - towards + 180) % 360 - 180) <= 45


inbound = join_with_rule(homes, stops, max_distance=2_000, rule=heads_to_centre,
                         right_cols=["CommonName", "Bearing"])
nearest = gpd.sjoin_nearest(homes, stops[["geometry"]], distance_col="d")
nearest = nearest[~nearest.index.duplicated()]["d"]
print(f"median distance: any stop {nearest.median():.0f} m, "
      f"an inbound stop {inbound['match_distance_m'].median():.0f} m")

2,607 of 2,607 matched
median distance: any stop 142 m, an inbound stop 214 m

Every home found an inbound stop within 2 km, but at a median of 214 m rather than the 142 m to the nearest stop of any direction: the nearest stop is often on the wrong side of the road.

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 every home-stop pair (2,607 × 16,450, about 43 million) to those within 2 km, and lexsort puts the nearest candidate first for each left row, so the loop moves on at 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.

The R-tree behind sjoin

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 indexing all 293,729 Leeds parcels took about a quarter of a second.

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 region you expand to must contain the distance, as a box does, never sit inside it, as a polygon buffer does.

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.

Positional results, stale indexes and predicate direction

• query returns positional indices. Use .iloc, or gdf.index[idx] to get labels.
• The index is built lazily and cached per object. A filter or copy produces 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 predicate reads from the query to the tree. tree.query(point, predicate="within") finds the polygons a point is within; "contains" there asks whether the point contains a polygon, and returns nothing.
• Use predicate="dwithin" for radius queries, not a buffer and intersects: a buffer polygon sits inside the true circle and misses features near the edge — 9% of the pairs in Example 2 with quad_segs=2.
• The vectorised query(geoseries) returns (2, n) positions, indexed frame second.
• nearest(..., return_all=False) keeps one result per input when distances tie.
• max_distance on nearest bounds the search and omits inputs with no neighbour inside it.
• sindex.query on an empty frame returns empty arrays rather than raising, so guard with len(...).
• Predicate cost scales with vertex count. A few enormous geometries dominate the exact phase; consider simplifying.
• shapely.STRtree is 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 — or let predicate="dwithin", distance=... do all three. For one point against the Leeds parcels it was about 40 times faster than testing every parcel, and the gap grows with the size of the layer.

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.