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 Drive-Time Isochrone Generation, within Healthcare Access & Network Analysis Automation.

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:

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.
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.
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 projection discipline enforced when you align WGS84 and UTM for county health data before any distance calculation, drawn from the broader Coordinate Reference Systems for Public Health method set.

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 — for example, the pedestrian catchments produced here become a supply-side service area inside a Two-Step Floating Catchment Area index.

FAQ

Why a convex hull instead of buffering the reachable nodes?

A convex hull wraps the reachable frontier in a single deterministic polygon with no tunable radius, which is exactly what a batch county-level run needs: it is fast, parameter-free, and conservative. Buffering individual nodes reintroduces a distance parameter and can fuse or over-extend the area across gaps the network does not actually connect. When you need the polygon to hug the network rather than over-claim, swap the hull for alpha shapes or network-constrained buffering and accept the higher per-clinic runtime.

Why penalize speed instead of just shrinking the time threshold?

The 900-second threshold is the policy definition of a “15-minute walk” and must stay fixed for the result to be comparable and defensible across jurisdictions. Terrain difficulty is a property of the edges, not of the policy, so it belongs in the per-edge travel_time_sec. Applying the 0.85 surface and 0.75 incline penalties at the edge level lets a paved corridor and an unpaved grade-4 track both be measured against the same 900-second budget while reaching different distances.

A clinic returns an empty isochrone but sits on a visible footpath — what happened?

Almost always topological fragmentation: the clinic snapped to a small connected component that is severed from the corridor you can see on the map, usually because of a missing or mistagged highway segment in the OpenStreetMap extract. Confirm with the per-component node-count diagnostic, check that the clinic node is inside the largest component, and either extend the custom_filter or stitch the components with a short snapping tolerance before re-routing.

Do I have to project the graph, or can I route in WGS84 degrees?

You must project. Summing edge length in WGS84 degrees mixes units and distorts distance with latitude, which silently invalidates the 900-second threshold. ox.project_graph resolves a local metric CRS (UTM or state plane) so every length is in meters and the time math is correct; record the resolved EPSG code in the run log for reproducibility.

Drive-Time Isochrone Generation — parent method set for vehicular and multimodal service areas.
Calculating the Two-Step Floating Catchment Area Index — the access metric that consumes these pedestrian catchments.
Spatial Equity Index Calculation — broader disparity scoring built on these isochrones.
How to Align WGS84 and UTM for County Health Data — the projection discipline that keeps the 15-minute threshold valid.