Generating 15-Minute Walk Isochrones for Rural Clinics

This guide solves one narrow but recurring problem: producing defensible 15-minute pedestrian service areas around rural clinics when the underlying walk network is sparse, fragmented, and unevenly tagged. It is part of the Drive-Time Isochrone Generation method set within the broader Healthcare Access & Network Analysis Automation pipeline.

Problem Context & Constraints

The naive approach — pull a network_type="walk" graph, assume a uniform 5.0 km/h pace, and buffer reachable nodes — works acceptably in dense urban grids and fails predictably in rural jurisdictions. Three constraints break the default:

1. Topological fragmentation. OpenStreetMap (OSM) extracts in rural areas often omit explicit footway tags, routing pedestrians onto highway=track, unclassified, or service corridors. The walk graph splinters into many weakly connected components, so a clinic snapped to a 4-node stub returns a near-empty isochrone while the real reachable area sits in a neighbouring component.
2. Heterogeneous edge friction. A single 5.0 km/h baseline overestimates the population reachable on foot. Unpaved surfaces, drainage ditches, and elevation gradients depress actual walking velocity, and a fixed pace silently inflates accessible-population counts in exactly the terrain where access is hardest.
3. Geographic-coordinate distance drift. Summing edge length in WGS84 degrees produces nonsense thresholds. Distance and time math is only valid after projecting to a meter-based CRS — the same discipline enforced for every distance calculation across Coordinate Reference Systems for Public Health.

The pipeline below addresses all three with explicit graph validation, terrain-adjusted friction weighting, deterministic node snapping, and audit-ready serialization for county-level clinic inventories.

Prerequisites

Pin the GIS stack so isochrone geometries are reproducible across agency machines:

# requirements (pinned):
#   osmnx==1.9.3
#   networkx==3.2.1
#   geopandas==0.14.4
#   shapely==2.0.4
#   pyproj==3.6.1
#   pyarrow==15.0.2   # GeoParquet writer

Input requirements:

• boundary_gdf — a single-row (or dissolvable) GeoDataFrame of the study jurisdiction, with a defined CRS. Used both to bound the OSM query and as the canonical output CRS.
• clinics_gdf — point geometries with a stable clinic_id column. A stable identifier is mandatory for deterministic, auditable output ordering.
• A meter-based projected CRS (local UTM zone or state plane) is selected automatically by osmnx; do not run any of the routing steps in geographic coordinates.

Step-by-Step Solution

1. Network acquisition and metric projection

Querying via osmnx requires explicit filter expansion to capture viable rural pedestrian corridors while excluding high-speed vehicular routes. Projection to a metric CRS happens immediately, before any edge-length is consumed.

import osmnx as ox
import networkx as nx
import geopandas as gpd
from shapely.geometry import Point
import logging

# Deterministic caching; route OSMnx output through Python logging
ox.settings.use_cache = True
ox.settings.log_console = False

def acquire_rural_walk_graph(boundary_gdf: gpd.GeoDataFrame) -> nx.MultiDiGraph:
    boundary_geom = (
        boundary_gdf.union_all() if hasattr(boundary_gdf, "union_all")
        else boundary_gdf.unary_union
    )
    G = ox.graph_from_polygon(
        boundary_geom,
        network_type="walk",
        custom_filter='["highway"~"footway|path|track|unclassified|residential|service|pedestrian"]',
    )
    # Immediate projection to a local metric CRS is mandatory for distance/time math.
    # Operating in WGS84 accumulates distortion during edge-length summation and
    # invalidates the 15-minute threshold.
    G_proj = ox.project_graph(G)
    return G_proj

Refer to the OSMnx documentation for the projection and coordinate-transformation contract.

2. Terrain-adjusted friction weighting

Friction is applied at the edge level, before any routing computation, producing a travel_time_sec attribute that networkx consumes directly. Penalties are conservative by design, biasing toward under-counting reachable population in degraded terrain.

def apply_rural_friction(G: nx.MultiDiGraph) -> nx.MultiDiGraph:
    for u, v, key, data in G.edges(data=True, keys=True):
        base_speed_kmh = 5.0

        surface = data.get("surface", "paved")
        if isinstance(surface, str):
            surface = surface.lower()

        track_type = data.get("tracktype", "")
        if isinstance(track_type, str):
            track_type = track_type.lower()

        incline = data.get("incline", 0.0)
        try:
            incline = float(incline)
        except (ValueError, TypeError):
            incline = 0.0

        # Surface degradation penalty
        if surface in ("unpaved", "gravel", "dirt", "earth") or "grade3" in track_type or "grade4" in track_type:
            base_speed_kmh *= 0.85

        # Elevation gradient penalty (Tobler's hiking function approximation)
        if abs(incline) > 0.08:
            base_speed_kmh *= 0.75

        length_m = data.get("length", 0.0)
        speed_ms = base_speed_kmh * 1000 / 3600
        data["travel_time_sec"] = length_m / speed_ms if speed_ms > 0 else float("inf")

    return G

3. Topology validation and clinic node snapping

Disconnected components and isolated clinic coordinates are the primary failure vectors in rural isochrone generation; unvalidated inputs trigger NetworkXNoPath during batch execution. Snapping is deterministic, distances are logged, and clinics outside the fallback radius are flagged rather than silently dropped.

def validate_and_snap_clinics(
    G: nx.MultiDiGraph,
    clinics_gdf: gpd.GeoDataFrame,
    max_snap_m: float = 500.0,
) -> dict:
    snapped_nodes = {}
    logging.info("Validating clinic network connectivity...")

    # Align clinics with the projected graph CRS so snap distances are in meters
    graph_crs = G.graph.get("crs", None)
    if graph_crs:
        clinics_gdf = clinics_gdf.to_crs(graph_crs)

    # Deterministic iteration order by stable clinic_id for audit reproducibility
    clinics_gdf = clinics_gdf.sort_values("clinic_id")

    for idx, row in clinics_gdf.iterrows():
        clinic_id = row.get("clinic_id", idx)
        clinic_x, clinic_y = row.geometry.x, row.geometry.y

        nearest_node = ox.distance.nearest_nodes(G, clinic_x, clinic_y, return_dist=False)
        node_data = G.nodes[nearest_node]
        dist_to_node = ((clinic_x - node_data["x"]) ** 2 + (clinic_y - node_data["y"]) ** 2) ** 0.5

        if dist_to_node > max_snap_m:
            logging.warning(
                f"Clinic {clinic_id} exceeds {max_snap_m}m snap radius "
                f"({dist_to_node:.1f}m). Flagged for manual review."
            )
            continue

        snapped_nodes[clinic_id] = nearest_node

    components = list(nx.weakly_connected_components(G))
    largest_cc = max(components, key=len)
    logging.info(f"Graph: {len(components)} components. Largest CC: {len(largest_cc)} nodes.")

    return snapped_nodes

4. Isochrone computation and polygon construction

With validated nodes and calibrated weights, the 15-minute (900-second) service area is a single-source shortest-path frontier. Reachable nodes are wrapped in a convex hull as a conservative, batch-efficient approximation.

def generate_walk_isochrones(
    G: nx.MultiDiGraph,
    snapped_nodes: dict,
    threshold_sec: float = 900.0,
) -> gpd.GeoDataFrame:
    isochrones = []
    graph_crs = G.graph.get("crs", None)

    for clinic_id, source_node in snapped_nodes.items():
        reachable_lengths = nx.single_source_dijkstra_path_length(
            G, source_node, weight="travel_time_sec", cutoff=threshold_sec
        )
        reachable_nodes = set(reachable_lengths.keys())

        if not reachable_nodes:
            logging.warning(f"No reachable nodes for clinic {clinic_id} within {threshold_sec}s.")
            continue

        coords = [(G.nodes[n]["x"], G.nodes[n]["y"]) for n in reachable_nodes]
        if len(coords) < 3:
            logging.warning(f"Clinic {clinic_id} reached <3 nodes; cannot form a polygon.")
            continue

        poly = gpd.GeoSeries([Point(c) for c in coords]).union_all().convex_hull
        isochrones.append(
            {"clinic_id": clinic_id, "geometry": poly, "threshold_sec": threshold_sec}
        )

    if not isochrones:
        return gpd.GeoDataFrame(
            {"clinic_id": [], "geometry": [], "threshold_sec": []}, crs=graph_crs
        )

    return gpd.GeoDataFrame(isochrones, crs=graph_crs)

The convex hull suits county-scale health-equity assessments. For higher precision, swap the hull for alpha shapes or network-constrained buffering — both increase per-clinic runtime substantially in batch workflows.

5. Batch wrapper and audit logging

The wrapper enforces CRS consistency, captures routing failures without aborting the inventory, and serializes audit-ready GeoParquet.

def run_rural_isochrone_pipeline(
    boundary_gdf: gpd.GeoDataFrame,
    clinics_gdf: gpd.GeoDataFrame,
    output_path: str,
):
    logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")

    try:
        G = acquire_rural_walk_graph(boundary_gdf)
        G = apply_rural_friction(G)
        snapped = validate_and_snap_clinics(G, clinics_gdf)
        isochrones_gdf = generate_walk_isochrones(G, snapped)

        # Enforce output CRS alignment with the input boundary
        if isochrones_gdf.crs and boundary_gdf.crs and isochrones_gdf.crs != boundary_gdf.crs:
            isochrones_gdf = isochrones_gdf.to_crs(boundary_gdf.crs)

        # GeoParquet (geopandas >= 0.12 + pyarrow)
        isochrones_gdf.to_parquet(output_path, compression="snappy")
        logging.info(
            f"Pipeline complete. {len(isochrones_gdf)} isochrones exported to {output_path}"
        )

    except Exception as e:
        logging.error(f"Pipeline aborted: {e}", exc_info=True)
        raise

Validation & Edge Cases

Three failure modes account for nearly every bad rural isochrone. Each surfaces in the structured log, so QA can triage from the log file alone.

1. Component fragmentation — a clinic returns a tiny or empty isochrone. When the largest connected component holds a small fraction of total nodes, clinics snapped to isolated stubs reach almost nothing. Watch the component summary line:

2026-06-25 09:14:02 | INFO | Graph: 47 components. Largest CC: 311 nodes.
2026-06-25 09:14:03 | WARNING | No reachable nodes for clinic RC-0188 within 900.0s.

A high component count relative to node count signals a missing or mistagged corridor in OSM. Diagnose with a per-component node count and confirm the clinic sits in the largest component:

sizes = sorted((len(c) for c in nx.weakly_connected_components(G)), reverse=True)
logging.info(f"Top-5 component sizes: {sizes[:5]} of {G.number_of_nodes()} nodes")

If fragmentation is an OSM artifact rather than real impassability, extend the custom_filter or stitch components with a short snapping tolerance before routing.

2. Snap radius exceeded — clinic excluded. A clinic more than max_snap_m from any node is flagged, not dropped silently:

2026-06-25 09:14:02 | WARNING | Clinic RC-0204 exceeds
    500.0m snap radius (842.6m). Flagged for review.

Seeing many of these warnings together usually means the clinic coordinates and the graph are in mismatched units or CRS — verify clinics_gdf.crs and that the projected graph CRS reports meters.

3. Degenerate polygon — fewer than three reachable nodes. Clinics on dead-end spurs reach one or two nodes and cannot form a hull; these emit a warning and are skipped rather than crashing the batch. Persistent cases indicate the clinic should be re-snapped to the nearest through-corridor.

Compliance Notes

For regulatory defensibility, the following must be captured in the run log or output schema so any isochrone can be reconstructed and defended:

• Parameters: threshold_sec (900), max_snap_m, the friction penalty factors (0.85 surface, 0.75 incline), and the resolved projected CRS / EPSG code.
• Determinism: clinics are processed in stable clinic_id order, and ox.settings.use_cache = True pins the OSM extract for a given run — record the cache timestamp or the OSM extract date alongside outputs.
• Transparency over silent loss: every excluded clinic (snap radius, no reachable nodes, degenerate polygon) is logged with its clinic_id. Excluding clinics silently would bias downstream access metrics; the warnings are the audit trail.
• Output schema: GeoParquet columns clinic_id, geometry, threshold_sec, serialized in the input boundary CRS so the isochrones overlay cleanly on the source jurisdiction.

These outputs feed directly into downstream Spatial Equity Index Calculation modules for disparity mapping and resource-allocation modeling.

Related Topics

• Drive-Time Isochrone Generation — parent method set for vehicular and multimodal service areas.
• Spatial Equity Index Calculation — consumes these isochrones for access-disparity scoring.
• Coordinate Reference Systems for Public Health — projection discipline that keeps the 15-minute threshold valid.