How to Measure Distance to the Grid for Renewable Siting

On this page

Problem statement

A field with excellent sunshine 15 km from the nearest substation is often a worse site than a mediocre one next to it. The connection line has to be built, permitted and paid for, and the substation has to have spare capacity at a voltage the plant can use. Distance to the grid is therefore one of the heaviest factors in any renewable suitability map — and one of the easiest to get wrong, by measuring to the nearest low-voltage line a plant cannot connect to, or to a line a plant cannot connect to without a substation.

OpenStreetMap maps power lines and substations with their voltages in much of the world. With them, distance to the grid becomes a raster on the same grid as the exclusion mask.

Quick answer

Fetch lines and substations, keep the voltages that matter, and compute a distance raster:

import numpy as np, pandas as pd, geopandas as gpd, osmnx as ox, rasterio
from rasterio import features
from scipy import ndimage

ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
WEST, SOUTH, EAST, NORTH = -4.30, 39.70, -3.90, 40.00
PAD = 0.15                                   # fetch well beyond the region

with rasterio.open("siting_mask.tif") as src:        # from the exclusion-mask guide
    mask, transform, crs = src.read(1), src.transform, src.crs
    RES = src.res[0]
inside, buildable = mask != 255, mask == 1

power = ox.features_from_bbox((WEST - PAD, SOUTH - PAD, EAST + PAD, NORTH + PAD),
                              {"power": ["line", "substation"]}).to_crs(crs)

def max_kv(v):
    try:
        return max(float(x) for x in str(v).split(";")) / 1000
    except ValueError:
        return np.nan

power["kv"] = power["voltage"].map(max_kv)
lines = power[power["power"] == "line"]
subs = power[power["power"] == "substation"].copy()
subs["geometry"] = subs.geometry.centroid
print(lines.kv.value_counts(dropna=False).sort_index().to_string())
print(f"{len(subs)} substations, {subs.kv.notna().sum()} with a voltage")

MARGIN = int(12_000 / RES)                   # cells: 12 km, inside the padded fetch

def distance_km(geoms):
    """Distance on the mask's grid widened by MARGIN, so features beyond the edge still count."""
    wide = (mask.shape[0] + 2 * MARGIN, mask.shape[1] + 2 * MARGIN)
    burnt = features.rasterize(((g, 1) for g in geoms), out_shape=wide,
                               transform=transform * transform.translation(-MARGIN, -MARGIN),
                               fill=0, all_touched=True)
    d = ndimage.distance_transform_edt(burnt == 0) * RES / 1000
    return d[MARGIN:-MARGIN, MARGIN:-MARGIN]

d_line = distance_km(lines[lines.kv >= 132].geometry)
d_sub = distance_km(subs[subs.kv >= 45].geometry)
print(f"\nbuildable land: median {np.median(d_line[buildable]):.1f} km to a ≥132 kV line, "
      f"{np.median(d_sub[buildable]):.1f} km to a ≥45 kV substation")

kv
20.0       2
45.0      48
132.0     46
220.0     19
400.0      2
1320.0     1
NaN       15
81 substations, 32 with a voltage

buildable land: median 5.5 km to a ≥132 kV line, 5.9 km to a ≥45 kV substation

Distances are measured to lines and substations outside the region too: the fetch is padded by 0.15° (13–17 km), and the distance is computed on the mask's grid widened by 12 km and then cropped back, so the edge cells are not penalised for a boundary drawn on a map. Burning onto the mask's own grid would silently drop every feature beyond its edge.

What the plant connects to, voltages, padded fetches

1. Decide what the plant connects to

Small plants (a few MW) connect to distribution lines at 15–45 kV; utility-scale plants need 132 kV or more, usually at a substation. The distance that matters is to the right voltage, not to the nearest wire.

2. Parse the voltage tag

OSM stores voltage in volts, and multi-circuit lines list several values separated by semicolons: 132000;20000. Take the maximum for the line's capability; treat unparseable or missing values as unknown, not zero.

3. Pad the fetch area

The nearest substation to a cell near the region's edge may be outside the region. Fetch at least the largest distance you intend to score.

4. Burn and distance-transform on the mask's grid, widened

Rasterise the chosen features onto the exclusion mask's transform, extended by a margin on every side, compute the Euclidean distance, and crop the margin off. The result is aligned with the mask by construction, and features outside the region still count — rasterising onto the mask's exact extent would clip them away, however far the fetch was padded.

5. Summarise the candidate land by distance band

A histogram of buildable area by distance to the grid is the first figure a developer asks for, and it shows immediately whether the good land is also connectable land.

6. Treat straight-line distance as a lower bound

Connection lines follow roads, avoid settlements and cannot cross protected areas freely. Straight-line distance is a useful score, not a route.

The nearest wire is rarely the one a utility-scale plant can use.

Buildable land by distance to substations

The data used on this page

The grid and the buildable land come from How to build a siting exclusion mask in Python: the region 4.30–3.90° W, 39.70–40.00° N south-west of Madrid, on a 25 m grid in ETRS89 / UTM 30N, with 500 m settlement, 100 m main-road and 50 m water setbacks, Natura 2000 sites, slopes over 10% and patches under 10 ha removed. If siting_mask.tif is not already in the working directory, the setup rebuilds it with exactly those rules. Power lines and substations are fetched from OpenStreetMap by the page's own code; the setup unpacks every OpenStreetMap answer as the Overpass mirror gave it on 9 October 2026.

import io
import os
import zipfile

import geopandas as gpd
import numpy as np
import osmnx as ox
import rasterio
import rasterio.merge
import requests
from rasterio import features
from rasterio.transform import from_origin
from rasterio.warp import reproject, Resampling
from scipy import ndimage
from shapely.geometry import box

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-distance-to-grid-python-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.requests_timeout = 600
ox.settings.max_query_area_size = 50e9          # one Overpass query for the padded box, not several

if not os.path.exists("siting_mask.tif"):
    W, S, E, N, CRS, RES = -4.30, 39.70, -3.90, 40.00, "EPSG:25830", 25.0
    region = gpd.GeoSeries([box(W, S, E, N)], crs=4326).to_crs(CRS)
    x0, y0, x1, y1 = (np.floor(b / RES) * RES for b in region.total_bounds)
    shape = (int((y1 - y0) / RES), int((x1 - x0) / RES))
    tr = from_origin(x0, y1, RES, RES)
    inside = features.rasterize(region, out_shape=shape, transform=tr) == 1
    burn = lambda gs, touched=True: features.rasterize(((g, 1) for g in gs), out_shape=shape, transform=tr,
                                                        fill=0, all_touched=touched).astype(bool)
    near = lambda m, d: ndimage.distance_transform_edt(~m) * RES < d
    osm = lambda tags: ox.features_from_bbox((W - 0.03, S - 0.03, E + 0.03, N + 0.03), tags).to_crs(CRS).geometry
    excluded = (near(burn(osm({"landuse": ["residential", "industrial", "commercial", "retail"]})), 500)
                | near(burn(osm({"highway": ["motorway", "trunk", "primary"]})), 100)
                | near(burn(osm({"natural": "water"})), 50))
    srcs = [rasterio.open(f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N39_00_{t}_00_DEM/"
                          f"Copernicus_DSM_COG_10_N39_00_{t}_00_DEM.tif") for t in ("W005", "W004")]
    mosaic, mt = rasterio.merge.merge(srcs, bounds=(W - 0.01, S - 0.01, E + 0.01, N + 0.01))
    dem = np.full(shape, np.nan, dtype="float32")
    reproject(mosaic[0], dem, src_transform=mt, src_crs=srcs[0].crs, dst_transform=tr, dst_crs=CRS,
              resampling=Resampling.bilinear)
    dz_dy, dz_dx = np.gradient(dem, RES)
    natura = []
    for layer in (0, 1):
        gj = requests.get("https://bio.discomap.eea.europa.eu/arcgis/rest/services/ProtectedSites/"
                          f"Natura2000Sites/MapServer/{layer}/query", params={
            "geometry": f"{W},{S},{E},{N}", "geometryType": "esriGeometryEnvelope", "inSR": 4326,
            "outSR": 4326, "spatialRel": "esriSpatialRelIntersects", "outFields": "SITECODE",
            "f": "geojson"}, timeout=120).json()
        natura += list(gpd.GeoDataFrame.from_features(gj["features"], crs=4326).to_crs(CRS).geometry)
    candidate = ~(excluded | burn(natura, touched=False) | (np.hypot(dz_dx, dz_dy) * 100 > 10) | ~inside)
    labels, _ = ndimage.label(candidate)
    keep = np.bincount(labels.ravel()) * RES ** 2 / 1e4 >= 10
    keep[0] = False
    with rasterio.open("siting_mask.tif", "w", driver="GTiff", dtype="uint8", count=1, width=shape[1],
                       height=shape[0], crs=CRS, transform=tr, nodata=255, compress="deflate") as dst:
        dst.write(np.where(inside, keep[labels], 255).astype("uint8"), 1)

with rasterio.open("siting_mask.tif") as src:
    m = src.read(1)
print(f"siting_mask.tif: {m.shape[1]} x {m.shape[0]}, buildable {(m == 1).sum() * 625 / 1e6:,.0f} km²")

siting_mask.tif: 1386 x 1349, buildable 593 km²

Example 1 — buildable land by distance band

bands = [0, 1, 2, 5, 10, 50]
km2 = lambda m: m.sum() * RES**2 / 1e6
print(f"{'distance':<10}{'to ≥132 kV line':>17}{'to ≥45 kV substation':>22}")
for lo, hi in zip(bands[:-1], bands[1:]):
    a = km2(buildable & (d_line >= lo) & (d_line < hi))
    b = km2(buildable & (d_sub >= lo) & (d_sub < hi))
    print(f"{lo:>2}–{hi:<2} km{a:>14.0f} km²{b:>18.0f} km²")

distance    to ≥132 kV line  to ≥45 kV substation
 0–1  km           106 km²                12 km²
 1–2  km            67 km²                28 km²
 2–5  km           108 km²               181 km²
 5–10 km           107 km²               304 km²
10–50 km           205 km²                69 km²

The good land and the grid are in different places. Only 173 of the 593 buildable km² lie within 2 km of a 132 kV line, and 40 km² within 2 km of a substation; 205 km² are more than 10 km from any line at 132 kV or above. The fields that survive the setbacks are the ones away from towns, and the substations are where the towns are. A plant can tap into a line only with a new substation, so the second column is usually the better cost proxy.

Example 2 — what the wrong voltage filter does

for label, sel in (("any line", lines.kv.notna() | lines.kv.isna()),
                   ("≥45 kV", lines.kv >= 45),
                   ("≥132 kV", lines.kv >= 132),
                   ("≥220 kV", lines.kv >= 220)):
    d = distance_km(lines[sel].geometry)
    print(f"{label:<9} {sel.sum():>4} lines, median distance from buildable land "
          f"{np.median(d[buildable]):>5.1f} km, share within 2 km "
          f"{(d[buildable] < 2).mean():>4.0%}")

any line   133 lines, median distance from buildable land   2.0 km, share within 2 km  50%
≥45 kV     116 lines, median distance from buildable land   2.0 km, share within 2 km  50%
≥132 kV     68 lines, median distance from buildable land   5.5 km, share within 2 km  29%
≥220 kV     22 lines, median distance from buildable land  10.8 km, share within 2 km   5%

Measuring to any line puts half of the buildable land within 2 km of the grid; measuring to 220 kV, 5% of it. The filter is the analysis: it encodes which plant size the map is for.

Example 3 — write the distance layer for the suitability model

profile = {"driver": "GTiff", "dtype": "float32", "count": 1, "width": mask.shape[1],
           "height": mask.shape[0], "crs": crs, "transform": transform,
           "nodata": -1, "compress": "deflate"}
with rasterio.open("dist_substation_km.tif", "w", **profile) as dst:
    dst.write(np.where(inside, d_sub, -1).astype("float32"), 1)
unknown = lines.kv.isna().sum()
print(f"written; {unknown} of {len(lines)} lines had no usable voltage "
      f"and were left out of every filter")

written; 15 of 133 lines had no usable voltage and were left out of every filter

Lines with no voltage tag cannot be classified, and they are silently dropped by any >= filter. Here that is 15 of 133 lines — and 49 of 81 substations, which is the bigger hole: the ≥45 kV substation distance was computed from 32 points. In a region with that many untagged features, the distance map is biased upwards; count them and say so.

Lines cross the region; substations are where a plant can actually connect.

Substations rather than lines

Why substations rather than lines

A plant cannot simply splice into a 132 kV line: it needs switchgear and protection — effectively a new substation — which is a large fixed cost and a long permitting process. Connecting to an existing substation with spare capacity avoids most of that. Distance to a line is a reasonable proxy for large plants that will build their own substation anyway.

Why capacity matters as much as distance

A substation that is already full cannot take another plant, whatever the distance. Grid operators publish connection capacity maps in many countries; where they exist, they should replace or weight the OSM-derived distance.

Why the distance transform is exact on a grid

distance_transform_edt computes the true Euclidean distance between cell centres, not a chamfer approximation. The only error is the half-cell from rasterisation — 12.5 m on a 25 m grid, negligible at kilometre scales.

Three decisions in five lines: several circuits, volts to kilovolts, and unknown kept as unknown.

Cables, planned lines and cross-border grids

Underground cables (power=cable) are often high voltage in cities; include them if the region is urban.
Substation polygons and nodes both occur; use centroids so each substation is one point.
Planned lines are tagged differently (construction:power=line or proposed:power=line) and are excluded by the query.
Cross-border grids: a line in another country is rarely a connection option; clip by jurisdiction where it matters.
Cost per km varies strongly with voltage and terrain; if you convert distance to money, make the rate an explicit parameter.

Internal links

How to build a siting exclusion mask in Python — the grid this distance lives on
How to build a solar suitability map — distance as a scored factor
Renewable siting constraints explained — exclusions against suitability
OpenStreetMap data model explained — tags, ways and relations
How to find the nearest point in GeoPandas — the vector alternative for points
How to measure distance to the nearest facility — distance along a network instead

FAQ

How do I calculate distance to the nearest power line in Python?

Rasterise the lines onto your analysis grid and compute scipy.ndimage.distance_transform_edt on the inverted mask, multiplied by the cell size.

Which voltage should I measure to?

The voltage the plant will connect at: distribution (15–45 kV) for small plants, transmission (132 kV and up) for utility-scale ones.

Is OSM power data complete?

In much of Europe the transmission grid is well mapped; distribution lines and voltage tags are patchier. Check coverage against the operator's maps.

Should I measure to lines or substations?

Substations for most plants; lines for large plants that will build their own substation.

Why is straight-line distance not enough?

Connection routes avoid settlements and protected land and often follow roads. Straight-line distance is a lower bound.

How do I handle lines with no voltage?

Count them and report them. Filters with >= drop them silently.