Building Isochrones from a Street Network

An isochrone answers "where can I get to in fifteen minutes", and the honest answer is a ragged, hole-punched shape dictated by which streets actually connect — never a circle. This guide is for analysts producing catchment or service-area polygons from an OSMnx graph: it sits under Street Network Analysis with OSMnx in Spatial Analysis & Advanced Query Techniques. The graph work is a handful of lines; the difficulty is entirely in the conversion from a set of reachable nodes to a polygon, which is where most published isochrones quietly become fiction.

Why This Approach / What Goes Wrong

nx.ego_graph gives you reachability for free: pass a weight of travel_time and a radius in seconds, and NetworkX runs a Dijkstra from the origin, keeping every node whose cumulative cost stays under the budget. That part is exact. What it returns, though, is a subgraph — a scatter of intersection nodes and the street segments between them — and a scatter of points is not a service area. Something has to fill the space between them, and the filling method decides whether the polygon is a measurement or a decoration.

The default reflex is convex_hull, and it is wrong almost everywhere. A convex hull is by definition the smallest shape with no dents, so every concavity a real street grid produces — a river with two bridges, a rail corridor, a golf course, a motorway with no crossing for a kilometre — gets paved over. A fifteen-minute walk isochrone drawn as a convex hull routinely claims 40–70% more area than the network supports, and every population figure derived from it inherits that inflation. It also lies in the other direction, because the hull's straight edges cut across the tips of long radial arms.

A concave hull (an alpha shape) fixes the gross over-coverage by letting the boundary dent inward, and GEOS 3.11 exposes it directly as shapely.concave_hull(geom, ratio=...). It is a real improvement, but the ratio parameter is a tuning knob with no physical meaning: too high and you are back to a convex hull, too low and the shape shatters into disconnected slivers around sparse arms. The tuned value never transfers between a dense city centre and a suburban fringe, so a batch of isochrones needs a batch of hand-tuned ratios.

The defensible method is to forget the nodes and polygonize the edges. Buffer every reachable street segment by a half-width that represents how far a person can stray from the centreline — 20–30 m for walking, roughly a block face — and dissolve the buffers into one geometry. The result covers exactly the streets you can travel and nothing else, keeps genuine holes where blocks are impassable, and its one parameter (the buffer distance) is a physical quantity you can defend in a report. The mechanics are ordinary proximity and buffer analysis; the difference is that the input is a reachable subgraph rather than the whole layer.

One reachable node set rendered three ways Three panels show the identical set of reachable intersection nodes and street segments, a ragged eight-armed star. The first panel wraps them in a convex hull, whose straight boundary shades over the whole area between the arms, including the blocks nobody can reach; the over-covered area is highlighted in warm red. The second panel wraps them in a concave hull, which dents inward between the arms but depends on a ratio parameter that has no physical meaning and does not transfer between neighbourhoods. The third panel buffers each reachable street segment by twenty-five metres and dissolves the buffers, producing a shape that covers only travelled streets, keeps the gaps between arms open, and has one physically defensible parameter. Same reachable nodes, three polygons, three different answers convex_hull concave_hull(ratio=0.15) buffered edges, dissolved shaded area = blocks nobody can reach +40 to +70% area on a real grid dents inward, but ratio is a free knob retune per neighbourhood density covers travelled streets only one parameter, in metres, defensible The node set is identical in all three panels — only the fill rule changes, and with it every area and population figure. Buffer half-width is a physical quantity (how far you stray from a centreline); a hull ratio is not.
The reachability computation is exact; the polygon you wrap around it is a modelling choice, and the convex hull is the one that inflates every downstream number.

Prerequisites

conda install -c conda-forge "osmnx=2.1.*" "networkx=3.3.*" "geopandas=1.0.*" "shapely=2.0.*" "pyproj=3.6.*"

Install from conda-forge in one solve: shapely.concave_hull raises UnsupportedGEOSVersionError when Shapely is linked against GEOS below 3.11, and a mixed pip/conda environment is the usual way to end up with an old GEOS under a new Shapely.

Step-by-Step Implementation

1. Download a graph with a margin, and project it before anything else.

Request more network than the band needs. If the download radius equals the reachable radius, the isochrone gets clipped flat against the bounding box and you will read the artefact as a real barrier. A rule of thumb: download at least 1.4× the straight-line distance the fastest mode could cover in the largest band.

import networkx as nx
import osmnx as ox

ox.settings.use_cache = True
ox.settings.cache_folder = "./osm_cache"

clinic = (45.0703, 7.6869)          # lat, lon — an urban clinic in Turin
walk_graph = ox.graph_from_point(
    clinic, dist=2000, dist_type="bbox", network_type="walk"
)
walk_graph = ox.truncate.largest_component(walk_graph, strong=True)
walk_graph = ox.projection.project_graph(walk_graph)     # local UTM zone, metres

print(walk_graph.graph["crs"])                       # EPSG:32632
print(walk_graph.number_of_nodes(), "nodes")         # 8342 nodes

Projection is not optional here even though length is already metric. Every other number in this pipeline — the buffer half-width, the simplification tolerance, the KD-tree distance inside nearest_nodes, the final area — is in CRS units, and in EPSG:4326 those are degrees. project_graph picks the UTM zone containing the graph's centroid. Do not substitute EPSG:3857: its scale factor at Turin's latitude inflates every distance by about 42%, so a 25 m buffer becomes 35 m on the ground and the reported area is out by roughly a factor of two.

2. Attach a travel time to every edge.

add_edge_travel_times divides length by speed_kph, so it needs a speed on every edge. For a driving graph, ox.routing.add_edge_speeds imputes one per highway class from the observed maxspeed tags. For a walking graph, do not use it: maxspeed describes cars, and imputing 50 km/h onto a footway produces a nine-kilometre "walk" isochrone. Set a constant pedestrian speed instead.

WALK_KPH = 4.5                       # ~1.25 m/s, a normal adult pace

for _, _, edge in walk_graph.edges(data=True):
    edge["speed_kph"] = WALK_KPH

walk_graph = ox.routing.add_edge_travel_times(walk_graph)   # writes travel_time, seconds

3. Snap the origin, then cut the reachable subgraph with ego_graph.

The graph is projected now, so the origin has to be projected too — and Transformer needs always_xy=True or PROJ hands back latitude first for EPSG:4326 and your clinic lands in the sea.

from pyproj import Transformer

to_metric = Transformer.from_crs("EPSG:4326", walk_graph.graph["crs"], always_xy=True)
origin_x, origin_y = to_metric.transform(clinic[1], clinic[0])    # lon, lat → x, y
origin_node = ox.distance.nearest_nodes(walk_graph, X=origin_x, Y=origin_y)

BUDGET_S = 15 * 60
reachable = nx.ego_graph(walk_graph, origin_node, radius=BUDGET_S, distance="travel_time")
print(reachable.number_of_nodes(), "intersections within 15 min")   # 1874 intersections within 15 min

radius is interpreted in the units of distance; omit distance= and NetworkX counts hops instead, silently returning a 900-edge neighbourhood. Note also the direction: on a MultiDiGraph, ego_graph follows out-edges, so this is "where can someone leave the clinic and reach". For an inbound catchment — "who can reach the clinic" — run it on walk_graph.reverse(copy=False).

The time-band pipeline from download to web map A flow in two rows that reverses direction between them. The top row runs left to right: graph_from_point downloads a walking network with a margin, project_graph converts it to a local UTM zone in metres, a constant speed plus add_edge_travel_times writes travel_time in seconds per edge, and single_source_dijkstra_path_length returns seconds from the origin for every node. The flow drops to the bottom row and runs right to left: partial edges are trimmed with substring to the remaining budget, the trimmed segments are buffered and dissolved with union_all, each band has the next smaller band differenced away so rings do not stack, and the result is reprojected to EPSG 4326 and written as GeoJSON for a web map. A label notes that the last four stages repeat once per time band. One band, eight stages — four of them repeat per band graph_from_point walk network + 1.4x margin project_graph local UTM · metres add_edge_travel_times seconds on every edge dijkstra_path_length seconds per node substring(0, frac) trim the boundary edges buffer + union_all 25 m half-width, dissolved difference(inner) rings, not stacked discs to_crs(4326) GeoJSON → web map repeat for band in [300 s, 600 s, 900 s] — the graph is built once Only the shaded cyan stages touch geometry; everything before them is graph arithmetic in seconds.
The download and projection happen once; the Dijkstra, the trim, the dissolve, and the difference run once per time band.

4. Trim the partial edges the subgraph cuts off.

ego_graph stops at nodes, so an edge whose far intersection sits one second past the budget is dropped whole. On a walking graph with 120 m blocks that removes up to a block face all the way around the boundary — a systematic under-count, visible as a boundary that stops short of every corner. The fix is to run the Dijkstra yourself, keep the per-node cost dictionary, and clip each boundary edge at the fraction of its length the remaining budget buys.

import geopandas as gpd
from shapely.geometry import Point
from shapely.ops import substring

times = nx.single_source_dijkstra_path_length(walk_graph, origin_node, weight="travel_time")
street_edges = ox.convert.graph_to_gdfs(walk_graph, nodes=False).reset_index()

def reachable_segments(budget_s: float) -> gpd.GeoSeries:
    """Whole edges inside the budget, plus the reachable fraction of boundary edges."""
    pieces = []
    for edge in street_edges.itertuples():
        t_u = times.get(edge.u)
        if t_u is None or t_u >= budget_s:
            continue                              # u itself is out of reach
        geom = edge.geometry
        u_pt = Point(walk_graph.nodes[edge.u]["x"], walk_graph.nodes[edge.u]["y"])
        if Point(geom.coords[0]).distance(u_pt) > Point(geom.coords[-1]).distance(u_pt):
            geom = geom.reverse()                 # orient the line so it starts at u
        if times.get(edge.v, float("inf")) <= budget_s:
            pieces.append(geom)                   # both ends inside — keep it whole
        elif edge.travel_time > 0:
            frac = min((budget_s - t_u) / edge.travel_time, 1.0)
            pieces.append(substring(geom, 0.0, frac, normalized=True))
    return gpd.GeoSeries(pieces, crs=walk_graph.graph["crs"])

Orienting the geometry matters: OSMnx stores one LineString per directed edge, and for the reciprocal edge of a two-way street the coordinate order can run v → u. Trimming an unoriented line takes the fraction from the wrong end and grows a spur outward from the boundary instead of inward.

Node-cut boundary versus edge-trimmed boundary Two panels show the same street segment crossing a fifteen-minute boundary. Node A is reached at eight hundred and forty seconds and node B at nine hundred and sixty, so B is outside the nine hundred second budget. In the left panel ego_graph drops the whole segment, and the isochrone boundary snaps back to node A, losing sixty seconds of walking, about seventy-five metres, all the way around the perimeter. In the right panel the remaining sixty seconds are divided by the edge travel time of one hundred and twenty seconds to give a fraction of one half, and shapely ops substring clips the line at its midpoint, so the boundary sits where the walker actually runs out of time. The boundary lands at a node, or where the clock actually runs out ego_graph alone · node cut substring trim · edge cut node A · 840 s node B · 960 s travel_time = 120 s whole edge dropped — 60 s of walking lost ≈75 m short, everywhere on the perimeter node A · 840 s node B · 960 s 900 s frac = (900 - 840) / 120 = 0.5 substring(geom, 0.0, 0.5, normalized=True) boundary sits at the real 15-minute point The error is one-sided: a node-cut isochrone is always too small, never too large.
Clipping boundary edges by the remaining budget removes a systematic shortfall of up to one block face around the entire perimeter.

5. Dissolve the segments into a polygon.

Buffer, dissolve, then close pinholes with a dilate–erode pass so that a bundle of parallel one-way pairs reads as one solid block rather than a lattice of hairline gaps.

def isochrone_polygon(budget_s: float, half_width_m: float = 25.0):
    segments = reachable_segments(budget_s)
    covered = segments.buffer(half_width_m, resolution=8).union_all()
    return covered.buffer(30).buffer(-30)      # close slivers, keep real holes

band_900 = isochrone_polygon(900)
print(f"{band_900.area / 1e6:.2f} km²")        # 1.94 km²

6. Build nested bands that do not stack.

Compute the largest band first, then union each band with the one inside it to guarantee strict nesting, and finally subtract the inner band from each ring. Without the subtraction the polygons overlap, and a semi-transparent web-map fill renders the inner bands three shades darker than the legend claims.

BANDS_S = [300, 600, 900]

filled = {}
for band in sorted(BANDS_S):
    poly = isochrone_polygon(band)
    inner = filled.get(band - 300)
    filled[band] = poly.union(inner) if inner is not None else poly   # enforce nesting

rings, previous = [], None
for band in sorted(BANDS_S):
    ring = filled[band].difference(previous) if previous is not None else filled[band]
    rings.append({"band_min": band // 60, "geometry": ring})
    previous = filled[band]

isochrones = gpd.GeoDataFrame(rings, crs=walk_graph.graph["crs"])
print(isochrones.assign(km2=isochrones.area / 1e6)[["band_min", "km2"]])
#    band_min       km2
# 0         5  0.243611
# 1        10  0.771044
# 2        15  0.928472
Nested five, ten and fifteen minute walk bands against a straight-line circle A map-style figure. Three nested ragged star polygons radiate from a clinic marker: a dark teal five-minute band, a cyan ten-minute band, and a pale fifteen-minute band, each one reaching further along the streets and dented inward between them. A dashed red circle of eleven hundred and twenty-five metres, the straight-line distance a walker covers in fifteen minutes at four and a half kilometres per hour, encloses all three and covers roughly four times the area of the fifteen-minute band. A legend on the right lists each band with its area in square kilometres and notes that the rings are differenced so the fills do not stack. Nested walk bands from one clinic — and the circle that would have lied dashed circle = 1125 m as the crow flies Walk band · 4.5 km/h 0–5 min 0.24 km² 5–10 min 0.77 km² 10–15 min 0.93 km² total reachable 1.94 km² 1125 m circle 3.98 km² rings are differenced — fills never stack The circle is a 2.1× over-estimate here, and it is wrong in a different direction in every neighbourhood you try it in.
Each band is differenced against the one inside it, so the three fills tile the service area exactly once — and the dashed circle shows what a straight-line buffer would have claimed.

7. Export for a web map.

Simplify while still in metres, then reproject: RFC 7946 GeoJSON is defined in WGS 84, and a simplify tolerance applied after reprojection would be in degrees.

web_layer = isochrones.copy()
web_layer["geometry"] = web_layer.geometry.simplify(8)      # 8 metres, not 8 degrees
web_layer = web_layer.to_crs("EPSG:4326")
web_layer.to_file("clinic_isochrones.geojson", driver="GeoJSON")

Draw the bands largest-first so the small ones are not hidden underneath; the style_function pattern for graduated fills is covered in Interactive Maps with Folium.

Verification

The strongest check is physical: no isochrone can escape the disc a walker could cover in a straight line at the same speed. If it does, a speed or a unit is wrong.

from shapely.geometry import Point

assert isochrones.crs.is_projected, "Areas and buffers are meaningless in degrees"
assert isochrones.geometry.is_valid.all(), "Invalid ring after difference()"

# Physical bound: 15 min at 4.5 km/h = 1125 m as the crow flies
crow_fly = Point(origin_x, origin_y).buffer(WALK_KPH * 1000 / 3600 * BUDGET_S)
assert filled[900].within(crow_fly.buffer(1.0)), "Isochrone escapes its own speed limit"

# Bands must nest, and each ring must add area
assert filled[300].within(filled[600]) and filled[600].within(filled[900])
assert isochrones.geometry.area.min() > 0

ratio = filled[900].area / crow_fly.area
print(f"network area {filled[900].area/1e6:.2f} km² · {ratio:.0%} of the circle")
# network area 1.94 km² · 49% of the circle

A ratio near 100% means the polygon has been over-filled — usually a convex hull, an oversized buffer half-width, or a graph so sparse that the dilate–erode pass welded separate arms together. Below about 20% the graph is probably fragmented; check that largest_component ran and that the origin snapped onto a routable street rather than a disconnected service lane.

Edge Cases & Debugging

Frequently Asked Questions

Which method should I ship: convex hull, concave hull, or buffered edges? Buffered edges, unless you have a specific reason not to. It is the only one whose parameter has a physical meaning, it preserves genuine barriers, and it degrades gracefully in both dense and sparse networks. Use concave_hull when a downstream consumer demands a single simply-connected polygon with no holes and you can tune ratio per area. Use convex_hull only for a rough visual index, and never for area or population figures.

How do I do this for driving instead of walking? Download with network_type="drive", then use ox.routing.add_edge_speeds(streets, fallback=40) instead of a constant speed, because maxspeed tags are meaningful for cars. Raise the buffer half-width to 40–60 m to reflect a road corridor, and expect strong asymmetry between inbound and outbound bands from one-way systems. Free-flow maxspeed also means the result is an off-peak isochrone; a peak-hour version needs observed speeds, which OSM does not carry. The routing side of that setup is covered in shortest path routing with NetworkX and OSMnx.

Why is my isochrone smaller than a competitor's for the same location? Three usual causes, in order of frequency: their polygon is a convex hull; they never trimmed partial edges, which under-counts, or they buffered by 100 m, which over-counts; or they used a faster assumed speed. Publish the speed, the buffer half-width, and the polygon method alongside the shape — two isochrones without those three numbers are not comparable.

Can I get population inside each band? Yes, and this is where the polygon method starts to matter financially. Overlay the differenced rings onto census or gridded-population polygons and area-weight the intersection, or sample a population raster per ring. Because a convex hull inflates the area by 40–70%, it inflates the served population by a similar margin — which is exactly the number that ends up in a funding application.

Related