Generating Isochrones with PySAL and GeoPandas

Isochrone mapping converts raw network travel times into precise spatial polygons that answer the question “where can I reach from here within N minutes?” The problem sits at the intersection of graph traversal and spatial analysis: a shortest-path algorithm produces discrete travel times at network nodes, but service-area planning requires continuous, exportable polygon boundaries. This page documents a reproducible, open-source pipeline that solves that gap using PySAL’s spatial analysis utilities and GeoPandas’ vector operations — part of the broader Python Routing Engines & Isochrone Mapping toolkit for logistics engineers and GIS developers who need auditable, locally computed service areas.

Prerequisites

Before building the pipeline, confirm your environment meets these requirements:

• Python 3.9+ inside an isolated virtual environment (venv or conda). C-extension wheels for GDAL and libspatialindex require careful platform alignment.
• Core libraries: geopandas, libpysal, networkx, osmnx, shapely>=2.0, scipy, pandas, numpy, scikit-image
• Routing backend: A graph source capable of returning travel-time edge weights. For offline capability and data-residency compliance, deploying OSRM with Docker for local routing eliminates external API dependencies entirely and provides a fast HTTP interface if you need to supplement OSMnx-based queries with pre-contracted road networks.
• Spatial data: OSM PBF extract or GraphML file for the target region, projected to a consistent metric CRS (EPSG:3857 or local UTM zone).

# Install the full isochrone stack
pip install geopandas libpysal networkx osmnx "shapely>=2.0" scipy numpy scikit-image

Confirm GeoPandas can load a test layer before proceeding — GDAL driver availability varies by platform and silent failures during CRS operations are common when system GDAL and Python bindings are mismatched.

Conceptual Architecture

The pipeline bridges two different representations of accessibility: the discrete graph domain (nodes with travel-time labels) and the continuous raster domain (a regular grid suitable for contouring). Three key decisions govern quality:

Impedance model. Node travel times are only as accurate as the edge weights feeding nx.single_source_dijkstra. OSMnx sets travel_time in seconds on each edge using maxspeed tags and free-flow assumptions. For freight or delivery scenarios, replace those weights with speed profiles calibrated to vehicle class — see the discussion of speed profile calibration for heavy vehicles for the tag-to-weight mapping that applies here.

Interpolation method. scipy.interpolate.griddata with method='linear' fits a piecewise linear surface through the scattered node samples. This works well for dense urban networks but degrades in low-density areas where node spacing is wide relative to grid resolution. Setting fill_value to cutoff_sec + 1 ensures cells outside the convex hull of sampled nodes are treated as unreachable rather than extrapolated.

Contour extraction. skimage.measure.find_contours operates directly on the numpy grid and returns closed contour paths as (row, col) arrays. This is the correct headless approach — matplotlib’s ContourSet.collections attribute was removed in matplotlib 3.10, making ax.contour() unreliable in server environments.

The diagram above shows how these three subsystems connect. Disconnected graph components must be filtered at step ② or they produce spurious islands in the interpolated grid at step ③.

Step-by-Step Implementation

Step 1: Graph Extraction and Impedance Assignment

# Requires: osmnx, networkx
import osmnx as ox
import networkx as nx

def build_weighted_graph(
    lat: float,
    lon: float,
    dist_m: float = 5000,
    network_type: str = "drive"
) -> nx.MultiDiGraph:
    """Extract an OSMnx drive network and assign travel-time edge weights."""
    G = ox.graph_from_point((lat, lon), dist=dist_m, network_type=network_type)
    G = ox.add_edge_speeds(G)          # fills maxspeed from OSM tags + heuristic fallback
    G = ox.add_edge_travel_times(G)    # computes travel_time = length / speed
    G_proj = ox.project_graph(G)       # project to UTM for metric operations
    return G_proj

ox.add_edge_speeds uses maxspeed tag values where present and applies road-class speed heuristics elsewhere. For accurate freight routing, override those heuristics with vehicle-class speed maps before calling ox.add_edge_travel_times. Retain the unprojected graph G separately — ox.distance.nearest_nodes requires longitude/latitude in WGS84 and must be called against the unprojected version.

Step 2: Single-Source Dijkstra Traversal

# Requires: networkx, osmnx
import numpy as np

def compute_travel_times(
    G_proj: nx.MultiDiGraph,
    G_wgs: nx.MultiDiGraph,
    origin_lon: float,
    origin_lat: float,
    cutoff_min: float
) -> dict:
    """Return {node_id: travel_time_sec} for nodes reachable within cutoff_min."""
    # nearest_nodes expects (lon, lat) on the *unprojected* WGS84 graph
    origin_node = ox.distance.nearest_nodes(G_wgs, origin_lon, origin_lat)
    cutoff_sec = cutoff_min * 60.0

    # single_source_dijkstra returns (lengths_dict, paths_dict)
    lengths, _ = nx.single_source_dijkstra(
        G_proj, origin_node, cutoff=cutoff_sec, weight="travel_time"
    )
    return {n: t for n, t in lengths.items() if t <= cutoff_sec}

Passing cutoff=cutoff_sec to single_source_dijkstra short-circuits traversal once the frontier exceeds the threshold, which cuts runtime significantly on large graphs. Always filter the result dict anyway — floating-point accumulation can leave a handful of nodes marginally above the threshold.

Step 3: Grid Interpolation of Scattered Travel Times

# Requires: numpy, scipy
from scipy.interpolate import griddata

def interpolate_travel_time_grid(
    G_proj: nx.MultiDiGraph,
    reachable: dict,
    cutoff_sec: float,
    resolution_m: int = 50
) -> tuple:
    """Interpolate node travel times onto a regular metric grid.

    Returns (grid_z, minx, miny, x_step, y_step) for downstream contour extraction.
    """
    x_coords = np.array([G_proj.nodes[n]["x"] for n in reachable])
    y_coords = np.array([G_proj.nodes[n]["y"] for n in reachable])
    times    = np.array(list(reachable.values()))

    minx, maxx = x_coords.min(), x_coords.max()
    miny, maxy = y_coords.min(), y_coords.max()

    grid_x, grid_y = np.meshgrid(
        np.arange(minx, maxx, resolution_m),
        np.arange(miny, maxy, resolution_m)
    )

    grid_z = griddata(
        (x_coords, y_coords),
        times,
        (grid_x, grid_y),
        method="linear",
        fill_value=cutoff_sec + 1.0   # mark extrapolated cells as unreachable
    )

    x_step = (maxx - minx) / grid_z.shape[1]
    y_step = (maxy - miny) / grid_z.shape[0]
    return grid_z, minx, miny, x_step, y_step

resolution_m controls the fidelity/performance trade-off. At 50 m in a 5 km radius study area, the grid contains roughly 40,000 cells — fast enough for synchronous service. For metropolitan-scale analysis (20+ km radius), raise resolution_m to 100–200 m or switch to dask-geopandas for tiled processing.

Step 4: Contour Extraction with skimage

# Requires: scikit-image, shapely
from skimage.measure import find_contours
from shapely.geometry import Polygon

def extract_isochrone_polygons(
    grid_z: np.ndarray,
    cutoff_sec: float,
    minx: float,
    miny: float,
    x_step: float,
    y_step: float
) -> list:
    """Convert grid contours at cutoff_sec into Shapely Polygons in projected coords."""
    # find_contours returns list of (N, 2) arrays in (row, col) order
    contour_paths = find_contours(grid_z, level=cutoff_sec)

    polygons = []
    for contour in contour_paths:
        # contour[:, 1] = column index → x axis
        # contour[:, 0] = row index    → y axis
        world_x = minx + contour[:, 1] * x_step
        world_y = miny + contour[:, 0] * y_step
        coords = list(zip(world_x, world_y))
        if len(coords) >= 4:          # Polygon requires ≥ 4 points (first = last)
            polygons.append(Polygon(coords))
    return polygons

find_contours uses a sub-pixel marching-squares algorithm that produces smoother boundaries than integer-grid approaches. The row/column-to-world-coordinate conversion is straightforward but easy to swap — verify by checking that the polygon centroid is close to the origin point after reprojecting to WGS84.

Step 5: Topology Validation and Export

# Requires: geopandas, shapely
import geopandas as gpd
from shapely.validation import make_valid

def build_isochrone_gdf(
    polygons: list,
    crs: str,
    max_time_min: float
) -> gpd.GeoDataFrame:
    """Merge and validate raw contour polygons into a clean GeoDataFrame."""
    if not polygons:
        return gpd.GeoDataFrame(columns=["geometry", "max_time_min"], crs=crs)

    gdf = gpd.GeoDataFrame(geometry=polygons, crs=crs)

    # Step 1: repair per-geometry issues (self-intersections, degenerate rings)
    gdf["geometry"] = gdf["geometry"].apply(lambda g: make_valid(g).buffer(0))

    # Step 2: dissolve all fragments into a single polygon (shapely 2.x API)
    merged = gdf.union_all()

    # Step 3: optionally replace convex envelope with concave hull for realism
    # merged = merged.concave_hull(ratio=0.3)   # shapely 2.x; tune ratio by network density

    return gpd.GeoDataFrame(
        {"geometry": [merged], "max_time_min": [max_time_min]},
        crs=crs
    )

buffer(0) corrects most self-intersections produced by interpolation noise. make_valid() handles edge cases that buffer(0) misses, including zero-area rings and repeated points. Apply both in sequence before union_all() — calling union_all() on an invalid geometry raises a TopologicalError that is difficult to trace back to its source.

Configuration Reference

OSM tag values that affect the impedance model: maxspeed, highway, access, oneway. When maxspeed is absent on highway=residential edges, OSMnx assigns 35 km/h by default. Override this in the speed fallback dictionary passed to ox.add_edge_speeds(G, fallback=...).

Production Optimization and Scaling

Precompute for fixed-origin grids. Logistics depots, hospitals, and transit hubs have stable origin points. Precomputing isochrones at multiple thresholds (5, 10, 15, 30 min) and storing them in PostGIS cuts response latency to a spatial lookup rather than a full graph traversal.

Cache with compound keys. Key a Redis cache on (geohash6_of_origin, cutoff_sec, network_type). A geohash at precision 6 resolves to ~1.2 km cells, which provides acceptable locality for service-area queries without over-fragmenting the cache. Serialise GeoDataFrames to GeoJSON or FlatGeobuf — GeoJSON is simpler; FlatGeobuf is faster for large polygons.

Chunked tiling for metro scale. For 20+ km radius analyses, partition the study area into tiles and run interpolation per-tile with a 20% overlap buffer to prevent edge discontinuities. Merge tile results with gpd.overlay using how='union' and dissolve on the threshold column.

Async graph traversal. If you serve dynamic origins over HTTP, run the Dijkstra step in a ProcessPoolExecutor with a pre-loaded projected graph object. Loading and projecting the OSMnx graph takes 2–8 seconds depending on area; the traversal itself is typically under 200 ms for a 5 km radius on a single core.

Replace griddata for large grids. scipy.interpolate.griddata builds a full Delaunay triangulation each call — it does not cache the triangulation. For repeated calls over the same network, switch to scipy.interpolate.LinearNDInterpolator, which separates the triangulation step (run once at startup) from the interpolation step (run per request):

# Requires: scipy, numpy
from scipy.interpolate import LinearNDInterpolator

# Build once at startup (expensive)
interp = LinearNDInterpolator(list(zip(x_coords, y_coords)), times, fill_value=cutoff_sec+1)

# Call cheaply per request
grid_z = interp(grid_x, grid_y)

Validation and Testing

After generating an isochrone, apply these checks before trusting the polygon for service-area decisions:

Centroid proximity. Reproject the merged polygon to WGS84 and confirm its centroid is within 500 m of the origin point. A centroid far from the origin indicates a grid coordinate transform error — usually a swapped lat/lon or a CRS mismatch between the snapped node and the grid origin.

Area plausibility. For a 15-minute drive-time isochrone in a medium-density urban network, expected area is typically 15–50 km². Values outside that range suggest a resolution or fill_value error.

Validity check. gdf.is_valid.all() must return True after topology repair. Any False entry will cause silent geometry failures in downstream PostGIS inserts or gpd.overlay calls.

Boundary smoothness. Plot the polygon exterior with gdf.boundary.plot(). Staircase artefacts along grid edges mean resolution_m is too coarse for the analysis scale. Jagged spikes indicate unreachable nodes that were not filtered before interpolation.

Multi-threshold consistency. Generate isochrones at 5, 10, and 15 minutes from the same origin. Each polygon must be a subset of the next — a 5-minute boundary that extends beyond the 10-minute boundary is a contour-level bug.

Regression on known routes. Select an origin-destination pair whose network travel time you can verify independently (e.g., from OSRM or Valhalla). Confirm the destination falls inside the matching threshold polygon. For 15-minute city isochrone analysis, this check directly validates the accessibility claims driving urban policy decisions.

Troubleshooting

Related

• Deploying OSRM with Docker for local routing — set up a self-hosted routing engine that eliminates external API dependencies for travel-time queries
• Creating 15-minute city isochrones in Python — apply this pipeline to urban accessibility metrics and policy-facing spatial outputs
• Valhalla configuration for multi-modal analysis — extend the impedance model to pedestrian, transit, and bicycle routing with per-mode cost functions
• Custom cost functions for routing solvers — design edge-weight functions that feed accurate travel times into the Dijkstra step
• NetworkX shortest-path algorithms for logistics — compare Dijkstra, A*, and Bellman-Ford for directed graph traversal in freight routing contexts