26 min read 5769 words

Grid Routing & Least-Cost Path Analysis

Least-cost path analysis is where a straight-line screen becomes a buildable corridor, and it is the stage that sits between proximity screening and a real interconnection estimate in the grid infrastructure and network proximity analysis pipeline. The failure mode it addresses is not that straight-line distance is wrong — it is a perfectly good lower bound — but that it is used as though it were the answer. Across sited interconnections the ratio of routed to straight-line length has a median near 1.28 and a tail beyond 1.9, and the projects in that tail are exactly the ones whose economics a screen declared healthy.

Routing turns three implicit assumptions into explicit ones. A cost surface says what the terrain actually costs per metre rather than treating all land as equal. A crossing penalty says what a river, a rail line or an interstate costs to span rather than pretending they are free. And a hard exclusion says where a route cannot go at all, which is the difference between an expensive corridor and a non-existent one.

Four inputs, one per-cell cost, and where the infinities go Four stacked input panels feeding one output. The first, land cover, lists relative multipliers: grassland 1.0, pasture 1.1, cropland 1.3, shrub 1.2, forest 2.4 to 2.6, developed 3.2 to 6.5, and water and wetlands as infinite. The second, slope, shows a quadratic multiplier rising from 1.0 at flat ground to 7.25 at 25 degrees, then infinity above the build limit. The third, vector exclusions, shows protected areas buffered by a 40 metre construction offset and burned in as impassable. The fourth, existing crossings, shows bridge and easement cells forced to a low cost of 1.5. All four combine into a single per-cell cost raster. The cost surface is the model — everything after it is mechanical land cover grass 1.0 · crop 1.3 · forest 2.4 water, wetland → ∞ slope 1 + (slope/10)² multiplier above 25° → ∞ vector exclusions buffered 40 m for working room burned in as ∞ existing crossings bridges, easements forced to 1.5 per-cell relative cost costs are multipliers, never currency — the $/km rate is applied to the routed length exclusions are ∞, not a large number: a finite cost lets the optimiser cross Every argument about a route is an argument about these weights.

Why the straight line and the route diverge

Three mechanisms produce the gap, and they compound rather than average out.

The first is avoidance. Protected areas, dense settlement and open water are not merely expensive — they are excluded, so a route that would cross them must go around, and the detour scales with the size of the obstacle rather than with its cost. A single designated corridor lying across the direct line can add several kilometres to a ten-kilometre route.

The second is crossing infrastructure. Rivers, railways and controlled-access highways can be crossed, but only at a cost and often only at specific points where an easement already exists. A routing model that treats a river as uniformly expensive produces a route that crosses at the cheapest cell; a model that knows about existing crossings produces the route an engineer would actually build, which is usually longer and cheaper.

The third is terrain. Slope raises construction cost non-linearly — access roads, pad preparation and structure spotting all get harder — and above the crane specification it stops being a cost and becomes an exclusion, exactly as it does in environmental constraint and exclusion screening.

Prerequisites and data requirements

The workflow assumes Python 3.11+ with rasterio>=1.3, numpy, scikit-image>=0.22 (for route_through_array) and geopandas>=0.14. Inputs are a DEM, a land-cover raster, vector layers for exclusions and crossings, and the origin and destination points — usually a project substation location and a point of interconnection.

Two requirements are structural. Every raster must share one grid: same CRS, same affine transform, same shape, because a cost surface built from bands that disagree by half a pixel produces routes that drift systematically toward the offset. And the working frame must be projected and metric, so that a cell’s cost can be expressed per metre and a route length can be read directly off the path.

Cell size is the parameter that decides everything else. A 30-metre grid over a 40-kilometre corridor is about 1.8 million cells and routes in under a second; a 5-metre grid over the same corridor is 64 million cells and needs tiling. Route coarse to find the corridor, then refine inside a buffer around it — the two-pass approach is both faster and more accurate than either resolution alone.

Core implementation: building the cost surface

The cost surface is the model. Everything downstream is mechanical, and every argument about a route is really an argument about the weights below.

python
import numpy as np
import rasterio
from rasterio.features import rasterize

# Cost multipliers are per-metre relative costs, not currencies: the absolute
# figure comes from the $/km estimate applied to the routed length afterwards.
LANDCOVER_COST = {
    11: np.inf,   # open water — excluded
    21: 3.2,      # developed, open space
    22: 6.5,      # developed, low intensity
    23: np.inf,   # developed, medium intensity — excluded in practice
    41: 2.4,      # deciduous forest — clearing cost
    42: 2.6,      # evergreen forest
    52: 1.2,      # shrub
    71: 1.0,      # grassland — the baseline
    81: 1.1,      # pasture
    82: 1.3,      # cultivated crops — easement cost
    90: np.inf,   # woody wetlands — excluded
    95: np.inf,   # emergent wetlands — excluded
}


def build_cost_surface(
    landcover: np.ndarray,
    slope_deg: np.ndarray,
    exclusions: list,
    crossings: list,
    transform,
    shape: tuple[int, int],
    *,
    max_slope_deg: float = 25.0,
) -> np.ndarray:
    """Per-cell relative cost of building a transmission line through that cell."""
    cost = np.ones(shape, dtype="float32")

    # 1. Land cover sets the base cost and the first set of exclusions.
    for code, factor in LANDCOVER_COST.items():
        cost[landcover == code] = factor

    # 2. Slope raises cost quadratically, then excludes above the build limit.
    slope_factor = 1.0 + (slope_deg / 10.0) ** 2
    cost *= slope_factor
    cost[slope_deg > max_slope_deg] = np.inf

    # 3. Vector exclusions are burned in as impassable.
    if exclusions:
        blocked = rasterize(
            [(geom, 1) for geom in exclusions],
            out_shape=shape, transform=transform, fill=0, dtype="uint8",
        )
        cost[blocked == 1] = np.inf

    # 4. Existing crossings punch cheap holes through linear barriers.
    if crossings:
        cheap = rasterize(
            [(geom, 1) for geom in crossings],
            out_shape=shape, transform=transform, fill=0, dtype="uint8",
        )
        cost[cheap == 1] = np.minimum(cost[cheap == 1], 1.5)

    return cost

Two design choices in that function are worth defending. Costs are relative multipliers rather than currency, because the dollar figure belongs at the end — applied to the routed length as a per-kilometre rate — and keeping it out of the surface stops a change in the cost estimate from requiring a re-route. And exclusions are np.inf rather than a large finite number, because a large finite cost lets the optimiser cross an exclusion when the detour is long enough, which produces a route that is optimal and unbuildable.

Running the route

With a surface in hand, the route itself is a few lines. skimage.graph.route_through_array runs Dijkstra over the array with optional diagonal movement and geometric weighting, which matters: without it a diagonal step is counted as one cell rather than 1.414, and routes acquire a characteristic staircase bias.

python
from skimage.graph import route_through_array
from shapely.geometry import LineString


def least_cost_route(
    cost: np.ndarray,
    transform,
    origin_rc: tuple[int, int],
    dest_rc: tuple[int, int],
) -> tuple[LineString, float, float]:
    """Return the route geometry, its length in metres and its accumulated cost."""
    finite = np.where(np.isinf(cost), np.float32(1e9), cost)
    indices, weight = route_through_array(
        finite, origin_rc, dest_rc, fully_connected=True, geometric=True
    )
    if weight >= 1e9:
        raise ValueError("no traversable route — origin and destination are separated by exclusions")

    xs, ys = zip(*[transform * (c + 0.5, r + 0.5) for r, c in indices])
    line = LineString(zip(xs, ys))
    return line, float(line.length), float(weight)

The 1e9 substitution is a deliberate compromise: route_through_array cannot handle infinities, so exclusions become a cost so large that any traversable alternative wins, and the returned weight is checked afterwards to distinguish “expensive route” from “no route”. Silently returning a route that crosses a wetland because the alternative was longer is the failure this check exists to prevent.

Straight line, naive route, and the route that knows where the bridge is A plan view with the project substation at the lower left and the point of interconnection at the upper right. A shaded wetland lies between them and a river runs north to south with a single marked bridge. Three paths are drawn: a dashed straight line of 8.2 kilometres crossing both the wetland and the river at an arbitrary point; a route of 15.9 kilometres that avoids the wetland but crosses the river away from the bridge; and a route of 14.6 kilometres that avoids the wetland and crosses at the bridge. A table beside the plan gives each path its length and circuity factor. One origin, one destination, three answers river protected wetland existing bridge project POI straight line 8.2 km circuity 1.00 not buildable — crosses both route, no crossing data 15.9 km circuity 1.94 fords the river at a random cell route with crossings 14.6 km circuity 1.78 uses the bridge — the real corridor A finite cost on the wetland would have produced the straight line again, at a price the optimiser was willing to pay — which is why exclusions are infinite and checked afterwards.

Error handling and edge cases

No traversable route. This is a real answer, not an error condition — a site behind a continuous exclusion has no corridor at the given constraints — and it should propagate as a result with the blocking geometry named, so a developer can see whether relaxing one constraint opens a path.

Origin or destination inside an exclusion. Common, and usually a data artefact: a substation polygon overlapping a developed land-cover class, or a point snapped to the wrong cell. Force the endpoint cells to a finite cost before routing and record that the override happened, rather than letting the route fail with an unhelpful message.

A route that hugs an exclusion boundary. Optimal and often unbuildable, because construction needs working room. Buffer exclusions by a construction offset — 30 to 50 metres is typical — before burning them into the surface, so the optimiser keeps its distance without any special-case logic in the routing step.

Diagonal staircase artefacts. If routes look like a flight of stairs rather than a corridor, the geometric weighting is off or the surface is too coarse relative to the cost variation. Simplify the resulting line with a tolerance of roughly one cell before reporting its length, but compute the cost from the unsimplified path.

Performance and scalability

Dijkstra over a raster is O(n log n) in the cell count, so cell count is the only lever that matters. Three techniques, in the order they pay off. Clip the surface to a corridor buffer around the straight line — three to five kilometres either side is generous for most interconnections — which typically removes 90 percent of the cells before any routing happens. Route coarse then refine, using a 30-metre pass to find the corridor and a 10-metre pass inside a buffer around it. And route each candidate independently in a worker pool, since routes share nothing but the read-only surface.

For a portfolio screen, the useful shortcut is not to route at all until the shortlist exists: run the straight-line screen over everything, take the top candidates by straight-line distance and capacity, and route those. Routing everything is the most common reason a screening pipeline that worked on one county does not finish on a portfolio.

Cell count is the only lever: clip, then refine A chart of four routing configurations by cell count and wall-clock time. Thirty metre cells over the full extent: 4.9 million cells, 6.2 seconds. Thirty metre cells clipped to a 4 kilometre corridor buffer: 0.6 million cells, 0.7 seconds. Ten metre cells over the full extent: 44 million cells, 71 seconds. Ten metre cells inside a buffer around the coarse route: 2.4 million cells, 3.1 seconds. The last configuration is marked as the working choice, being twenty times faster than the full fine pass at the same effective resolution where it matters. Two passes beat one fine grid, by a factor of twenty 30 m · full extent 6.2 s · 4.9M cells 30 m · corridor buffer 0.7 s · 0.6M cells 10 m · full extent 71.0 s · 44.0M cells 10 m · around coarse route 3.1 s · 2.4M cells wall-clock to route one 40 km corridor Clip first: a corridor buffer removes about 90% of the cells before routing starts Then refine: the fine pass only needs to run where the coarse route already goes

Validation and audit trail

A route is a claim about buildability, so the record has to carry enough to defend it: the cost surface version and the weights that produced it, the exclusion layers and their buffers, the cell size of each pass, the routed length, the accumulated cost, the crossings used, and the straight-line length for comparison. The ratio between the last two is the circuity factor, and publishing it makes the difference between the screen and the route explicit rather than surprising.

Three assertions belong in CI. The routed length must be greater than or equal to the straight-line length, which catches a coordinate or transform error. The route must not intersect any exclusion geometry, which catches the finite-cost-substitution bug directly. And the accumulated cost must be finite, which catches the case where the only path found crosses a barrier.

Reading the accumulated-cost surface

Dijkstra produces more than a path. The accumulated-cost array it fills in on the way is a complete answer to “what would it cost to reach every cell from this origin”, and reading it directly is often more useful than the single route extracted from it.

Three questions fall out of that array for free. Comparing several candidate points of interconnection needs one pass rather than one per destination, because the cost to each is simply the value of the array at that cell. Drawing iso-cost contours over the array shows the shape of the accessible region — where the cheap corridors run, and where an obstacle splits the surface into two basins that only connect a long way round. And subtracting the straight-line distance from the accumulated cost highlights exactly where the terrain is expensive, which is the map a routing engineer wants when deciding whether to challenge a constraint.

The array is also the honest place to express uncertainty. Running the same origin over two plausible weightings and differencing the two accumulated-cost surfaces shows which parts of the study area are robustly cheap and which are cheap only under one set of assumptions. A route drawn through the second kind of region is a route that will be re-litigated.

Should the accumulated-cost surface be published with the route?

For any study that supports a decision, yes. The route is one line through a surface, and a reviewer who disagrees with it is really disagreeing with the surface. Publishing both — as a raster with its weights recorded — turns an argument about a line into an argument about assumptions, which is the one that can actually be resolved.

Frequently asked questions

Should the cost surface include existing corridors as cheap cells?

Yes, and it is one of the highest-value additions. Co-locating with an existing transmission or pipeline corridor avoids new easements, follows land already cleared, and is usually favoured by permitting authorities. Give existing corridors a cost below the grassland baseline and the optimiser will follow them wherever the detour is modest — which is what a routing engineer does by hand.

How sensitive is the route to the weights?

Very, in the corridor it chooses, and much less in the length it reports. Doubling the forest multiplier can move a route several kilometres sideways while changing the routed length by a few percent, because the alternatives are similar in length and different in character. This is why the weights belong in versioned configuration and why a sensitivity run — the same route under two plausible weightings — is worth more than a single precise-looking answer.

Can least-cost paths handle multiple destinations?

Yes, and it is cheaper than it looks: a single Dijkstra from one origin yields the accumulated-cost surface to every cell, so the cost to each of several candidate points of interconnection comes from one pass. Use that when comparing interconnection options for one project. Routing many origins to one destination is the same problem reversed, and the same trick applies.

What resolution should the final route use?

Fine enough to resolve the crossings that decide it, which in practice means 10 metres or better near rivers, railways and highways, and 30 metres elsewhere. A uniform fine grid over the whole corridor buys almost nothing: the route’s length is insensitive to resolution in open terrain and highly sensitive to it where a crossing point is being chosen.

How does routed distance feed the interconnection cost estimate?

As the length term in a per-kilometre rate, with the crossings and terrain classes carried through as adders rather than folded into the length. A route reported as “14.6 kilometres, two river crossings, 3.1 kilometres above 15 percent slope” supports a cost estimate that a reviewer can challenge line by line; a route reported as “14.6 kilometres” invites a single blended rate that hides all of it. The straight-line figure from proximity and distance calculations remains the screen; the routed figure is what the estimate uses.

Does the same approach work for collection systems inside a project?

Partly. The cost-surface idea carries over directly, but the problem changes from a single path to a minimum-cost tree connecting many turbines to one substation, which is a Steiner-tree problem rather than a shortest-path one. A practical approximation is to route each turbine to the substation over the same surface and then merge shared segments, which is not optimal but is close enough for layout screening and uses exactly the machinery above.

How should the route handle land the developer does not control?

As a cost, not an exclusion, unless the landowner has formally refused. Easement acquisition is expensive and slow, and a routing model that treats every uncontrolled parcel as impassable will report that no corridor exists for almost every project. The workable encoding is a parcel-level multiplier that rises with the number of distinct owners crossed, which favours routes along fewer, larger holdings — the same preference a land agent expresses.

Does the cost surface need to be rebuilt for every project?

The physical layers do not; the weights and exclusions usually do. Terrain, land cover and hydrography are regional and can be prepared once and cached, which is most of the build time. What changes per project is the slope limit set by the crane, the buffer applied for working room, and which constraint classes are treated as hard. Separating the two halves — a cached physical surface and a per-project weighting pass — keeps a re-route to seconds rather than minutes.

What happens when two candidate routes are within a few percent of each other?

Report both. A cost surface is a model with uncertain weights, and two routes separated by less than the weight uncertainty are not distinguishable by the model, however precisely the optimiser ranks them. Presenting the pair with the properties that differ — number of owners, crossings, forest kilometres — moves the decision to the criteria a routing engineer would use, which is where it belongs.

Can the surface encode permitting difficulty rather than construction cost?

Yes, and it often should. Nothing in the method requires the multiplier to represent dollars: a surface weighted by expected permitting duration produces the route that reaches energisation soonest, which is not always the cheapest one to build. Running both and comparing the corridors is one of the more useful outputs this stage can produce, because the two objectives diverge most exactly where a project is most at risk.

Explore this section

Back to Grid Infrastructure & Network Proximity Analysis