25 min read 5414 words

Building a Transmission Cost Surface Raster in NumPy

The scenario: a routing run produces a corridor that crosses a wetland, and the modeller’s first instinct is to raise the wetland cost. It gets raised from 50 to 500, the route still crosses, and at 5,000 the route finally goes around — through a residential subdivision. The cost surface, not the optimiser, is where routing goes wrong, and this page builds one that behaves. It is the input stage for grid routing and least-cost path analysis.

Root-cause analysis

Three modelling errors produce the behaviour above, and each has a specific fix.

  1. Exclusions encoded as large finite costs. Any finite cost is a price the optimiser is willing to pay when the detour is long enough. A wetland at 5,000 is not excluded — it is expensive, and the model will cross it rather than take a 40-kilometre detour. Exclusions have to be infinite, with a separate check afterwards that distinguishes “no route” from “expensive route”.
  2. Costs conflated with currency. Baking a dollar rate into the surface means every change to the cost estimate requires a re-route, and it hides the fact that the ratios between cells are what the optimiser actually uses. Relative multipliers, with the per-kilometre rate applied to the routed length afterwards, separate the two cleanly.
  3. Grids that do not align. A cost surface assembled from a 30-metre land-cover raster and a 10-metre DEM with different origins is a surface where each cell means two slightly different places. Routes then drift systematically toward the offset, which looks like a modelling preference and is an alignment bug.
Aligned cells versus a half-cell offset Two panels. The left panel shows a land-cover grid and a slope grid whose cell boundaries coincide exactly, with a routed path following the cheap cells. The right panel shows the same two grids offset by half a cell, so each combined cell draws its land cover from one place and its slope from another; the routed path drifts consistently toward the offset direction. Annotations give the offset as 15 metres on a 30 metre grid and note that the drift is systematic rather than random. Two rasters, one grid — or two grids pretending to be one aligned: same origin, same transform offset by half a cell (15 m on a 30 m grid) land cover cells slope cells The drift is systematic, so it survives averaging and looks like terrain following. Assert CRS, transform and shape across every band before a single cost is computed.

Pre-flight validation

Alignment is the check worth running first, because everything downstream inherits it.

python
import numpy as np
import rasterio


def assert_grids_align(paths: list[str], *, tol: float = 1e-6) -> dict:
    """Every band in a cost surface has to describe the same cells."""
    profiles = []
    for p in paths:
        with rasterio.open(p) as src:
            profiles.append(
                {"path": p, "crs": src.crs, "transform": src.transform,
                 "shape": (src.height, src.width), "nodata": src.nodata}
            )
    ref = profiles[0]
    for prof in profiles[1:]:
        if prof["crs"] != ref["crs"]:
            raise ValueError(f"{prof['path']}: CRS {prof['crs']} != {ref['crs']}")
        if prof["shape"] != ref["shape"]:
            raise ValueError(f"{prof['path']}: shape {prof['shape']} != {ref['shape']}")
        if not np.allclose(np.array(prof["transform"]), np.array(ref["transform"]), atol=tol):
            raise ValueError(f"{prof['path']}: affine transform differs from the reference grid")
    return ref

Fix implementation

python
import numpy as np
from rasterio.features import rasterize

BASELINE = 1.0   # grassland: every other multiplier is relative to this

LANDCOVER_COST = {
    11: np.inf, 21: 3.2, 22: 6.5, 23: np.inf, 31: 1.4,
    41: 2.4, 42: 2.6, 43: 2.5, 52: 1.2, 71: 1.0,
    81: 1.1, 82: 1.3, 90: np.inf, 95: np.inf,
}


def build_cost_surface(
    landcover: np.ndarray,
    slope_deg: np.ndarray,
    *,
    exclusions=(),
    crossings=(),
    transform=None,
    max_slope_deg: float = 25.0,
    slope_scale_deg: float = 10.0,
    crossing_cost: float = 1.5,
    construction_offset_m: float = 40.0,
) -> np.ndarray:
    """Per-cell relative build cost. Infinite where a line cannot go."""
    shape = landcover.shape
    cost = np.full(shape, BASELINE, dtype="float32")

    for code, factor in LANDCOVER_COST.items():
        cost[landcover == code] = factor

    # Slope raises cost quadratically: doubling the angle quadruples the multiplier.
    cost *= 1.0 + (slope_deg / slope_scale_deg) ** 2
    cost[slope_deg > max_slope_deg] = np.inf

    if len(exclusions):
        buffered = [g.buffer(construction_offset_m) for g in exclusions]
        blocked = rasterize([(g, 1) for g in buffered], out_shape=shape,
                            transform=transform, fill=0, dtype="uint8")
        cost[blocked == 1] = np.inf

    if len(crossings):
        cheap = rasterize([(g, 1) for g in crossings], out_shape=shape,
                          transform=transform, fill=0, dtype="uint8")
        cost = np.where(cheap == 1, np.minimum(cost, crossing_cost), cost)

    return cost

The crossing step deserves attention: it uses np.minimum rather than assignment, so an existing bridge over a river makes those cells cheap without also making them cheap where the bridge crosses a wetland. Assignment would punch a hole straight through an exclusion, which is the most common way a “corrected” surface produces an unbuildable route.

Relative build-cost multipliers by land cover A horizontal bar chart of relative cost multipliers against a grassland baseline of 1.0: existing corridor 0.7, grassland 1.0, pasture 1.1, shrub 1.2, cropland 1.3, barren 1.4, deciduous forest 2.4, evergreen forest 2.6, developed open space 3.2 and low-intensity development 6.5. Three classes are drawn as infinite rather than as bars: open water, wetlands and medium-intensity development. A note explains that the multipliers are relative, so the dollar rate is applied to the routed length afterwards. Multipliers, not currency — the rate is applied afterwards existing corridor 0.7× grassland (baseline) 1.0× pasture 1.1× shrub 1.2× cropland 1.3× barren 1.4× deciduous forest 2.4× evergreen forest 2.6× developed, open space 3.2× developed, low intensity 6.5× open water — ∞ wetlands — ∞ developed, medium+ — ∞ Infinity, not a large number: any finite cost is a price the optimiser will pay when the detour is long enough — which is exactly how an "excluded" wetland ends up with a line across it.

Fallback routing and performance tuning

  • Store the surface as float32, not float64. Routing is memory-bound on large corridors, and the precision beyond float32 is meaningless for a relative multiplier.
  • Clip to a corridor buffer before assembling. A national surface is never needed; a buffer three to five kilometres either side of the straight line removes about 90 percent of the cells.
  • Keep the physical layers cached and the weights per project. Land cover, slope and hydrography are regional and slow to prepare; the multipliers and exclusions are per project and fast to apply.
  • Represent infinity honestly in storage. GeoTIFF cannot hold np.inf in every dtype — write a companion uint8 exclusion mask and reconstruct the infinities on read.
  • Run a sensitivity pass, not a single surface. Two plausible weightings that produce the same corridor are a robust answer; two that diverge tell you which weight the study actually hinges on.

Downstream validation

python
def assert_cost_surface(cost: np.ndarray, *, origin_rc, dest_rc) -> None:
    """Four properties a usable cost surface must have."""
    assert cost.dtype == np.float32, "use float32 — routing is memory bound"
    assert np.isfinite(cost).any(), "every cell is excluded — the mask swallowed the corridor"
    assert np.nanmin(cost[np.isfinite(cost)]) > 0, "zero or negative cost lets a route loop for free"
    for name, rc in (("origin", origin_rc), ("destination", dest_rc)):
        if not np.isfinite(cost[rc]):
            raise ValueError(f"{name} cell is excluded — snap it or override before routing")
Which weights the answer actually depends on A table of four weight perturbations with their effect on the routed corridor. Doubling the forest multiplier moves the corridor 3.1 kilometres and changes routed length by 4.0 percent. Halving the crossing cost moves it 0.0 kilometres and changes length by 0.3 percent. Removing the existing-corridor discount moves it 6.8 kilometres and changes length by 9.2 percent. Raising the slope exponent moves it 0.4 kilometres and changes length by 1.1 percent. The two perturbations that move the corridor are flagged as the weights the study depends on, and the other two as noise. Perturb each weight; watch the corridor, not the cost perturbation corridor moves length changes forest multiplier 2.4 → 4.8 3.1 km 4.0% crossing cost 1.5 → 0.75 0.0 km 0.3% existing-corridor discount removed 6.8 km 9.2% slope exponent 2 → 2.5 0.4 km 1.1% A weight whose corridor does not move is a weight the study does not depend on — report it and stop arguing about it. Two weightings that produce different corridors both belong in the deliverable.

Storing and versioning the surface

A cost surface is an artefact, not a scratch array, and it should be written with enough metadata to be re-used and challenged. Three things belong in the file: the weight table that produced it, the source layers with their vintages, and the exclusion buffer applied. GeoTIFF and Zarr both carry arbitrary key-value metadata, so none of this needs a sidecar that can be separated from the raster.

Infinity is the one storage awkwardness. Most raster formats cannot hold np.inf in a float32 band in a way every reader honours, so the durable pattern is two bands: a finite cost band with exclusions written as the maximum finite value, and a uint8 exclusion mask. On read, the mask restores the infinities. Doing it the other way round — writing a sentinel like −9999 and hoping every consumer knows — is how an exclusion silently becomes a cheap cell in someone else’s pipeline.

Version the surface by content rather than by date. Hashing the weight table together with the source layer vintages gives a short identifier that changes exactly when something that matters changed, and recording that identifier on every route makes a corridor traceable to the surface that produced it. Two routes with different surface identifiers are not comparable, however similar they look.

Frequently asked questions

Why must the minimum cost be strictly positive?

Because a zero-cost cell is free to traverse, and a connected region of them lets the optimiser wander at no cost — which produces routes with pointless meanders that all have the same total cost. A baseline of 1.0 for the cheapest land keeps every step priced and makes the shortest of several equal-cost routes the one that wins.

Should slope raise cost linearly or quadratically?

Quadratically, because construction cost does. Access-road switchbacks, pad cut-and-fill and structure spotting all get disproportionately harder as the ground steepens, and a linear multiplier under-prices the steep ground that actually decides the route. The exact exponent matters less than the shape; what matters most is the hard cut-off at the crane or construction limit.

How should water crossings be priced when there is no existing bridge?

As a finite, large adder rather than an exclusion, encoded as a narrow band of expensive cells across the water rather than as a blanket cost on the whole water body. That way the optimiser chooses the narrowest sensible crossing, which is what a routing engineer does, instead of treating the entire river as uniformly expensive and crossing at an arbitrary point.

Can the surface include a preference for existing corridors?

Yes, and it is one of the highest-value weights available. Give existing transmission and pipeline corridors a multiplier below the grassland baseline — 0.6 to 0.8 is a common range — and routes will follow them wherever the detour is modest, which reflects both the easement saving and the permitting preference.

What resolution should the surface be?

Thirty metres for the corridor search and ten metres or better where crossings are chosen. A uniform fine grid buys almost nothing in open terrain and costs an order of magnitude in cells; the two-pass approach described in the parent page gets the accuracy where it matters at a fraction of the run time.

How do I know a weight change actually mattered?

Run both and difference the accumulated-cost surfaces, not just the routes. Two weightings that produce visibly different corridors but nearly identical costs are within the model’s own noise, and the honest report presents both. A weight that moves the total cost by more than the difference between the top two candidate routes is a weight the study depends on, and it belongs in the sensitivity table.

Should the surface be rebuilt when a new DEM is published?

Only with a comparison. A newer DEM is usually finer and more accurate, and both properties change the slope band — sometimes enough to move a corridor. Rebuild, route both surfaces, and report the difference rather than silently replacing the old answer, because a route that appears in a study is a claim tied to the data that produced it.

How should nodata cells be treated?

As excluded, and counted. A nodata hole in the land-cover raster is not cheap land, and leaving it at the baseline multiplier is exactly how a route ends up running through the one area nobody has mapped. Count the nodata cells inside the corridor and report the fraction: above a percent or two, the surface needs a better input rather than a better weight.

Related articles

Back to Grid Routing & Least-Cost Path Analysis