Evacuation Routing & Road Network Analysis

Problem Framing

A fast-moving grass fire jumps a ridge at 14:20 and the sheriff orders a mandatory evacuation of three zones on the wildland-urban interface. Incident command needs to push turn-by-turn guidance to residents and to the deputies running door-to-door, and it needs to do so against a road network that is changing minute by minute: a low-water bridge is already under water, a two-lane county road has been converted to contraflow to double outbound capacity, and the fire perimeter itself is advancing across the only paved route out of the northern zone. If the routing engine sends a column of vehicles toward a shelter that sits on the wrong side of the fire, or up a one-way street into oncoming apparatus, the map has actively made the incident more dangerous. Evacuation routing exists to make that class of error structurally impossible: it computes least-time paths on a graph that encodes real, currently-open road segments, applies live hazard closures the instant a perimeter update lands, and records every decision so the routes remain defensible in the after-action review. This page specifies that engine as runnable Python and connects it to the wider Incident Mapping & Multi-Agency Sync Workflows contract that feeds it live perimeters and shelter status.

Prerequisites

This workflow assumes a senior engineer is comfortable with graph algorithms and the Python geospatial stack, and that the following preconditions hold before the first route is computed:

• Packages: networkx >= 3.0 for the directed graph and shortest-path search, shapely >= 2.0 for geometry and the STRtree spatial index, and pyproj >= 3.4 for reprojection. igraph or pandana are drop-in alternatives when the network reaches millions of edges and pure-Python Dijkstra becomes the bottleneck; the modelling here maps cleanly onto either.
• A projected, metre-based CRS. Every length, hazard buffer, and isochrone radius in this page is expressed in metres, which is only meaningful in a projected reference system. Route in the local Universal Transverse Mercator zone and treat geographic degrees as a display format only. The datum-aware reprojection is owned by the Coordinate Reference Systems for Disaster Zones workflow; this stage assumes the network and every hazard polygon already share one projected CRS.
• A clean road-segment source. Each segment must carry a stable identifier, from-node and to-node, a projected LineString, a length in metres, a free-flow speed, and a oneway flag. Whether the source is OpenStreetMap via osmnx, a county centreline file, or an agency roads geodatabase, the topology must be noded — segments that visually cross without sharing a node are two separate, unconnected edges and will silently break routing.
• A live hazard feed. Fire and flood perimeters arrive as polygons on a message backbone; the transport and delivery guarantees are handled by the WebSocket & MQTT for Live Incident Feeds layer. This stage consumes a validated Polygon and applies it; it does not manage the subscription.
• A shelter and assembly-point catalogue. Destinations must be snapped to graph nodes and carry live open/closed and capacity status, aligned to the shelter capacity and resource tracking schemas so the router never sends evacuees to a facility that is already full.

Routing Architecture

Evacuation routing is a staged transformation from a static base network to a hazard-aware, capacity-aware set of routes and reachability surfaces. The base graph is built once and kept resident in memory. Each perimeter update from the hazard feed is applied as a fast spatial overlay that closes or penalizes only the intersecting edges, leaving the rest of the multi-megabyte graph untouched. Routing then runs a single least-time search from each evacuation-zone origin to a virtual sink attached to every open shelter, so the algorithm — not a hand-coded rule — selects the destination that is genuinely fastest to reach on the roads that are still open. The same graph answers a second question the incident commander always asks: not just which way out, but which areas can clear within the window, which is an isochrone computed by a truncated shortest-path search. The diagram below traces a single origin-to-shelter route around a flood perimeter that has closed the three central segments.

Step-by-Step Implementation

Step 1 — Build a routable directed road graph

Routing correctness begins with the graph. A two-way street is two directed edges; a one-way street is one; and a contraflow conversion temporarily opens the reverse lane of a one-way segment at a governed speed. Modelling these as distinct directed edges — rather than an undirected graph with side-flags — means the shortest-path search physically cannot traverse a one-way road backwards. Travel time, not raw length, is the edge weight, because a short segment on a 25 mph residential street is slower than a long segment on a highway and the evacuation cares about time.

from __future__ import annotations

import logging
from dataclasses import dataclass
from typing import Iterable

import networkx as nx
from shapely.geometry import LineString

logger = logging.getLogger("incidentgis.evac.graph")


class RoadGraphError(RuntimeError):
    """Raised when the input segments cannot produce a routable graph."""


@dataclass(frozen=True)
class RoadSegment:
    seg_id: str
    u: int                    # from-node id
    v: int                    # to-node id
    geometry: LineString      # projected CRS (metres), not degrees
    length_m: float
    speed_mps: float          # free-flow speed
    oneway: bool
    contraflow_reversible: bool = False


def build_directed_graph(
    segments: Iterable[RoadSegment],
    contraflow_active: bool = False,
    contraflow_factor: float = 0.6,
) -> nx.DiGraph:
    """Build a directed, travel-time-weighted road graph.

    Two-way segments add both directions. One-way segments add only u->v,
    unless contraflow is active and the segment is flagged reversible, in
    which case the reverse lane opens at ``contraflow_factor`` of free-flow.
    """
    graph = nx.DiGraph()
    added = 0
    for seg in segments:
        try:
            if seg.length_m <= 0 or seg.speed_mps <= 0:
                raise RoadGraphError(
                    f"segment {seg.seg_id} has non-positive length or speed"
                )
            forward_time = seg.length_m / seg.speed_mps
            graph.add_edge(
                seg.u, seg.v, seg_id=seg.seg_id, length_m=seg.length_m,
                travel_time_s=forward_time, geometry=seg.geometry, direction="forward",
            )
            if not seg.oneway:
                graph.add_edge(
                    seg.v, seg.u, seg_id=seg.seg_id, length_m=seg.length_m,
                    travel_time_s=forward_time, geometry=seg.geometry, direction="reverse",
                )
            elif contraflow_active and seg.contraflow_reversible:
                # Contraflow opens the reverse lane at a governed speed only.
                graph.add_edge(
                    seg.v, seg.u, seg_id=seg.seg_id, length_m=seg.length_m,
                    travel_time_s=seg.length_m / (seg.speed_mps * contraflow_factor),
                    geometry=seg.geometry, direction="contraflow",
                )
            added += 1
        except RoadGraphError as exc:
            # A single malformed segment must not abort the whole build.
            logger.error("skip_segment", extra={"seg_id": seg.seg_id, "err": str(exc)})
            continue
    if graph.number_of_edges() == 0:
        raise RoadGraphError("no routable edges produced from input segments")
    logger.info(
        "graph_built",
        extra={"segments": added, "nodes": graph.number_of_nodes(),
               "edges": graph.number_of_edges()},
    )
    return graph

Build the graph deterministically: iterate segments in a stable, sorted order so that any ties in the later shortest-path search are broken identically on every run. Reproducibility is not a nicety here — a route that changes between two identical runs cannot be defended in an after-action review.

Step 2 — Apply live hazard closures as a spatial overlay

When a new fire or flood perimeter arrives, do not rebuild the graph. Index the edge geometries once in a shapely STRtree and query it with the hazard polygon; only the handful of edges that actually intersect the perimeter need to change. The choice between removing an edge and penalizing it matters operationally: removal is correct for a flooded road that is physically impassable, while a large penalty is correct for a smoke-degraded corridor that should be used only if every clean route is gone. The live perimeters themselves come off the message backbone described in WebSocket & MQTT for Live Incident Feeds; this function assumes the polygon is already validated and in the routing CRS.

from __future__ import annotations

import logging

import networkx as nx
from shapely.geometry import LineString, Polygon
from shapely.strtree import STRtree

logger = logging.getLogger("incidentgis.evac.hazard")


class ClosureError(RuntimeError):
    """Raised when a hazard polygon cannot be resolved against the graph."""


def apply_hazard_closures(
    graph: nx.DiGraph,
    hazard: Polygon,
    buffer_m: float = 0.0,
    penalize_only: bool = False,
    penalty_factor: float = 1000.0,
) -> int:
    """Close or penalize every edge whose geometry intersects a hazard polygon.

    Returns the count of affected directed edges. When ``penalize_only`` is
    True the edges remain with an inflated travel time (last-resort routing);
    otherwise they are removed outright.
    """
    if hazard.is_empty or not hazard.is_valid:
        raise ClosureError("hazard polygon is empty or invalid")
    zone = hazard.buffer(buffer_m) if buffer_m else hazard

    edges = list(graph.edges(data=True))
    geoms: list[LineString] = [data["geometry"] for _, _, data in edges]
    tree = STRtree(geoms)  # bounding-box index for sub-linear candidate lookup

    affected: list[tuple[int, int]] = []
    # shapely >= 2.0: query() returns integer positions into ``geoms``.
    for pos in tree.query(zone):
        u, v, data = edges[int(pos)]
        # The bbox query is a coarse filter; confirm true intersection.
        if data["geometry"].intersects(zone):
            affected.append((u, v))

    for u, v in affected:
        if penalize_only:
            graph[u][v]["travel_time_s"] *= penalty_factor
            graph[u][v]["closed"] = "penalized"
        else:
            graph.remove_edge(u, v)

    logger.warning(
        "hazard_closures_applied",
        extra={"affected": len(affected),
               "mode": "penalize" if penalize_only else "remove",
               "buffer_m": buffer_m},
    )
    return len(affected)

Apply closures against a working copy of the graph (graph.copy()) if the base network must survive the incident untouched; removal is otherwise destructive and the base topology is gone. For a multi-day incident where perimeters advance and recede, a penalize-and-restore model on a persistent graph is usually cleaner than repeated copies.

Step 3 — Compute a capacity-aware route to the nearest safe destination

With the graph hazard-adjusted, routing is a single shortest-path search — but “shortest” must mean least time under load, not least distance. A road at capacity moves slowly, so the edge weight follows a Bureau of Public Roads (BPR) style congestion curve that inflates travel time as the volume-to-capacity ratio climbs. Rather than routing to each shelter separately and comparing, attach a zero-cost virtual sink to every open destination and run one search; Dijkstra then selects the fastest-reachable shelter itself. Feed it only destinations that are actually accepting evacuees, per the shelter capacity contract.

from __future__ import annotations

import logging
from typing import Callable, Hashable, Optional

import networkx as nx

logger = logging.getLogger("incidentgis.evac.route")


class NoRouteError(RuntimeError):
    """Raised when no open path exists from the origin to any safe node."""


def capacity_weight(
    congestion: dict[str, float],
    capacity_exponent: float = 4.0,
) -> Callable[[Hashable, Hashable, dict], float]:
    """Return a BPR-style edge-weight function.

    weight = free_flow_time * (1 + volume_capacity_ratio ** exponent)

    Congestion ratios arrive keyed by segment id from the live feed; a missing
    key is treated as free-flowing (0.0), never as an error.
    """
    def _weight(u: Hashable, v: Hashable, data: dict) -> float:
        base = float(data["travel_time_s"])
        ratio = max(0.0, float(congestion.get(data.get("seg_id", ""), 0.0)))
        return base * (1.0 + ratio ** capacity_exponent)
    return _weight


def route_to_safety(
    graph: nx.DiGraph,
    origin: Hashable,
    safe_nodes: set[Hashable],
    congestion: Optional[dict[str, float]] = None,
) -> list[Hashable]:
    """Least-time path from ``origin`` to the nearest open safe node.

    A virtual super-sink is attached to every safe node so a single Dijkstra
    call chooses the best destination. The sink is always removed afterwards
    so the graph is left pristine for the next query.
    """
    if origin not in graph:
        raise NoRouteError(f"origin node {origin!r} absent from graph")
    reachable_sinks = sorted(n for n in safe_nodes if n in graph)
    if not reachable_sinks:
        raise NoRouteError("no safe destination present in graph")

    weight_fn = capacity_weight(congestion or {})
    super_sink = "__SAFE_SINK__"
    graph.add_node(super_sink)
    for n in reachable_sinks:  # sorted -> deterministic tie-breaking
        graph.add_edge(n, super_sink, travel_time_s=0.0, seg_id="__sink__")
    try:
        path = nx.shortest_path(graph, origin, super_sink, weight=weight_fn)
    except nx.NetworkXNoPath as exc:
        raise NoRouteError(
            f"origin {origin!r} is isolated from every safe node"
        ) from exc
    finally:
        graph.remove_node(super_sink)  # never leave the scaffold behind

    route = path[:-1]  # drop the virtual sink from the returned path
    logger.info(
        "route_computed",
        extra={"origin": origin, "dest": route[-1], "hops": len(route)},
    )
    return route

The finally block is load-bearing: if the sink node were left attached after an exception, the next route query would traverse phantom zero-cost edges and return nonsense. Removing it unconditionally keeps the resident graph reusable across thousands of route requests during a surge.

Step 4 — Derive evacuation-zone isochrones

Incident command needs reachability, not just routes: which neighbourhoods can clear within fifteen minutes on the roads still open? An isochrone answers that. Run a single-source shortest-path search truncated at the time budget, then buffer the reachable nodes into a coverage surface. Reverse the graph before the search and the same code answers the mirror question — who can still reach this shelter — which drives demand estimates for staffing and supplies.

from __future__ import annotations

import logging
from typing import Hashable

import networkx as nx
from shapely.geometry import MultiPoint, Point, Polygon
from shapely.ops import unary_union

logger = logging.getLogger("incidentgis.evac.isochrone")


class IsochroneError(RuntimeError):
    """Raised when a reachability surface cannot be computed."""


def evacuation_isochrone(
    graph: nx.DiGraph,
    source: Hashable,
    cutoff_s: float,
    node_points: dict[Hashable, Point],
    buffer_m: float = 60.0,
) -> Polygon:
    """Reachability polygon: everywhere reachable from ``source`` within cutoff.

    Uses single-source Dijkstra truncated at the time budget, then buffers the
    reachable node points into a coverage surface. ``node_points`` maps node id
    to a projected Point so the buffer distance is a true metre radius.
    """
    if source not in graph:
        raise IsochroneError(f"source node {source!r} absent from graph")
    if cutoff_s <= 0:
        raise IsochroneError("cutoff must be a positive number of seconds")
    try:
        lengths = nx.single_source_dijkstra_path_length(
            graph, source, cutoff=cutoff_s, weight="travel_time_s"
        )
    except nx.NetworkXError as exc:
        raise IsochroneError("reachability search failed") from exc

    pts = [node_points[n] for n in lengths if n in node_points]
    if len(pts) < 3:
        raise IsochroneError("insufficient reachable nodes for a surface")

    surface = unary_union([p.buffer(buffer_m) for p in pts])
    logger.info(
        "isochrone_built",
        extra={"source": source, "reachable": len(pts), "cutoff_s": cutoff_s},
    )
    if isinstance(surface, Polygon):
        return surface
    # Disconnected buffers -> return the convex hull as a coverage envelope.
    return MultiPoint(pts).convex_hull

For production isochrones over a dense network, buffer-and-union of node points is a fast approximation; where a smoother boundary matters, sample points along each reachable edge or feed the reachable subgraph to a concave-hull routine. The approximation is deliberate — during an active evacuation, a surface computed in tens of milliseconds and refreshed every perimeter update beats a perfect boundary that arrives after the fire has moved.

Configuration Reference

Tune these per deployment; a dense metropolitan graph and a rural county network will not share thresholds.

Verification & Smoke Test

Run these assertions on a staging node before promoting any change to the routing engine. They confirm that a hazard closure forces a detour, that routing is deterministic, and that the graph survives a route query intact.

import logging

from shapely.geometry import LineString, Polygon

logger = logging.getLogger("incidentgis.evac.smoke")


def smoke_test() -> None:
    # A unit square: node 0 top-left, 1 top-right, 2 bottom-left, 3 bottom-right.
    segs = [
        RoadSegment("s01", 0, 1, LineString([(0, 100), (100, 100)]), 100.0, 13.9, False),
        RoadSegment("s13", 1, 3, LineString([(100, 100), (100, 0)]), 100.0, 13.9, False),
        RoadSegment("s02", 0, 2, LineString([(0, 100), (0, 0)]), 100.0, 13.9, False),
        RoadSegment("s23", 2, 3, LineString([(0, 0), (100, 0)]), 100.0, 13.9, False),
    ]
    graph = build_directed_graph(segs)

    # 1. A hazard over the northern edge forces the southern detour 0 -> 2 -> 3.
    hazard = Polygon([(40, 92), (110, 92), (110, 108), (40, 108)])
    closed = apply_hazard_closures(graph, hazard)
    assert closed >= 1, "hazard must close the northern edge"
    route = route_to_safety(graph, origin=0, safe_nodes={3})
    assert 2 in route and 1 not in route, "route must detour south, not through the closure"

    # 2. Determinism: identical inputs yield an identical route.
    g2 = build_directed_graph(segs)
    apply_hazard_closures(g2, hazard)
    assert route == route_to_safety(g2, origin=0, safe_nodes={3}), \
        "routing must be deterministic for a fixed graph and hazard"

    # 3. The virtual sink is always cleaned up.
    assert "__SAFE_SINK__" not in graph, "super-sink must not survive a route query"

    logger.info("smoke test passed")


smoke_test()

A one-line CLI check confirms the stack is wired before the engine is deployed to a field node:

python -c "import networkx, shapely, pyproj; print('routing stack ok')"

Integration With Adjacent Workflows

Evacuation routing sits downstream of nearly every other sync workflow. Its hazard closures are driven by perimeters delivered over the WebSocket & MQTT for Live Incident Feeds backbone, so a dropped or replayed perimeter message translates directly into a wrong route — the delivery guarantees there are part of this engine’s correctness envelope. Every length, buffer, and isochrone radius is metres, which is only coherent once the network and hazards share the projected reference frame that Coordinate Reference Systems for Disaster Zones establishes; route in degrees and the BPR weights and isochrones silently corrupt. Destinations come from the shelter capacity and resource tracking schemas, so the router only ever offers open, under-capacity facilities. Because National Incident Management System (NIMS) doctrine and Federal Emergency Management Agency (FEMA) after-action review both require that operational decisions be reconstructable, the engine must emit a route record — origin, chosen destination, the perimeter version in force, the graph revision, and the ordered edge list — into the incident audit log, and it can publish reachability surfaces to an Open Geospatial Consortium (OGC) API – Features service for the common operating picture. When a perimeter advances mid-evacuation and invalidates a route already issued, the incremental reroute path is the subject of rerouting around dynamically closed roads during flooding.

Troubleshooting

Symptom: the route runs straight through the hazard polygon. Almost always the network is in geographic degrees, so a buffer_m in metres is meaningless and the closure barely touches the edges, or the graph and the hazard are in different CRSes and never truly intersect. Reproject both into the shared projected CRS before calling apply_hazard_closures, and confirm hazard.intersects returns True for a known-closed edge in a unit test.

Symptom: travel times and isochrones are absurdly small or large. Edge weights were computed from a length_m that is actually in degrees, so a 100 m block reads as 0.001 “metres”. Verify the segment geometries are projected and that length_m matches geometry.length to within rounding before trusting any route.

Symptom: evacuees are directed the wrong way up a one-way street. A segment reached build_directed_graph with oneway=False because the source attribute was missing or parsed as a string "no". Validate the oneway flag as a real boolean at ingest, and confirm the reverse edge only ever appears for two-way roads or an active contraflow-reversible segment.

Symptom: NoRouteError is raised even though roads out of the zone clearly exist. Either the hazard buffer is large enough to close the last remaining corridor, or the origin snapped to a node on a disconnected graph component. Fall back to penalize_only=True so the corridor stays usable at high cost, and validate that the origin is in the graph’s largest strongly-connected component before routing.

Symptom: the same request returns a different route on successive runs. Segments or safe nodes are being fed from an unordered set or dict, so tie-breaking in Dijkstra varies. Sort the segment iterable before build_directed_graph and pass safe_nodes through the sorted handling already in route_to_safety; identical inputs must always produce byte-identical routes for the audit trail to hold.

Frequently Asked Questions

Why route on a directed road graph instead of straight-line distance to the nearest shelter? Straight-line distance ignores one-way streets, rivers, closed bridges, and the hazard perimeter itself, so the nearest shelter by Euclidean distance is frequently unreachable or on the wrong side of the fire. A directed graph encodes travel time along real, currently-open segments, respects one-way and contraflow rules, and lets a single least-time search pick the shelter that is genuinely fastest to reach. It is also the only representation that can be audited edge by edge after the incident.

How should live hazard closures be applied without rebuilding the whole graph? Keep the base road graph in memory and apply closures as a fast overlay: index every edge geometry in an STRtree, query it with the incoming hazard polygon, and remove or penalize only the intersecting edges. Rebuilding a metropolitan graph on every perimeter update costs seconds you do not have during a surge; an indexed overlay closes the affected edges in milliseconds and is trivially reversible when the hazard recedes.

Which coordinate reference system should evacuation routing run in? A projected, metre-based system such as the local UTM zone, never geographic degrees. Edge weights, hazard buffers, capacity math, and isochrone radii are all expressed in metres, and those values are meaningless in EPSG:4326 where a degree of longitude shrinks toward the poles. Reproject the network and every hazard polygon into the same projected CRS before routing, then transform results back to WGS 84 only for display.

Related

• Rerouting Around Dynamically Closed Roads During Flooding
• WebSocket & MQTT for Live Incident Feeds
• Coordinate Reference Systems for Disaster Zones
• Shelter Capacity & Resource Tracking Schemas

Up: Incident Mapping & Multi-Agency Sync Workflows