Proximity Distance Calculations #

Proximity distance calculations are the primary feasibility filter for renewable interconnection projects, and they sit at the analytical core of the Grid Infrastructure & Network Proximity Analysis pipeline. The specific failure mode this page addresses is pairwise O(N×M) proximity scaling: the moment a screening workflow tries to answer “how far is each candidate site from the nearest grid asset?” by looping every one of N candidate generation sites against every one of M transmission features, the run time and memory footprint explode. Fifty thousand sites against a few hundred thousand line segments is two-and-a-half billion distance evaluations — a calculation that finishes in a notebook demo at toy scale and never returns on a continental portfolio. The naive script does not raise an error; it simply hangs, gets killed by the out-of-memory reaper, or quietly returns distances that are wrong by the cosine of the latitude because the geometries were never projected.

This page builds a deterministic proximity-scoring workflow that turns raw asset and candidate geometries into audit-ready feasibility scores. It follows the order the data actually travels: inputs are forced into a projected coordinate frame and topologically validated, the search space is pruned with an R-tree spatial index before any precise geometric operation runs, distances are computed in bounded memory chunks, network-constrained corridors are resolved asynchronously where straight-line distance is meaningless, and every output row carries the capacity and regulatory flags an interconnection queue submission needs. Raw Euclidean metrics are only the starting point — terrain, right-of-way limits, and regulatory setbacks all bend the real interconnection cost away from the straight line.

Why Naive Distance Calculations Fail #

The brute-force approach fails for four compounding reasons, and only one of them reliably raises an exception at the point of error.

First, projected-distance error. Geographic coordinate systems such as EPSG:4326 express position in decimal degrees, and Shapely’s geometry.distance() operates in planar Cartesian space. Call site.distance(line) on unprojected lon/lat and you get a number in degrees that mixes a longitudinal axis whose metric value collapses toward the poles with a latitudinal axis that does not — the result is not off by a rounding error, it is off by a latitude-dependent scale factor. Every distance must be computed in a projected frame whose units are meters, which is exactly the coordinate reference system alignment discipline the rest of the pipeline depends on.

Second, quadratic search-space blow-up. A direct double loop is

distance evaluations. An R-tree spatial index reduces a nearest-feature query to roughly

by pruning every grid feature whose bounding box cannot possibly contain the nearest geometry before a single exact distance is measured. On a 50,000 × 300,000 problem that is the difference between billions of operations and tens of millions.

Third, memory spike. Materializing a dense N×M distance matrix — or calling unary_union on an entire national transmission layer at once — allocates gigabytes that the host never reclaims inside a long-running batch. Bounded chunking with explicit cleanup keeps the resident set flat regardless of portfolio size.

Fourth, async latency on network-constrained legs. When a straight line is meaningless — a candidate separated from the grid by a ridge, a protected wetland, or a missing right-of-way — the real distance comes from a routing service or a cost-surface solver. Issuing those calls synchronously, one site at a time, serializes thousands of independent I/O waits into an unusable wall-clock time.

Prerequisites & Data Requirements #

This workflow assumes the following inputs and environment:

• Candidate sites as a GeoDataFrame of Point (parcel centroids or proposed array centers) with a defined, non-null CRS. Land-cover and parcel preprocessing should already have passed through spatial data quality validation.
• Grid assets as a GeoDataFrame of LineString/MultiLineString conductors and Point substations, with voltage_kv and available_capacity_mw attributes confirmed by Network Attribute Validation. The geometry itself should originate from validated Transmission Line & Substation Mapping so distances reference real conductor corridors, not simplified centerlines.
• A projected target CRS chosen for the region of interest. For a single UTM zone, EPSG:32610 (UTM 10N) or the appropriate state plane code preserves meter-level distance; a multi-state portfolio needs an equal-area or equidistant conic projection — see the projection-choice walkthrough in aligning EPSG:4326 and EPSG:3857 for solar site mapping.
• Library versions: geopandas >= 0.14, shapely >= 2.0 (vectorized predicates and union_all), pyproj >= 3.4. The sindex.query two-array return signature used below requires Shapely 2.x.

CRS choice is not cosmetic. Using the web-mercator EPSG:3857 frame for distance work introduces scale error that grows with latitude — acceptable for tiles, unacceptable for a 500 m regulatory setback test.

Core Implementation #

The happy path is two stages: normalize and validate every input to the projected frame, then score proximity in memory-bounded chunks using the spatial index to prune before measuring.

import geopandas as gpd
import numpy as np
from shapely.validation import make_valid


def normalize_and_validate(
    gdf: gpd.GeoDataFrame, target_epsg: int = 32610
) -> gpd.GeoDataFrame:
    """
    Validate topology and transform a GeoDataFrame to a projected CRS (meters).
    Repairs invalid rings before projection so a single bad geometry cannot
    abort an overnight batch with no traceable cause.
    """
    if gdf.crs is None:
        raise ValueError("Input GeoDataFrame must have a defined CRS (e.g. EPSG:4326).")

    gdf = gdf.copy()
    # Repair self-intersections and invalid rings prior to projection
    gdf.geometry = gdf.geometry.apply(
        lambda geom: make_valid(geom) if not geom.is_valid else geom
    )

    # Explicit transformation into a projected (meter) frame
    gdf = gdf.to_crs(epsg=target_epsg)

    # Drop empties/invalids introduced by repair or transformation
    valid_mask = gdf.geometry.is_valid & ~gdf.geometry.is_empty
    return gdf[valid_mask].reset_index(drop=True)

With both layers in the same projected frame, the scoring loop queries the R-tree index for each chunk of sites, dissolves only the pruned grid subset, and measures exact distances against that small union:

import pandas as pd
from typing import Generator


def chunked_proximity_scores(
    sites_gdf: gpd.GeoDataFrame,
    grid_gdf: gpd.GeoDataFrame,
    chunk_size: int = 5000,
    search_radius_m: float = 25_000.0,
) -> Generator[pd.DataFrame, None, None]:
    """
    Yield nearest-grid distances in memory-bounded chunks, pruning the search
    space with an R-tree before any exact distance is computed.

    search_radius_m bounds the candidate set per site; sites with no grid
    feature inside the radius are reported as +inf rather than forcing a
    full-layer union.
    """
    grid_gdf = grid_gdf.copy()
    grid_gdf["geometry"] = grid_gdf.geometry.buffer(0)  # cheap topology fix
    grid_sindex = grid_gdf.sindex

    for start in range(0, len(sites_gdf), chunk_size):
        chunk = sites_gdf.iloc[start:start + chunk_size].copy()

        # Bounding-box pre-filter: query each site's search envelope (Shapely 2.x)
        envelopes = chunk.geometry.buffer(search_radius_m)
        site_pos, grid_pos = grid_sindex.query(envelopes, predicate="intersects")

        distances = np.full(len(chunk), np.inf, dtype="float64")
        # Group candidate grid features per site position, then measure exactly
        for s in np.unique(site_pos):
            site_geom = chunk.geometry.iloc[s]
            candidate_idx = grid_pos[site_pos == s]
            candidates = grid_gdf.geometry.iloc[candidate_idx]
            distances[s] = candidates.distance(site_geom).min()

        yield pd.DataFrame({
            "site_id": chunk.index.to_numpy(),
            "nearest_grid_distance_m": distances,
        })

        del chunk, envelopes, site_pos, grid_pos  # keep the resident set flat

The search_radius_m envelope is what keeps the pruned subset small: it bounds the bounding-box query so dense corridors do not degenerate back toward the pairwise case, and it makes “no grid within reach” an explicit inf result rather than an exception.

Error Handling & Edge Cases #

Each of the failure modes named above has a concrete guard.

Unprojected or mismatched CRS. The single most common silent corruption. Assert that both layers share the projected target frame before scoring — never compute distance across a CRS boundary:

def assert_projected_meters(*gdfs: gpd.GeoDataFrame, expected_epsg: int = 32610) -> None:
    for g in gdfs:
        if g.crs is None:
            raise ValueError("Geometry has no CRS; distances would be undefined.")
        if g.crs.to_epsg() != expected_epsg:
            raise ValueError(
                f"CRS {g.crs.to_epsg()} != target EPSG:{expected_epsg}; "
                "distances must be computed in a single projected (meter) frame."
            )
        if g.crs.is_geographic:
            raise ValueError(
                "Geographic CRS detected — distance() would return degrees, not meters."
            )

Sites with no grid feature in range. Returning inf (as the scorer does) is correct, but downstream code must treat it as “infeasible,” not coerce it to a real distance. Tag and partition these rather than dropping them silently — an interconnection screen needs to report why a site failed.

Obstructed straight-line paths. Where a candidate is separated from the grid by terrain or an exclusion zone, the Euclidean nearest distance understates the true interconnection length. Flag any site whose nearest asset lies across a known barrier layer for re-routing in the network-constrained stage below, and never let a straight-line distance silently stand in for a routed one in the feasibility score.

Performance & Scalability — Network-Constrained Routing #

For obstructed legs, the real distance comes from a routing service or a Dijkstra solve over a rasterized impedance surface. These are I/O- and compute-bound and must run concurrently. The pattern below dispatches routing requests asynchronously while validating each spatial input before the call, and preserves order so results align with the input sites:

import asyncio
import aiohttp
from shapely.geometry import Point
from typing import List, Tuple


async def resolve_network_distances(
    site_coords: List[Tuple[float, float]],
    routing_endpoint: str,
    session: aiohttp.ClientSession,
    max_concurrency: int = 32,
) -> List[float]:
    """
    Concurrently fetch network-constrained distances for obstructed sites.
    Validates each geometry before dispatch and bounds concurrency so a large
    portfolio does not exhaust the routing service or local sockets.
    """
    sem = asyncio.Semaphore(max_concurrency)

    async def _fetch(site: Point) -> float:
        if not site.is_valid:
            raise ValueError(f"Invalid site geometry: {site.wkt}")
        payload = {"origin": [site.x, site.y], "mode": "grid_tie"}
        async with sem:
            async with session.post(routing_endpoint, json=payload) as resp:
                resp.raise_for_status()
                data = await resp.json()
                return float(data.get("distance_m", float("inf")))

    tasks = [_fetch(Point(x, y)) for x, y in site_coords]
    results = await asyncio.gather(*tasks, return_exceptions=True)

    # Failed legs degrade to +inf (infeasible) while preserving input order
    return [r if isinstance(r, float) else float("inf") for r in results]

Additional scaling levers for continental runs:

• Build the index once. Construct grid_gdf.sindex a single time and reuse it across every chunk; rebuilding per chunk reintroduces the cost the index was meant to remove.
• Tune chunk_size to the host. Larger chunks amortize per-call overhead but raise the peak resident set; size it against the worker’s memory budget, not a fixed default.
• Bound concurrency, not just parallelism. The Semaphore ceiling protects the routing endpoint and local socket pool — unbounded gather over 50,000 sites is its own outage.
• Spatially partition the portfolio. Process geographically contiguous tiles so each chunk’s pruned grid subset stays small and cache-local.

Validation & Audit Trail #

A distance is only a feasibility input once it is reconciled against capacity and regulatory constraints. The final stage cross-references each distance against Grid Capacity Buffer Analysis thresholds and applies the minimum environmental setback, then emits a bounded score plus the flags an interconnection or permitting reviewer needs:

def apply_compliance_filters(
    proximity_df: pd.DataFrame,
    capacity_threshold_km: float = 15.0,
    regulatory_setback_m: float = 500.0,
) -> pd.DataFrame:
    """
    Reconcile raw proximity against capacity reach and regulatory setback,
    returning a 0-100 feasibility score with explicit, auditable flags.
    """
    df = proximity_df.copy()
    reach_m = capacity_threshold_km * 1000.0

    # Within interconnection reach of a viable asset?
    df["capacity_viable"] = df["nearest_grid_distance_m"] <= reach_m
    # Beyond the minimum regulatory/environmental setback?
    df["regulatory_compliant"] = df["nearest_grid_distance_m"] >= regulatory_setback_m

    # Distance-efficiency score, zeroed when either constraint is violated
    df["feasibility_score"] = np.where(
        df["capacity_viable"] & df["regulatory_compliant"],
        100.0 * (1.0 - (df["nearest_grid_distance_m"] / reach_m)),
        0.0,
    ).clip(0.0, 100.0)

    # Lineage so a reviewer can reproduce the verdict
    df["capacity_threshold_km"] = capacity_threshold_km
    df["regulatory_setback_m"] = regulatory_setback_m
    df["audit_timestamp"] = pd.Timestamp.utcnow().isoformat()
    return df

The capacity_threshold_km, regulatory_setback_m, and audit_timestamp columns are not decorative — they are the provenance that lets a screening result be independently re-run and arrive at the same verdict. A feasibility score without the thresholds and timestamp that produced it is a number a reviewer has no basis to trust.

At production scale, treat the whole sequence as a deterministic, auditable pipeline rather than an ad-hoc script: enforce schema validation on incoming GeoJSON/Parquet payloads, emit structured logs for every CRS assertion and inf-distance partition, and containerize the async workers to isolate routing-I/O bottlenecks. That discipline is what lets energy developers and GIS engineers scale interconnection feasibility studies across multi-state portfolios while staying inside regional grid codes and environmental permitting standards.

Frequently asked questions #

Is geodesic distance worth the cost over projected distance? #

Rarely, inside a study area. A well-chosen projected frame agrees with the geodesic answer to within a few parts per ten thousand over the tens of kilometres an interconnection study spans, and the projected calculation is vectorised while the geodesic one is not. Geodesic distance earns its cost when the pairs span continental distances or cross zone boundaries — a portfolio-level “nearest substation anywhere in the country” query, not a county screen.

How should obstructed routes be handled in a first-pass screen? #

By flagging them, not by routing them. Compute straight-line distance for everything, then mark the legs whose straight line crosses a river, a protected corridor or a rail line, and carry a provisional circuity multiplier for those. Routing is an order of magnitude more expensive per leg and is the right tool for the shortlist, not the population — and a flagged leg is far more useful to a developer than a silently optimistic straight line.

What distance should a screen use when the point of interconnection is unknown? #

The distance to the nearest point on the nearest suitable circuit, with the suitability filter stated. That is not the same as the distance to the nearest substation, and the difference decides whether the project is a tap or a node interconnection — a distinction with a large cost consequence. Where neither is known, report both and let the ranking carry both columns rather than collapsing them into one number.

Why does the screen return a nearest substation that is obviously wrong? #

Almost always because the query ran against a reference set that was not filtered to serviceable assets: a decommissioned yard, a distribution-class node below the project’s voltage, or a duplicate record from a second dataset. The geometry is doing exactly what it was asked. Filter the reference set by status and voltage class before building the index, and the same query returns a defensible answer.

How should distances be reported — to the asset or to its boundary? #

To the boundary of a polygon asset and to the point for a node, stated explicitly either way. A substation mapped as a yard polygon is several hundred metres across, so “distance to the substation” differs by that much depending on which convention is used — enough to reorder a shortlist of similar candidates. Fix the convention once, name it in the output column, and apply it everywhere.

Related #

• Grid Infrastructure & Network Proximity Analysis — the parent pipeline this proximity-scoring stage belongs to.
• Grid Capacity Buffer Analysis — the capacity-reach thresholds the feasibility filter reconciles against.
• Transmission Line & Substation Mapping — the validated asset geometry these distances are measured to.
• Network Attribute Validation — schema enforcement for the voltage and capacity attributes the screen keys off.
• Coordinate Reference Systems for Energy Projects — the projected-frame selection every distance calculation depends on.
• Calculating 5 km Proximity Buffers Around Substations in Shapely — the single-asset walkthrough of the projected-distance failure mode.