13 min read 2788 words

Calculating 5km proximity buffers around substations in Shapely

Scenario / symptom: you call substation.buffer(5000) on a shapely.geometry.Point in EPSG:4326 expecting a 5 km exclusion circle, and instead buffer.area prints 78539816.34 — a number in square degrees, not square metres. Plotted, the buffer swallows several states; passed into a siting query it qualifies thousands of phantom parcels. No exception is raised. This failure lands in the buffer-generation stage of the screening pipeline, and it is the single-asset root case of the additive over-allocation dissected in the parent workflow, Grid Capacity Buffer Analysis. Because Shapely and the underlying GEOS engine operate purely in planar Cartesian space, the literal 5000 is interpreted in whatever units the coordinates carry — and geographic coordinates carry degrees.

This page isolates the compounding causes, gives a pre-flight guard that surfaces the bug before a single buffer is generated, then builds a projected, memory-safe pipeline that produces audit-ready 5 km zones suitable for interconnection routing and environmental screening. The fix is the same discipline applied across the site: enforce a projected coordinate reference system before any distance call, buffer in metres, and gate the output against an equal-area area check.

Root-cause analysis

The distorted buffer is not one bug — it is up to four causes that each pass silently on their own:

  1. Planar-vs-geodetic mismatch. GEOS treats buffer(5000) as 5000 map units. In EPSG:4326 those units are decimal degrees, so the result is a ~5000-degree disc that wraps the globe rather than a 5 km circle. At 34°N one degree of longitude is ≈ 92 km, so the radius is overstated by roughly five orders of magnitude.
  2. Latitude-dependent distortion. Even a “small” degree-based buffer is an anisotropic ellipse, not a circle: a degree of longitude shrinks with cos(latitude) while a degree of latitude stays roughly constant, so the east–west and north–south reach diverge and the exclusion zone is wrong by the cosine of the latitude.
  3. Static UTM-zone assignment. Hard-coding one projected CRS for a national substation set introduces >1% linear distortion for assets far from the central meridian, quietly invalidating a compliance setback near a zone edge.
  4. Invalid input geometry. Duplicate vertices, self-intersecting rings, or NaN coordinates raise a GEOSException mid-batch — or worse, produce an invalid buffer polygon that downstream union/overlay calls silently drop.

The relationship between radius and screened area is quadratic, , which is why a units error of this magnitude does not merely shift the answer — it detonates it. A correct 5 km radius screens (≈ 78.54 km²); the degree-based buffer screens an area larger than many countries.

Degree-based buffer versus a projected metre-based buffer The broken path takes a Point in EPSG:4326 straight into buffer(5000); GEOS reads 5000 as map units, which are decimal degrees, so the result is a roughly 600,000 km blob that swallows several states and raises no exception. The corrected path reprojects the point to its local UTM zone, buffers 5000 in metres for a true 5 km radius, then reprojects back to EPSG:4326 with a crs_source lineage tag, producing an audit-ready 78.54 km exclusion zone. BROKEN — RADIUS TREATED AS DEGREES CORRECTED — RADIUS IN METRES Point geometry EPSG:4326 · units = ° .buffer(5000) 5000 read as degrees ≈ 600,000 km² blob swallows several states no exception is raised Point geometry EPSG:4326 Transform → UTM always_xy=True .buffer(5000) metres → true 5 km Transform → 4326 store + crs_source tag ≈ 78.54 km² zone audit-ready

Pre-flight validation

Catch the bug before it propagates. This guard inspects the source CRS and the magnitude of a trial buffer’s area, raising a precise error instead of letting a degree-based blob flow downstream. Run it once per input layer at the head of the pipeline.

python
import math
from pyproj import CRS
from shapely.geometry import Point

def assert_buffer_is_metric(
    substation_pt: Point,
    source_epsg: int,
    buffer_meters: float = 5000.0,
) -> dict:
    """Surface the planar/geodetic mismatch before generating real buffers.

    Raises ValueError if the working CRS is geographic (degrees) or if a
    trial buffer's area is wildly off the expected pi*r^2 metric target.
    """
    crs = CRS.from_epsg(source_epsg)
    audit = {"source_epsg": source_epsg, "is_geographic": crs.is_geographic}

    if crs.is_geographic:
        raise ValueError(
            f"EPSG:{source_epsg} is geographic (units=degrees). "
            f"buffer({buffer_meters}) would be interpreted as degrees — "
            "reproject to a projected metre-based CRS (e.g. a local UTM) first."
        )

    expected_area_m2 = math.pi * buffer_meters**2
    trial_area = substation_pt.buffer(buffer_meters).area
    ratio = trial_area / expected_area_m2
    audit["area_ratio"] = round(ratio, 4)
    if not 0.99 <= ratio <= 1.01:
        raise ValueError(
            f"Buffer area {trial_area:.1f} deviates from expected "
            f"{expected_area_m2:.1f} m^2 (ratio {ratio:.3f}) — likely a unit "
            "or projection-distortion error."
        )
    return audit

The naive failure it exists to stop is one line:

python
from shapely.geometry import Point

substation = Point(-118.2437, 34.0522)   # EPSG:4326, Los Angeles
distorted = substation.buffer(5000)       # interpreted as 5000 DEGREES
print(f"{distorted.area:.2f}")            # -> 78539816.34  (square degrees!)

Fix implementation

The corrected pipeline derives a UTM zone per asset from its own longitude and latitude, transforms into that metre-based frame with an explicit pyproj.Transformer (always_xy=True to lock lon/lat order), buffers in metres, repairs any invalid output, then transforms back to EPSG:4326 for storage. It streams features in bounded chunks so a national substation set never materialises in memory at once, and it tags every output with the CRS it was computed in for lineage.

python
import logging
from typing import Any, Iterator

import pyproj
from shapely.geometry import Point, mapping
from shapely.ops import transform
from shapely.validation import make_valid

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
logger = logging.getLogger("substation_buffer_pipeline")


def utm_epsg_for(lon: float, lat: float) -> int:
    """Return the EPSG code of the UTM zone containing (lon, lat)."""
    zone = int((lon + 180) / 6) + 1
    return (32600 if lat >= 0 else 32700) + zone  # 326xx north, 327xx south


def generate_substation_buffers(
    substations: Iterator[dict[str, Any]],
    buffer_meters: float = 5000.0,
    chunk_size: int = 2500,
) -> Iterator[list[dict[str, Any]]]:
    """Yield validated GeoJSON features of 5 km buffers, computed in metres.

    Input features carry EPSG:4326 (lon, lat) coordinates; output geometry is
    returned in EPSG:4326 with a `crs_source` UTM tag for audit reproducibility.
    """
    chunk: list[dict[str, Any]] = []

    for idx, sub in enumerate(substations):
        coords = sub.get("geometry", {}).get("coordinates", [None, None])
        if any(c is None for c in coords):
            logger.warning("record %s skipped: missing coordinates", idx)
            continue

        lon, lat = coords[0], coords[1]
        target_epsg = utm_epsg_for(lon, lat)
        to_utm = pyproj.Transformer.from_crs(4326, target_epsg, always_xy=True)
        to_wgs84 = pyproj.Transformer.from_crs(target_epsg, 4326, always_xy=True)

        try:
            utm_point = transform(to_utm.transform, Point(lon, lat))
            utm_buffer = utm_point.buffer(buffer_meters)        # metres -> true 5 km
            if not utm_buffer.is_valid:
                utm_buffer = make_valid(utm_buffer)
            wgs84_buffer = transform(to_wgs84.transform, utm_buffer)
        except Exception as exc:  # GEOSException, transform failure
            logger.error("record %s failed during buffer/transform: %s", idx, exc)
            continue

        chunk.append({
            "type": "Feature",
            "properties": {
                **sub.get("properties", {}),
                "buffer_m": buffer_meters,
                "crs_source": f"EPSG:{target_epsg}",
            },
            "geometry": mapping(wgs84_buffer),
        })

        if len(chunk) >= chunk_size:
            yield chunk
            chunk = []

    if chunk:
        yield chunk

Explicit parameter choices, justified for grid-GIS use: the per-asset UTM zone keeps linear distortion under ~0.04% near the central meridian (far tighter than the 1% a static zone risks); always_xy=True eliminates the legacy lon/lat axis-swap that silently mirrors geometry in pyproj 6+; make_valid is a deterministic repair rather than a buffer(0) heuristic; and chunk_size=2500 bounds peak heap during the GEOS C-extension calls and any downstream serialisation. The buffer radius itself should not stay a hard-coded scalar for production capacity work — derive it per asset from voltage class and thermal rating as the parent Grid Capacity Buffer Analysis workflow does, and feed it geometry that has already passed spatial data quality validation.

Streaming per-asset buffer pipeline with a bounded-memory chunk buffer EPSG:4326 features stream through five per-record stages: utm_epsg_for picks the local UTM zone, a Transformer projects the point and buffers 5000 metres for a true 5 km radius, make_valid repairs any invalid polygon, and a second Transformer reprojects back to EPSG:4326 with a crs_source tag. Each finished feature is appended to a chunk buffer that lives inside a dashed bounded-memory region capped at chunk_size 2500. When the chunk fills it is yielded to downstream proximity work and reset, so the full national set never materialises in memory at once. Stream features utm_epsg_for() to_utm + buffer make_valid() to_wgs84 EPSG:4326 lon,lat per-asset zone 5000 m → 5 km repair output back to 4326 BOUNDED HEAP chunk_size=2500 Chunk buffer accumulate yield chunk len ≥ 2500 → downstream proximity reset []

Fallback routing & performance tuning

  • Polar and offshore assets (>84° latitude): UTM is undefined beyond the standard zones, so fall back to a global equal-area frame (EPSG:6933) or a regional state-plane CRS, and log the substitution explicitly rather than letting utm_epsg_for return a meaningless zone.
  • Zone-boundary substations: for a point within 5 km of a UTM zone edge the buffer straddles two zones; reproject all affected buffers to one shared frame and shapely.union_all() adjacent zones before export to keep exclusion polygons contiguous.
  • Batch scale (>100k assets): push the per-asset transform into a vectorised geopandas.GeoSeries.to_crs() grouped by UTM zone, or distribute chunks with dask-geopandas, so transformer construction is amortised across the group instead of rebuilt per record.
  • Spatial-index reuse: when buffers feed a proximity query, build an STRtree once over the buffered set rather than per query — see proximity & distance calculations for the indexing pattern.
  • CI/CD memory ceiling: if GEOSException persists under tight runners, lower chunk_size to 500 and confirm the GEOS C-extension is linked against the expected libgeos; pin versions so the deterministic gate runs against the same engine as production.

Downstream validation

Before buffers reach a permitting or routing engine, assert their integrity in a form a CI/CD gate can run. This audit checks validity, emptiness, the metric-area sanity bound, and that every feature carries its crs_source lineage tag.

python
import math
from shapely.geometry import shape

def audit_buffer_features(
    features: list[dict],
    buffer_meters: float = 5000.0,
    tolerance: float = 0.02,
) -> dict:
    """Assert buffer outputs are valid, non-empty, correctly sized and tagged."""
    expected_area_km2 = math.pi * (buffer_meters / 1000.0) ** 2  # ~78.54 km^2
    report = {"checked": 0, "failures": []}

    for feat in features:
        report["checked"] += 1
        geom = shape(feat["geometry"])
        props = feat.get("properties", {})

        if geom.is_empty or not geom.is_valid:
            report["failures"].append((props.get("id"), "invalid_or_empty"))
        if "crs_source" not in props:
            report["failures"].append((props.get("id"), "missing_crs_lineage"))

        # area sanity in an equal-area frame (EPSG:6933), reported in km^2
        from pyproj import Geod
        area_m2 = abs(Geod(ellps="WGS84").geometry_area_perimeter(geom)[0])
        ratio = (area_m2 / 1e6) / expected_area_km2
        if abs(ratio - 1.0) > tolerance:
            report["failures"].append((props.get("id"), f"area_ratio={ratio:.3f}"))

    assert not report["failures"], f"buffer audit failed: {report['failures'][:5]}"
    return report

A 5 km buffer that audits to ≈ 78.54 km² with a valid geometry and a crs_source tag is reproducible: an interconnection study or environmental reviewer can confirm exactly which projection produced each zone, which is what turns a map artefact into evidence.

Related articles

Back to Grid Capacity Buffer Analysis