17 min read 3818 words

Computing Least-Cost Interconnection Routes with scikit-image

The scenario: route_through_array returns a route, its length is plausible, and a reviewer notices it clips the corner of a designated wetland for 180 metres. The optimiser did exactly what it was asked. This page covers the mechanics that sit between a cost surface and a defensible corridor, and it is the execution half of grid routing and least-cost path analysis.

Root-cause analysis

Four mechanical faults account for nearly every bad route that comes out of an otherwise correct surface.

  1. Infinity substitution without a check. route_through_array cannot traverse np.inf, so the surface is usually converted to a large finite value. Without checking the returned weight afterwards, “there is no route” becomes “here is a route through the exclusion”.
  2. Geometric weighting left off. With geometric=False a diagonal step costs the same as an orthogonal one, so diagonal movement is under-priced by a factor of 1.414 and routes acquire a staircase bias that is easy to mistake for terrain following.
  3. Endpoints inside excluded cells. A substation polygon overlapping a developed land-cover class puts the destination in an infinite-cost cell, and the failure message is unhelpful.
  4. Cell-centre geometry. Converting indices to coordinates without the half-cell offset shifts the whole route by half a pixel, which is invisible at national scale and matters when the route is compared against a parcel boundary.
Diagonal steps priced at 1, or at the square root of two Two grids with the same origin and destination. In the first, with geometric weighting off, the route takes a staircase of alternating orthogonal steps whose total is reported as 12 cells while the true ground distance is 17 cell-lengths, a 41 percent understatement. In the second, with geometric weighting on, diagonal steps are priced at 1.414 and the route runs straight to the destination, reporting the true distance. A note records that the staircase is easy to mistake for terrain following. geometric=False under-prices every diagonal by 41% geometric=False — staircase 12 cells reported · 17 cell-lengths walked geometric=True — true diagonal 8.5 cell-lengths, reported correctly The staircase looks like terrain following and is a pricing artefact. fully_connected=True allows the diagonal; geometric=True is what makes it cost what it should.

Pre-flight validation

python
import numpy as np


def prepare_for_routing(
    cost: np.ndarray,
    origin_rc: tuple[int, int],
    dest_rc: tuple[int, int],
    *,
    blocked_value: float = 1e9,
) -> tuple[np.ndarray, dict]:
    """Substitute infinities, rescue endpoints, and report what was changed."""
    report = {"blocked_cells": int(np.isinf(cost).sum()), "endpoint_overrides": []}
    finite = np.where(np.isinf(cost), np.float32(blocked_value), cost).astype("float32")

    for name, rc in (("origin", origin_rc), ("destination", dest_rc)):
        if finite[rc] >= blocked_value:
            # Rescue rather than fail: take the cheapest finite cell in a 5x5 window.
            r, c = rc
            window = finite[max(r - 2, 0):r + 3, max(c - 2, 0):c + 3]
            if not np.isfinite(window).any() or window.min() >= blocked_value:
                raise ValueError(f"{name} is inside a large excluded region — move it or relax a constraint")
            finite[rc] = float(window.min())
            report["endpoint_overrides"].append(name)
    return finite, report

Fix implementation

python
import numpy as np
from shapely.geometry import LineString
from skimage.graph import route_through_array


def least_cost_route(
    cost: np.ndarray,
    transform,
    origin_rc: tuple[int, int],
    dest_rc: tuple[int, int],
    *,
    blocked_value: float = 1e9,
    cell_size_m: float | None = None,
) -> dict:
    """Route, then prove the route is buildable before returning it."""
    finite, report = prepare_for_routing(cost, origin_rc, dest_rc, blocked_value=blocked_value)

    indices, weight = route_through_array(
        finite, origin_rc, dest_rc,
        fully_connected=True,     # allow diagonal movement
        geometric=True,           # and price it at sqrt(2), not 1
    )
    if weight >= blocked_value:
        raise ValueError("no traversable route: every path crosses an excluded region")

    # Half-cell offset puts the vertex at the cell centre, where the cost applies.
    coords = [transform * (c + 0.5, r + 0.5) for r, c in indices]
    line = LineString(coords)

    cs = cell_size_m or abs(transform.a)
    return {
        "geometry": line,
        "length_m": float(line.length),
        "accumulated_cost": float(weight),
        "cells": len(indices),
        "mean_cost_per_cell": float(weight) / max(len(indices), 1),
        "straight_line_m": float(
            LineString([coords[0], coords[-1]]).length
        ),
        "circuity": float(line.length) / max(
            LineString([coords[0], coords[-1]]).length, 1e-9
        ),
        "cell_size_m": cs,
        **report,
    }

Returning the circuity factor alongside the length is what connects this stage back to the straight-line screen in proximity and distance calculations: a route with a circuity of 1.9 is a project whose screen was optimistic, and that fact should travel with the number.

One accumulated-cost pass serves every candidate destination A diagram of one origin and five candidate points of interconnection over a shaded accumulated-cost surface, with iso-cost contours radiating from the origin. Each candidate is annotated with the cost read directly from the array: 14.6, 18.2, 21.9, 24.4 and 31.7 kilometres of cost-weighted distance. A comparison beside it gives the wall clock: five independent routing runs take 28.4 seconds, while one find_costs pass plus five array lookups takes 6.1 seconds and produces the same ranking. find_costs once, read every destination project 14.6 km 18.2 km 21.9 km 24.4 km 31.7 km five independent routes 28.4 s one find_costs + 5 lookups 6.1 s Identical answers — the array is the same computation, kept instead of discarded Trace back only the routes you intend to draw: the ranking comes from the array, the geometry from the traceback.

Fallback routing and performance tuning

  • One Dijkstra, many destinations. route_through_array discards the accumulated-cost array, but skimage.graph.MCP_Geometric exposes it: run find_costs once from the origin and read the cost to each candidate point of interconnection, then traceback only the ones worth drawing.
  • Route coarse, then refine. A 30-metre pass to find the corridor and a 10-metre pass inside a buffer around it is roughly twenty times faster than a single fine pass and more accurate where it matters.
  • Simplify the output line, not the cost. A raster route has a vertex per cell; simplifying with a tolerance of about one cell removes the staircase without changing the corridor. Compute the cost from the unsimplified path.
  • Watch memory on wide corridors. The MCP structures are several arrays the size of the surface; at 44 million cells that is gigabytes, which is the practical reason for the corridor clip.

Downstream validation

python
import geopandas as gpd


def assert_route_is_buildable(route: dict, exclusions: gpd.GeoSeries, *, tol_m: float = 1.0) -> None:
    """The four assertions that separate an optimal route from a buildable one."""
    line = route["geometry"]
    assert route["accumulated_cost"] < 1e9, "route traverses an excluded region"
    assert route["length_m"] >= route["straight_line_m"] - tol_m, (
        "routed length below the straight line — a transform or CRS error"
    )
    hit = exclusions[exclusions.intersects(line)]
    assert hit.empty, f"route intersects {len(hit)} exclusion geometries despite a finite cost"
    assert route["circuity"] < 4.0, (
        f"circuity {route['circuity']:.2f} — the corridor is almost certainly blocked, not merely expensive"
    )
Four assertions and the bug each one catches A four-row table pairing an assertion with the failure it catches. Accumulated cost below the blocked sentinel catches a route that crossed an excluded region. Routed length at least the straight-line length catches a coordinate or affine transform error. No intersection with any exclusion geometry catches a sentinel value set too low relative to the corridor length. A circuity factor under four catches a destination that is unreachable under the current constraints rather than simply far away. The route is a claim about buildability — check it accumulated_cost < blocked_value a route straight through an exclusion length_m >= straight_line_m a coordinate or transform error not route.intersects(exclusions) a sentinel value set too low circuity < 4.0 a destination that is unreachable, not far The third assertion is the one that matters most: it checks the geometry, not the arithmetic that produced it.

What the accumulated-cost array is worth on its own

route_through_array throws away the most useful thing it computes. MCP_Geometric.find_costs keeps it: an array holding the cost of reaching every cell from the origin, which answers several questions the route alone cannot.

The first is comparison. Five candidate points of interconnection cost five array lookups rather than five routing runs, and the ranking is exact rather than approximate because all five costs come from the same pass. The second is shape. Contouring the array shows where the cheap corridors run and where an obstacle splits the surface into two basins that only connect a long way round — which is the map a routing engineer wants when deciding whether a constraint is worth challenging. The third is uncertainty: differencing two accumulated-cost arrays computed under two plausible weightings shows which parts of the study area are robustly cheap and which are cheap only under one set of assumptions.

Keeping the array costs nothing beyond memory, since it was computed either way. Publishing it alongside the route turns an argument about a line into an argument about the surface, which is the one that can be settled.

Frequently asked questions

Why does the route hug the edge of an exclusion?

Because the cheapest traversable cells are the ones immediately outside it, and nothing in the model says a line needs working room. Buffer the exclusions by a construction offset before rasterising — 30 to 50 metres is typical — and the optimiser keeps its distance without any special-case logic.

Should fully_connected ever be False?

Only when the movement model genuinely forbids diagonals, which for a transmission line it does not. With diagonals disabled every route becomes a staircase of orthogonal steps whose length is systematically overstated by up to 41 percent on diagonal runs.

How do I route to several substations at once?

Use MCP_Geometric.find_costs from the project location, which fills the accumulated-cost array for the whole surface in one pass, then read the cost at each substation cell. Only trace back the routes you intend to draw. For five candidate points of interconnection this is roughly five times faster than five independent routes, and the comparison is exact rather than approximate.

What does a very high circuity factor mean?

Usually that the destination is effectively unreachable under the current constraints rather than merely expensive. A circuity above about three says the optimiser is taking a long way round an obstacle that spans the direct line, and the useful output is the name of that obstacle — which comes from intersecting the route’s bounding corridor with the exclusion layers rather than from the route itself.

Should the routed geometry be smoothed?

Simplified, not smoothed. Simplification with a one-cell tolerance removes the raster staircase while keeping every vertex on the routed path; smoothing with a spline moves vertices off it, which can push the line into a cell the model excluded. Simplify for presentation, keep the raw path for the cost.

Can this handle a route that must pass through a waypoint?

Yes — route origin to waypoint and waypoint to destination, then concatenate. The concatenation is exact because the accumulated cost is additive, and it is the standard way to honour a landowner agreement or a mandated crossing point without distorting the cost surface to force the outcome.

How long should a route take to compute?

Under a second for a clipped 30-metre corridor, and a few seconds for a refined 10-metre pass inside a buffer. A route that takes minutes is almost always running over an unclipped surface, and the fix is the corridor buffer rather than a faster machine. Wall-clock is a useful smoke test for exactly that reason: a sudden increase usually means the clip stopped working, not that the terrain changed.

Can the same code route a distribution feeder or an access road?

Yes — only the weights change. An access road cares about slope far more and about land cover far less, and a distribution feeder can use narrower corridors and cross land a transmission line cannot. The machinery is identical, which is a good argument for keeping the weights in configuration where a second profile is a file rather than a fork.

Related articles

Back to Grid Routing & Least-Cost Path Analysis