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.

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.