Spatial Epi Public Health GIS

Modeling Hospital Bed Capacity with Network Analysis

This guide solves a single operational problem: assigning surge-period patient origins to hospital beds using impedance-weighted road routing that reacts to live occupancy, rather than the static Euclidean catchments that silently misallocate during a saturation event. It is part of the Facility Capacity Allocation Models collection within Healthcare Access & Network Analysis Automation.

Problem Context & Constraints

The naive approach — assign each patient origin to the nearest facility by straight-line distance, then check whether that facility has a free bed — fails for three reasons that compound during a disease surge.

First, Euclidean distance ignores the road network. Two origins equidistant from a hospital as the crow flies can have travel times that differ by a factor of three once rivers, limited-access corridors, and one-way restrictions are accounted for. During a surge the cheapest geographic assignment is frequently not the cheapest reachable one.

Second, nearest-facility assignment is occupancy-blind. A static catchment keeps routing arrivals to a hospital that is already past its diversion threshold, manufacturing queueing where the model reports adequate access. Capacity has to enter the cost function itself, not a post-hoc filter — once the nearest facility is full, every downstream reassignment must re-rank the reachable set, which a distance lookup cannot do.

Third, the network is non-stationary. Edge speeds shift with time-of-day traffic and incident closures, so a travel-time matrix computed once at the start of a shift drifts out of calibration within an hour. The matrix must be recomputed on a cadence, and every recomputation must be reproducible from a recorded network snapshot for later defensibility.

The technique below replaces the distance lookup with a routing graph whose terminal-node impedance is a function of live bed occupancy. When a facility crosses its clinical diversion threshold, its incoming edges are penalized by a configurable multiplier, so the shortest-path solver naturally spills demand onto the next-reachable hospital with headroom. This is the mechanism that lets the broader Facility Capacity Allocation Models operate against live state instead of a nightly snapshot. The same isochrone fabric produced by Drive-Time Isochrone Generation feeds the travel-time edges this routine consumes.

Impedance-Weighted Bed-Capacity Routing Three inputs — de-identified patient origins, the projected road graph, and a live bed-occupancy feed — converge on an impedance-weighted shortest-path router. When a facility reaches or exceeds 85 percent occupancy its incoming edges are penalized by the impedance multiplier, so demand spills from the saturated hospital to the next reachable hospital that still has headroom. INPUTS Patient origins de-identified points Road graph travel_time edges Live occupancy current / total beds Impedance- weighted router shortest path on travel_time occupancy ≥ 85% → edges × multiplier REACHABLE FACILITIES Hospital A 92% — saturated, penalized spill Hospital B 61% — has headroom RUN_MANIFEST snapshot hash · CRS · diversion params Capacity enters the cost function, so saturated facilities divert demand to the next reachable hospital

Prerequisites

Pin the GIS stack so impedance arithmetic and nearest-node snapping are reproducible across runs:

# requirements.txt — pinned for reproducible surge runs
osmnx==1.9.3
networkx==3.3
geopandas==1.0.1
shapely==2.0.5
pyproj==3.6.1
pandas==2.2.2
scipy==1.14.0

Input data and CRS state required before the routine runs:

  • Road graph — an OSMnx GraphML export with travel_time (seconds) and length (metres) on every edge. If travel_time is absent, derive it from length and maxspeed with ox.add_edge_speeds / ox.add_edge_travel_times before serializing.
  • Patient origins — a GeoDataFrame of point geometries that has already been de-identified: aggregate to census block-group centroids, apply spatial k-anonymity, or inject calibrated Laplace noise. Raw protected-health-information coordinates must never reach this stage. Aggregation policy is governed by the Compliance Mapping Frameworks.
  • Capacity table — a DataFrame keyed by facility node_id with current_beds and total_beds, refreshed at the same cadence as the routing cycle.
  • A single metric CRS. All inputs must share one projected, metre-based coordinate reference system — a UTM zone (EPSG:326xx) or a state-plane projection. Mixing geographic and projected layers produces cumulative distance error across thousands of origin-facility pairs. CRS selection rationale is covered in Coordinate Reference Systems for Public Health.

Step-by-Step Solution

The routine projects every layer to one CRS, keeps the largest weakly connected component so unreachable subgraphs cannot raise no-path errors, penalizes the incoming edges of saturated facilities, then routes each origin to its least-impeded reachable hospital. Every consequential decision is logged so the run can be reconstructed from the audit trail.

# Hospital bed capacity routing with dynamic impedance.
# osmnx==1.9.3, networkx==3.3, geopandas==1.0.1, pandas==2.2.2, pyproj==3.6.1
import hashlib
import json
import logging
from datetime import datetime, timezone

import geopandas as gpd
import networkx as nx
import osmnx as ox
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")


def build_dynamic_capacity_routing(
    network_path: str,
    capacity_df: pd.DataFrame,
    origins_gdf: gpd.GeoDataFrame,
    saturation_threshold: float = 0.85,
    impedance_multiplier: float = 3.0,
    target_crs: str = "EPSG:32148",  # WA state plane (metres); set per study area
) -> pd.DataFrame:
    """Assign de-identified origins to the least-impeded reachable facility.

    capacity_df must contain: node_id, current_beds, total_beds.
    Deterministic: origins and facilities are processed in stable sorted order.
    """
    run_id = datetime.now(timezone.utc).isoformat(timespec="seconds")

    # 1. CRS enforcement — never route across mixed CRS.
    if str(origins_gdf.crs) != target_crs:
        origins_gdf = origins_gdf.to_crs(target_crs)
        logging.info("Reprojected origins to %s", target_crs)

    # 2. Load and project the road graph to match the origins.
    G = ox.load_graphml(network_path)
    G = ox.project_graph(G, to_crs=target_crs)

    # Record an immutable snapshot hash for regulatory traceability.
    snapshot_hash = hashlib.sha256(
        f"{network_path}:{G.number_of_nodes()}:{G.number_of_edges()}".encode()
    ).hexdigest()[:16]

    # 3. Keep the largest weakly connected component to avoid no-path errors
    #    from isolated import artefacts.
    largest_cc = max(nx.weakly_connected_components(G), key=len)
    G = G.subgraph(largest_cc).copy()
    logging.info("Retained %d nodes in largest WCC (snapshot %s)",
                 G.number_of_nodes(), snapshot_hash)

    # 4. Apply diversion impedance to the incoming edges of saturated facilities.
    #    Sort for deterministic, reproducible application order.
    diverted = []
    for _, row in capacity_df.sort_values("node_id").iterrows():
        node_id = row["node_id"]
        if node_id not in G.nodes:
            logging.warning("Facility node %s not in graph; skipping.", node_id)
            continue
        occupancy = row["current_beds"] / row["total_beds"]
        if occupancy >= saturation_threshold:
            for _u, _v, data in G.in_edges(node_id, data=True):
                base = data.get("travel_time", data.get("length", 1.0))
                data["travel_time"] = base * impedance_multiplier
            diverted.append(node_id)
            logging.info("Diversion x%.1f applied to node %s (occupancy %.1f%%)",
                         impedance_multiplier, node_id, occupancy * 100)

    facility_nodes = [n for n in capacity_df["node_id"].tolist() if n in G.nodes]

    # 5. Route each origin to its least-impeded reachable facility.
    results = []
    for idx in origins_gdf.sort_index().index:
        geom = origins_gdf.loc[idx, "geometry"]
        try:
            src = ox.distance.nearest_nodes(G, geom.x, geom.y)
            best_time, best_facility = None, None
            for fnode in facility_nodes:
                try:
                    t = nx.shortest_path_length(G, src, fnode, weight="travel_time")
                except nx.NetworkXNoPath:
                    continue
                if best_time is None or t < best_time:
                    best_time, best_facility = t, fnode
            results.append({
                "origin_id": idx,
                "assigned_facility": best_facility,
                "travel_time_min": round(best_time / 60, 2) if best_time is not None else None,
                "diverted_from_saturation": best_facility not in diverted if best_facility else None,
            })
        except Exception as exc:  # isolate per-origin failures
            logging.warning("Routing failed for origin %s: %s", idx, exc)
            results.append({"origin_id": idx, "assigned_facility": None,
                            "travel_time_min": None, "diverted_from_saturation": None})

    # 6. Emit a machine-readable run manifest for the audit trail.
    manifest = {
        "run_id": run_id,
        "snapshot_hash": snapshot_hash,
        "target_crs": target_crs,
        "saturation_threshold": saturation_threshold,
        "impedance_multiplier": impedance_multiplier,
        "diverted_facilities": diverted,
        "origins_routed": len(results),
    }
    logging.info("RUN_MANIFEST %s", json.dumps(manifest, sort_keys=True))

    return pd.DataFrame(results)

Persist the returned frame to partitioned GeoParquet keyed by run_id, and persist the logged RUN_MANIFEST alongside it so each allocation cycle is reconstructable from a single key. The batch wrapper, retry, and isolation patterns this routine assumes are documented in Batch Routing & Error Handling.

Validation & Edge Cases

Three failure modes account for nearly every silent miscount in surge routing. Each has a cheap diagnostic.

1. Saturated-facility starvation. If every facility within reach of a region is past threshold, the multiplier inflates all candidate paths equally and the solver still picks one — there is nowhere to spill demand. Detect it by checking whether the assigned travel time is implausibly large:

overflow = result[result["travel_time_min"] > 90]
if not overflow.empty:
    logging.error("No-headroom region: %d origins routed > 90 min", len(overflow))

A nonzero count means real-world diversion or a temporary facility is required; it is not a code bug, but it must surface rather than hide inside an inflated cost.

2. Off-network origins. nearest_nodes always returns a node, even for an origin snapped kilometres away — for example a centroid that landed in water after de-identification noise. Guard the snap distance explicitly:

# expected log when an origin snaps too far from any road
# WARNING | Origin 4471 snapped 5128 m to node 32004187 (tolerance 1500 m)
node, dist = ox.distance.nearest_nodes(G, geom.x, geom.y, return_dist=True)
if dist > 1500:
    logging.warning("Origin %s snapped %.0f m to node %s (tolerance 1500 m)", idx, dist, node)

3. Missing travel_time edges. A graph exported without travel times falls back to raw length (metres), so impedance is computed in mixed units and travel_time_min is meaningless. Assert the attribute exists before routing:

missing = [1 for _u, _v, d in G.edges(data=True) if "travel_time" not in d]
assert not missing, f"{len(missing)} edges missing travel_time; run ox.add_edge_travel_times first"

Before deployment, audit assignments for spatial bias: compute a Two-Step Floating Catchment Area or weighted travel-time accessibility score across demographic strata and confirm the diversion logic does not concentrate longer travel times in any single group. That equity audit is the subject of Spatial Equity Index Calculation.

Compliance Notes

For an allocation run to be regulatorily defensible, the following must be recorded and retained with the output:

  • Network snapshot hash (snapshot_hash) and the GraphML source path — proves which road state produced the assignment.
  • CRS identifier (target_crs) — proves distances were computed in a metric projection, not on mixed or geographic coordinates.
  • Diversion parameters (saturation_threshold, impedance_multiplier) and the explicit list of diverted_facilities — proves why demand was redistributed.
  • Run timestamp (run_id, UTC) and origin count — anchors the cycle in the surge-response timeline.
  • De-identification provenance — the aggregation or geomasking method applied to origins upstream, so no raw protected-health-information coordinate is ever implied in the audit chain.

Emit these as the single RUN_MANIFEST JSON line shown above and store it with the GeoParquet output. Deterministic execution — sorting facilities and origins by stable identifier — guarantees that re-running the same inputs reproduces the same assignment, which is the baseline a regulator expects when validating a constrained-capacity decision.