Vectorized Nearest-Substation Search with a KDTree #

Answering “which substation is closest to each of my candidate sites, and how far?” for tens of thousands of points is the highest-frequency query in an interconnection screen, and the naive loop that calls site.distance(substation) for every pair is the wrong tool for it. This page builds a vectorized nearest-substation lookup on scipy.spatial.cKDTree — the point-to-point workhorse behind the proximity distance calculations workflow — and eliminates the four failure modes that make a first draft return distances that are wrong, slow, or quietly missing the substation that actually matters. A KDTree turns a nearest-neighbour query from a pairwise scan into a logarithmic tree descent, but only when the tree is built over the right coordinates, in the right units, exactly once.

A cKDTree indexes points. It answers point-to-point nearest-neighbour queries in Euclidean space and nothing else — it has no concept of a line, a polygon, or a geographic coordinate system. That single constraint is the source of every failure below: feed it degrees and it measures in degrees; feed it a substation footprint polygon and it silently uses the centroid; ask it for the single nearest node and it hands you a substation with zero available headroom.

Root-cause analysis #

Four compounding causes turn a one-line tree.query() into a corrupted screen, and each maps to a distinct fix in the corrected function below.

1. Tree built on lon/lat degrees. cKDTree computes straight-line Euclidean distance over whatever numbers you hand it. Build the tree from an EPSG:4326 coordinate array and every returned “distance” is in degrees — a unit that mixes a longitudinal axis whose ground length collapses toward the poles with a latitudinal axis that does not. The query still returns a nearest index, so nothing raises; the ranking is simply wrong wherever the north–south and east–west scale factors diverge, and the reported distance is meaningless to a setback test.
2. k=1 misses the usable substation. The single nearest asset is frequently the wrong answer for interconnection: the closest substation may be capacity-saturated, at the wrong voltage class, or already queued out. A k=1 query has no fallback — it returns one node and stops. You need the k nearest so downstream capacity logic can walk outward to the first substation with real headroom.
3. Tree rebuilt on every query. Constructing the KDTree is the expensive step. Rebuilding it inside a per-site loop or per-chunk call throws away the entire advantage of the structure and reintroduces the quadratic cost the index exists to remove.
4. Point index against line/polygon geometry. cKDTree needs an (M, 2) array of coordinates. If your substations are stored as polygon footprints, or you conflate them with transmission LineString geometry, there is no single coordinate to index — the nearest point on a line is not a tree lookup. That case belongs to geopandas.sjoin_nearest, not a cKDTree.

cKDTree vs. sjoin_nearest — pick by geometry #

The decision is entirely about what geometry you are measuring to.

If the target is a substation point and both layers sit in a projected metre frame, cKDTree is the fast path. If you are measuring to a conductor LineString or a substation footprint polygon, use geopandas.sjoin_nearest, which computes true point-to-geometry distance and returns the matched row — at the cost of the geometry-aware distance evaluation the tree avoids.

Pre-flight validation #

Surface the two silent corruptions — a geographic CRS and non-finite coordinates — before the tree is ever built. A NaN or inf in the coordinate array does not raise on construction; it poisons the partition and returns garbage neighbours.

import geopandas as gpd
import numpy as np


def preflight_kdtree_inputs(
    substation_gdf: gpd.GeoDataFrame,
    sites_gdf: gpd.GeoDataFrame,
    target_epsg: int = 32614,
) -> None:
    """Raise on the exact root cause before any KDTree is constructed."""
    for name, gdf in (("substation_gdf", substation_gdf), ("sites_gdf", sites_gdf)):
        if gdf.crs is None:
            raise ValueError(f"{name} has no CRS; KDTree distances would be undefined.")
        if gdf.crs.is_geographic:
            raise ValueError(
                f"{name} is in geographic CRS {gdf.crs.to_epsg()} (degrees). "
                f"Reproject to a projected metre frame (e.g. EPSG:{target_epsg}) first."
            )
        if gdf.crs.to_epsg() != target_epsg:
            raise ValueError(
                f"{name} is EPSG:{gdf.crs.to_epsg()}, expected EPSG:{target_epsg}; "
                "both layers must share one projected frame."
            )
        if not (gdf.geom_type == "Point").all():
            raise TypeError(
                f"{name} holds non-point geometry; cKDTree indexes points only — "
                "use geopandas.sjoin_nearest for line/polygon targets."
            )
        coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
        if not np.isfinite(coords).all():
            raise ValueError(f"{name} contains non-finite coordinates (NaN/inf).")

The is_geographic check is the load-bearing one: it is the difference between a distance in metres and a distance in degrees, and it is the coordinate reference system alignment discipline every distance calculation on this site depends on. EPSG:32614 (UTM zone 14N) is chosen here for a mid-continent portfolio; pick the UTM zone or state-plane code that covers your region.

Fix implementation #

The corrected function builds the tree once over the substation coordinates, then issues a single vectorized query for all site points asking for the k nearest. It returns a tidy frame carrying nearest_substation_id and distance_m — the identifier, not just an index, so the result survives a reindex downstream.

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


def nearest_substations(
    sites_gdf: gpd.GeoDataFrame,
    substation_gdf: gpd.GeoDataFrame,
    id_col: str = "substation_id",
    k: int = 3,
    distance_upper_bound_m: float = 25_000.0,
    target_epsg: int = 32614,
) -> pd.DataFrame:
    """
    Vectorized k-nearest-substation lookup on a single cKDTree.

    Builds the tree once over substation points and queries every site at once.
    Returns one row per (site, neighbour rank) with the substation id and the
    Euclidean distance in metres, so downstream capacity logic can walk outward
    from k=0 to the first substation with real headroom.
    """
    preflight_kdtree_inputs(substation_gdf, sites_gdf, target_epsg)

    # Build the index ONCE over the substation coordinate array (M, 2)
    sub_coords = np.column_stack([substation_gdf.geometry.x, substation_gdf.geometry.y])
    tree = cKDTree(sub_coords)

    # Single vectorized query for ALL sites; k neighbours each
    site_coords = np.column_stack([sites_gdf.geometry.x, sites_gdf.geometry.y])
    distances, positions = tree.query(
        site_coords, k=k, distance_upper_bound=distance_upper_bound_m
    )

    # scipy returns 1-D arrays when k == 1; normalise to 2-D for uniform handling
    if k == 1:
        distances = distances[:, None]
        positions = positions[:, None]

    sub_ids = substation_gdf[id_col].to_numpy()
    site_ids = sites_gdf.index.to_numpy()

    records = []
    for rank in range(k):
        pos = positions[:, rank]
        dist = distances[:, rank]
        # positions == len(tree.data) marks "no neighbour within the bound"
        found = pos < len(sub_ids)
        records.append(pd.DataFrame({
            "site_id": site_ids,
            "neighbour_rank": rank,
            "nearest_substation_id": np.where(found, sub_ids[pos % len(sub_ids)], None),
            "distance_m": np.where(found, dist, np.inf),
        }))

    return pd.concat(records, ignore_index=True).sort_values(
        ["site_id", "neighbour_rank"]
    ).reset_index(drop=True)

The parameter choices are deliberate. k=3 gives capacity logic two fallback substations to reach past a saturated nearest node without a second query. distance_upper_bound_m=25_000 prunes the search so a site with no substation within interconnection reach returns inf rather than a spurious match hundreds of kilometres away — scipy encodes that “not found” case by returning an index equal to the tree size, which the pos < len(sub_ids) mask converts to an explicit infeasible result. The neighbour distance measured here is a straight-line lower bound; reconcile it against real headroom via Grid Capacity Buffer Analysis thresholds before treating any substation as usable.

The Euclidean distance the tree minimises is the ordinary planar norm in the projected frame,

which is only a true ground distance because both coordinate arrays are in metres. The reason to prefer the tree is asymptotic: a pairwise scan of $N$ sites against $M$ substations costs

distance evaluations, while a balanced KDTree answers each nearest query by descending the tree, reducing the batch to

after an $O(M\log M)$ one-time build. For 40,000 sites against 6,000 substations that is roughly a fifty-fold reduction in comparisons — realised only if the tree is built once and reused.

Fallback routing & performance tuning #

• Build once, query in bulk. Construct the cKDTree a single time and pass the entire site coordinate array to one tree.query call. scipy vectorizes the batch internally; a Python-level per-site loop is strictly slower and defeats the structure.
• Raise k when the nearest is saturated. If capacity screening rejects the rank=0 substation, walk to rank=1, rank=2, and so on. Size k to how deep your headroom fallback realistically goes — k=3 to k=5 covers dense grids; going higher wastes query time and memory.
• Tune distance_upper_bound to the reach you screen for. A tight bound prunes far substations and makes “infeasible” explicit; too tight and every site returns inf. Match it to the maximum economic tie-line length, not an arbitrary default.
• Chunk the site array for huge portfolios. For millions of sites, query in slices of a few hundred thousand to cap the peak result-array size, reusing the same tree across slices — never rebuild it per chunk.
• Pass workers=-1 for multi-core queries. tree.query(..., workers=-1) parallelises the batch across all cores; the build stays single-threaded but the query dominates at scale.
• Switch to sjoin_nearest for geometry targets. The moment you must measure to a conductor centreline or a substation footprint edge rather than a representative point, the KDTree is the wrong structure — see the spatial index and proximity quick reference for the full index-selection matrix.

Downstream validation #

Gate the result before it feeds a capacity screen. This assertion catches the three regressions that a KDTree change tends to introduce: non-finite or negative distances, out-of-range substation identifiers, and a silently dropped neighbour rank. It is written to fail a CI/CD build.

import numpy as np
import pandas as pd


def assert_nearest_integrity(
    result: pd.DataFrame,
    substation_gdf: gpd.GeoDataFrame,
    id_col: str = "substation_id",
    expected_k: int = 3,
) -> None:
    """CI/CD gate: fail the build if the nearest-substation table is unsound."""
    finite = result["distance_m"].replace(np.inf, np.nan).dropna()
    assert (finite >= 0).all(), "negative distance present — coordinate or CRS corruption"
    assert np.isfinite(finite).all(), "non-finite distance where a match was reported"

    matched = result.loc[result["nearest_substation_id"].notna(), "nearest_substation_id"]
    valid_ids = set(substation_gdf[id_col])
    assert set(matched).issubset(valid_ids), "matched id absent from substation table"

    ranks_per_site = result.groupby("site_id")["neighbour_rank"].nunique()
    assert (ranks_per_site == expected_k).all(), (
        f"every site must carry {expected_k} neighbour ranks; "
        "a dropped rank means the k-query or concat regressed"
    )

Asserting that every matched identifier exists in the substation table is what turns the %-guarded index arithmetic above into a provable invariant rather than an assumption. Pin scipy, geopandas, and shapely in pyproject.toml so a change to the query return signature or the distance_upper_bound sentinel cannot silently shift results between runs — the same lineage discipline that keeps the broader proximity distance calculations workflow reproducible.

Related #

• Proximity Distance Calculations — the parent workflow this nearest-substation lookup plugs into.
• Automating Interconnection Queue Screening with Async Proximity Scoring — where the k-nearest results drive concurrent capacity scoring.
• Spatial Index & Proximity Quick Reference — the full comparison of cKDTree, R-tree, and sjoin_nearest by geometry and query type.
• Coordinate Reference Systems for Energy Projects — the projected-frame enforcement that makes KDTree distances mean metres.