Spatial Index & Proximity Quick Reference #

Every proximity question in Grid Infrastructure & Network Proximity Analysis — nearest substation, sites within a clearance buffer, the closest energized corridor to fifty thousand candidate parcels — reduces to one decision made early and often: which spatial index, and which proximity method, for this query shape? Get it wrong and the run is either quadratically slow, silently wrong (distances measured in degrees), or memory-bound. This page is the lookup table for that decision. It sits alongside the deep walkthrough in proximity distance calculations and the projection and CRS quick reference, and it assumes the one precondition every method here shares: geometries are already in a projected, metric coordinate reference system such as EPSG:32610 (UTM 10N) or EPSG:5070 (CONUS Albers) before a single distance is measured.

Use these tables as anchors. Cross-links from the rest of the site point here when they need to justify an index choice without re-deriving it in place.

Spatial index types at a glance #

The index is the data structure that prunes the search space before any exact geometry math runs. Pick it by the query you actually issue and the geometry type you hold.

The R-tree sindex is the default for corridor and substation work because it queries true geometry envelopes — a LineString conductor, not a centroid approximation. A cKDTree is faster and lighter but only understands points, so it answers “nearest substation location” cleanly and “nearest line” only via midpoint or vertex proxies. STRtree is the right call when the reference layer is static (a fixed transmission network queried by many candidate batches), since it is immutable and cheap to reuse. H3 and geohash trade exactness for near-constant-time bucketing — ideal as a first-pass tile key across a national portfolio, not as a final distance. PostGIS GiST moves the whole problem server-side when the grid dataset outgrows a worker’s RAM.

Choosing a proximity method #

The method is the operation you call on top of the index. Each maps to a distinct query shape.

For point-to-point screening — “which substations are near this site?” — cKDTree.query with k>1 returns a ranked shortlist in one vectorized call. For the authoritative geometry-to-geometry answer that feeds an interconnection screen, sjoin_nearest on the sindex is correct because it measures to the conductor, not a proxy. When the question is membership rather than ranking — every asset inside a right-of-way or environmental setback — buffer then sjoin on intersects is the idiom, and it is exactly how grid capacity buffer analysis tests clearance. Network routing is the fallback reserved for legs where a ridge, wetland, or missing right-of-way makes the Euclidean distance a lie.

Complexity cheat-sheet #

The entire reason an index exists is to move the dominant term from a product to a logarithm. For $N$ candidate sites screened against $M$ grid features:

The pairwise cost

is what kills desktop workflows: it does not raise an error, it simply never returns, or is killed by the out-of-memory reaper. An index replaces it with

by discarding every feature whose bounding box cannot contain the nearest geometry before one exact distance is computed. On a 50,000 × 300,000 problem that is roughly the difference between $1.5\times 10^{10}$ and $9\times 10^5$ operations — four to five orders of magnitude, which is the gap between an overnight job and a sub-second query.

Decision matrix #

Read left to right: the query shape you hold determines the index, which determines the method to call.

cKDTree vs sindex.query in practice #

The two workhorse indexes answer the same question — nearest grid asset — with different trade-offs. A cKDTree over substation coordinates is the fastest possible point-to-point screen but ignores line geometry; the GeoPandas sindex with sjoin_nearest measures true distance to conductors at slightly higher cost. The snippet below runs both against the same inputs so the difference is concrete. Both require the layers to already share a projected frame such as EPSG:32610.

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

TARGET_EPSG = 32610  # UTM 10N — metric, distances in metres


def nearest_substation_kdtree(
    sites_gdf: gpd.GeoDataFrame, substation_gdf: gpd.GeoDataFrame
) -> gpd.GeoDataFrame:
    """Fast point-to-point: nearest substation LOCATION via a KD-tree."""
    assert sites_gdf.crs.to_epsg() == TARGET_EPSG, "Sites must be projected metric"
    assert substation_gdf.crs.to_epsg() == TARGET_EPSG, "Substations must be projected"

    sub_xy = np.column_stack((substation_gdf.geometry.x, substation_gdf.geometry.y))
    site_xy = np.column_stack((sites_gdf.geometry.x, sites_gdf.geometry.y))

    tree = cKDTree(sub_xy)                      # build: O(M log M)
    dist_m, idx = tree.query(site_xy, k=1)      # query: O(N log M)

    out = sites_gdf.copy()
    out["nearest_sub_id"] = substation_gdf.iloc[idx]["substation_id"].to_numpy()
    out["kdtree_distance_m"] = dist_m
    return out


def nearest_substation_sindex(
    sites_gdf: gpd.GeoDataFrame, substation_gdf: gpd.GeoDataFrame
) -> gpd.GeoDataFrame:
    """Authoritative geometry-to-geometry distance via the R-tree sindex."""
    # sjoin_nearest builds and queries substation_gdf.sindex internally
    joined = gpd.sjoin_nearest(
        sites_gdf, substation_gdf[["substation_id", "geometry"]],
        distance_col="sindex_distance_m", how="left",
    )
    # Collapse ties (a site equidistant to two assets) to the first match
    return joined[~joined.index.duplicated(keep="first")]


# Contrast on the same inputs
kd = nearest_substation_kdtree(sites_gdf, substation_gdf)
sj = nearest_substation_sindex(sites_gdf, substation_gdf)
delta = (kd["kdtree_distance_m"] - sj["sindex_distance_m"]).abs()
print(f"max |Δ| between methods: {delta.max():.2f} m")  # ~0 for point layers

For a point substation layer the two agree to floating-point noise, and the KD-tree wins on speed and memory. The moment the reference layer becomes lines or polygons — real transmission corridors — the KD-tree can only see midpoints or vertices, and sjoin_nearest on the sindex becomes the correct choice because it measures perpendicular distance to the conductor itself.

Guidance notes #

• Project first, always. Every method here assumes a metric CRS. Running any of them on EPSG:4326 returns degrees, not metres — the canonical silent bug. Enforce the reprojection covered in the projection and CRS quick reference before indexing.
• Points → KD-tree, geometries → R-tree. Reach for cKDTree when the reference layer is genuinely point-like (substations, met masts) and you want ranked k-nearest. Reach for sindex / sjoin_nearest when distance-to-line or distance-to-polygon must be exact.
• Build the index once. Construct substation_gdf.sindex or the cKDTree a single time and reuse it across every candidate chunk; rebuilding per chunk reintroduces the very cost the index removes.
• Bound the search with max_distance. Passing a max_distance to sjoin_nearest (or a distance_upper_bound to cKDTree.query) caps work per query and makes “no asset within reach” an explicit null instead of a spurious far match.
• Use H3 for a first pass, not the verdict. Cell bucketing is a near-constant-time way to shard a continental portfolio into tiles; refine within each tile with an exact R-tree query rather than trusting the cell distance.
• Push to PostGIS GiST when RAM runs out. When the grid layer no longer fits a worker, the <-> KNN operator over a GiST index keeps the join out-of-core and server-side.

Worked example: sizing an index for a national screen #

A concrete workload makes the trade-offs legible. Screening 42,000 candidate parcels against 8,600 substations, 61,000 transmission ways and 4,200 constraint polygons involves three different query shapes, and each wants a different structure.

The substation query is point-to-point nearest, so a cKDTree built on projected coordinates is the right structure: 0.4 seconds to build, roughly 3 microseconds per query, and an exact answer as long as the coordinates are metric. Because the tree is built once and queried 42,000 times, the build cost is irrelevant — it is amortised on the first few hundred queries.

The transmission query is point-to-line nearest, which a KD-tree cannot answer directly: a tree over line vertices returns the nearest vertex, not the nearest point on the line, and the two differ by up to half a segment length. Here the STRtree behind gdf.sindex is correct — query for candidate geometries by bounding box, then call shapely.distance on the handful that survive. The exact predicate runs on tens of candidates instead of tens of thousands of lines.

The constraint query is point-in-polygon over a modest polygon set with expensive geometries, which is the case prepared geometry was built for. Preparing each constraint polygon once and testing candidates against the prepared version turns a per-test edge walk into an indexed lookup, and the preparation pays for itself after roughly twenty tests against the same polygon.

The fourth structure, an H3 cell hash, answers none of these correctly — it answers a different question very fast. Hashing every parcel and every substation to resolution 8 cells and joining on the cell identifier finds candidates in microseconds, but the answer is only exact to the cell size, which at resolution 8 is about 0.74 square kilometres. That is fine for a portfolio rollup and useless for a setback, so the workable pattern is to use the hash to find candidates and geometry to measure them.

Frequently asked questions #

Why is sindex.query returning features that do not intersect? #

Because it is a bounding-box query by design: it returns candidates whose envelopes overlap, and the exact predicate is the caller’s job. That two-step shape is the whole point — the cheap test prunes the population, the expensive test decides. A workflow that treats the candidate list as the answer over-selects by whatever the difference between the envelopes and the geometries happens to be, which for long diagonal lines is very large.

Does building an index help for a single query? #

No. One query against an unindexed frame is a linear scan; one query against a freshly built index is a linear-time build plus a fast lookup, which is strictly slower. Indexes pay off when the reference set is reused, which is the normal case in screening and the abnormal case in an interactive notebook — where the index is often rebuilt implicitly on every call because the frame was copied in between.

Should the index be rebuilt after a filter? #

Yes, if the filter removed a meaningful share of the reference set, and GeoPandas will do it lazily on first access to sindex after a copy. The subtle failure is the opposite: holding a reference to an index built over the unfiltered frame and querying it with positional indices that now refer to different rows. Always query the index attached to the frame you are indexing into.

How does H3 resolution map to grid work? #

Resolution 6 averages about 36 square kilometres per cell and suits balancing-area rollups; resolution 7 averages 5.2 and suits county-scale capacity aggregation; resolution 8 averages 0.74, close to a utility-scale solar block; resolution 9 averages 0.105, roughly a substation yard. Pick the resolution from the thing being counted, and remember that a cell join is exact only to the cell size.

Is a spatial index useful for temporal filtering too? #

Not directly — but the same principle applies, and the two compose. Filter on time first when the temporal predicate is selective, because dropping rows before a spatial query shrinks both the index build and the candidate set. In a partitioned store the temporal filter is usually a partition prune, which costs nothing at all, and the spatial index then runs over a fraction of the data.

Why does sjoin_nearest return more rows than the left frame? #

Because ties are returned in full by default: when two reference geometries are exactly equidistant, both survive the join, and a downstream aggregation then double-counts that row. Exact ties are common in gridded or snapped data, where several candidates sit at identical rounded distances. Resolve them deterministically — lowest identifier, highest voltage, whatever the domain justifies — rather than letting row order decide.

What is the cheapest way to speed up a slow spatial join? #

Reduce the candidate population before the join rather than optimising the join itself. Filtering the reference set to serviceable assets, projecting to a metric frame once instead of per call, and dropping columns that are not needed downstream routinely produce a larger speed-up than any index change, because they shrink both sides of the operation. Reach for the index next, and for a distributed scheduler last.

Does a spatial index help with contains as well as intersects? #

Yes — every binary predicate benefits from the same candidate-then-refine pattern, because the bounding-box test is a necessary condition for all of them. The refinement step differs, and contains is the more expensive refinement, which makes the pruning more valuable rather than less. Prepared geometry compounds the gain when one side is reused across many tests.

How large can a reference set get before an index stops helping? #

The index keeps helping; what stops scaling is holding the whole reference set in one process. A tree over a few million points is unremarkable, and the query cost grows only logarithmically, so the practical ceiling is memory rather than algorithmic. Past that point the answer is to partition the reference set spatially — by state, by balancing area, by H3 cell — and index each partition, rather than to abandon indexing for a distributed scan.

Should distances be cached between runs? #

Cache the pairings, not the distances. Which substation is nearest to a given site changes rarely, while the distance to it may be recomputed cheaply once the pairing is known — and a cached distance becomes wrong silently when either geometry is edited. Storing the nearest-asset identifier with the inputs that produced it gives the speed-up without the staleness class of bug.

Related #

• Grid Infrastructure & Network Proximity Analysis — the parent pipeline whose proximity stage these indexes power.
• Proximity Distance Calculations — the full chunked, indexed scoring workflow this reference distils.
• Grid Capacity Buffer Analysis — where the buffer-plus-sjoin method decides clearance and setback membership.
• Projection and CRS Quick Reference — the metric-CRS choice every method here depends on.
• Coordinate Reference Systems for Energy Projects — why an unprojected index returns degrees, not metres.