Creating 15-Minute City Isochrones in Python

Generating a 15-minute city isochrone requires combining OpenStreetMap network extraction, time-weighted graph traversal, and spatial polygon generation — the precise workflow covered by generating isochrones with PySAL and GeoPandas within the broader Python Routing Engines & Isochrone Mapping discipline. This page focuses on the specific variant where the goal is a single, accurate walk, bike, or transit accessibility polygon centred on one coordinate — the exact output needed for urban planning reports, logistics facility siting, and compliance mapping, where the polygon boundary itself is the deliverable rather than a ranked list of routes.

Unlike raster flood-fill approaches, the graph-traversal method described here follows actual road topology. The resulting boundary respects one-way streets, access restrictions, and real path connectivity, making it suitable for regulatory submissions and service-area modelling.

When to Use This Approach

Use the graph-traversal method when:

• The boundary polygon is the deliverable. Urban planning submissions, accessibility audits, and site selection reports need a single authoritative polygon, not a distance matrix or a list of candidate stops.
• Routing fidelity matters more than computation speed. Euclidean buffers or Voronoi-based approximations overestimate reachability by 20–40% in fragmented street grids. Dijkstra on an OSM-derived graph captures dead ends, waterways, rail barriers, and one-way constraints that geometric methods miss.
• Multiple travel modes are required. The same Dijkstra kernel handles walk, bike, and drive graphs by switching the network_type argument to graph_from_point — no separate algorithm is needed.
• Output feeds downstream spatial joins. A GeoDataFrame result integrates directly with census block or parcel datasets without intermediate format conversion.

This technique is less suitable when you need isochrones for hundreds of origins simultaneously (use a distance-matrix engine such as deploying OSRM with Docker for local routing for that scale) or when real-time traffic data must be reflected in the boundary.

Implementation

The function below implements the six-stage pipeline with modern shapely 2.x and geopandas 0.14+ APIs. It projects to a metric CRS for accurate buffering before returning a WGS84 GeoDataFrame.

# requires: osmnx>=1.9, networkx>=3.2, geopandas>=0.14, shapely>=2.0
import osmnx as ox
import networkx as nx
import geopandas as gpd
import shapely
from shapely.geometry import Point


def generate_15min_isochrone(
    lat: float,
    lon: float,
    travel_mode: str = "walk",
    speed_kmh: float = 5.0,
    buffer_m: float = 50.0,
    cutoff_s: int = 900,
) -> gpd.GeoDataFrame:
    """
    Generate a single isochrone polygon for a given coordinate.

    Returns a GeoDataFrame with one row in EPSG:4326.
    Raises ValueError if the graph has fewer than 4 reachable nodes.
    """
    # ① Download and simplify the street network
    G = ox.graph_from_point(
        (lat, lon),
        dist=3000,           # metres — sufficient for 15-min walk/bike
        network_type=travel_mode,
        simplify=True,
    )

    # ② Convert edge lengths (metres) to travel time (seconds) — vectorised
    speed_ms = (speed_kmh * 1000) / 3600
    nx.set_edge_attributes(
        G,
        {
            (u, v, k): data["length"] / speed_ms
            for u, v, k, data in G.edges(keys=True, data=True)
            if "length" in data
        },
        name="travel_time",
    )

    # ③ Snap origin coordinate to the nearest valid node
    origin_node = ox.distance.nearest_nodes(G, lon, lat)

    # ④ Single-source Dijkstra with temporal cutoff
    reachable: dict = nx.single_source_dijkstra_path_length(
        G, origin_node, weight="travel_time", cutoff=cutoff_s
    )

    if len(reachable) < 4:
        raise ValueError(
            f"Only {len(reachable)} nodes reachable. "
            "Check coordinates, network_type, or increase dist."
        )

    # ⑤ Extract node geometries and project for metric buffering
    node_data = G.nodes(data=True)
    points_gdf = gpd.GeoDataFrame(
        geometry=[Point(node_data[n]["x"], node_data[n]["y"]) for n in reachable],
        crs="EPSG:4326",
    ).to_crs(epsg=3857)

    buffered = points_gdf.geometry.buffer(buffer_m)
    merged = shapely.union_all(buffered)

    # ⑥ Concave hull removes buffer artifacts; ratio 0.0 = convex, 1.0 = tightest
    isochrone_metric = shapely.concave_hull(merged, ratio=0.9)

    # Reproject result to WGS84 for standard GIS interoperability
    iso_gdf = gpd.GeoDataFrame(
        {"travel_time_max_s": [cutoff_s], "travel_mode": [travel_mode]},
        geometry=[isochrone_metric],
        crs="EPSG:3857",
    ).to_crs(epsg=4326)

    return iso_gdf


# --- Stacked isochrones (5 / 10 / 15 minutes) ---
# Run Dijkstra once, partition the result — avoids three graph traversals.
def generate_stacked_isochrones(
    lat: float,
    lon: float,
    thresholds_s: list[int] = [300, 600, 900],
    travel_mode: str = "walk",
    speed_kmh: float = 5.0,
    buffer_m: float = 50.0,
) -> gpd.GeoDataFrame:
    G = ox.graph_from_point((lat, lon), dist=3000, network_type=travel_mode, simplify=True)
    speed_ms = (speed_kmh * 1000) / 3600
    nx.set_edge_attributes(
        G,
        {(u, v, k): d["length"] / speed_ms for u, v, k, d in G.edges(keys=True, data=True) if "length" in d},
        name="travel_time",
    )
    origin = ox.distance.nearest_nodes(G, lon, lat)
    all_costs = nx.single_source_dijkstra_path_length(G, origin, weight="travel_time")

    node_data = G.nodes(data=True)
    bands = []
    for cutoff in thresholds_s:
        nodes_in_band = [n for n, t in all_costs.items() if t <= cutoff]
        pts = gpd.GeoDataFrame(
            geometry=[Point(node_data[n]["x"], node_data[n]["y"]) for n in nodes_in_band],
            crs="EPSG:4326",
        ).to_crs(epsg=3857)
        merged = shapely.union_all(pts.geometry.buffer(buffer_m))
        hull = shapely.concave_hull(merged, ratio=0.9)
        bands.append({"travel_time_max_s": cutoff, "geometry": hull})

    return gpd.GeoDataFrame(bands, crs="EPSG:3857").to_crs(epsg=4326)

The nx.set_edge_attributes call with a dictionary comprehension is meaningfully faster than a Python-level loop over G.edges(data=True) for graphs with tens of thousands of edges, because the attribute assignment is handled in a single C-level dict update rather than one Python call per edge.

Key Parameters and Tuning

For topography-sensitive projects, replace speed_kmh with a variable speed function derived from the OSM ele node attribute or a DEM raster lookup. An 8% gradient typically reduces effective walking speed by 25–30%, which shifts the isochrone boundary inward by 150–250 m for hilly urban terrain.

Integration Points

The GeoDataFrame produced by this function is immediately usable in several downstream workflows:

Spatial joins with census or parcel data. A gpd.sjoin against a population or land-use GeoDataFrame quantifies the population within each isochrone ring without any additional geometry construction.

Overlay with GTFS stop buffers. When modelling transit access, dissolve the OSM walk isochrone with transit catchment polygons derived from a GTFS feed. This is the foundation of multi-modal accessibility analysis as covered in implementing multi-modal transit layers.

Custom cost functions. Replacing the uniform speed_kmh impedance with a per-edge cost table (e.g. surface quality, elevation, or cargo weight) is the same operation as configuring edge weights for freight logistics — the only difference is the weight attribute name passed to single_source_dijkstra_path_length.

Routing engine handoff. Export the polygon as GeoJSON or GeoPackage and load it into OSRM or Valhalla as a permitted service area filter, or use it as the geographic extent for a turn-restriction graph subset. For OSRM-specific integration, see integrating custom traffic weights into OSRM.

Validation Checklist

After generating an isochrone, verify the following before publishing or committing the polygon to a downstream system:

1. Topology validity. assert iso_gdf.geometry.iloc[0].is_valid must pass. Invalid geometries (self-intersections, dangling edges) will silently corrupt spatial joins. Run shapely.make_valid(geom) to repair if needed.
2. Area reasonableness. In a dense urban grid, a 15-minute walk isochrone covers approximately 0.8–1.5 km². Project to a local equal-area CRS and check: iso_gdf.to_crs(epsg=6933).geometry.area.iloc[0] / 1e6 should fall in that range. Values above 2.5 km² usually indicate an overly large buffer_m or incorrect speed.
3. No clipping at download boundary. Verify iso_gdf.geometry.iloc[0].bounds does not coincide with the graph’s bounding box. If it does, increase dist and re-run.
4. Origin containment. The origin coordinate must lie within the polygon: assert iso_gdf.geometry.iloc[0].contains(Point(lon, lat)). If it does not, the origin node snap failed or the polygon is fragmented.
5. CRS correctness. assert iso_gdf.crs.to_epsg() == 4326 before any downstream export or API call expecting WGS84.
6. Stacked isochrone containment. For multi-ring outputs, each inner ring must be fully contained within the next: assert iso_5min.within(iso_10min).all(). Violations indicate buffering inconsistencies between bands.

Related

• Generating Isochrones with PySAL and GeoPandas — matrix-based accessibility pipeline
• NetworkX Dijkstra vs A* for route calculation — algorithm selection for logistics graphs
• Valhalla cost matrix generation for urban planners — multi-modal travel time matrices
• Implementing multi-modal transit layers in routing graphs