Source code: Street network analysis with OSMnx: routes, isochrones and walking accessibility

"""Street network analysis with OSMnx - TechnicalGIS tutorial.

Downloads Sheffield's walkable street network from OpenStreetMap, then
computes a shortest path, walking isochrones and a network-based
accessibility surface (minutes to the nearest supermarket).
Run with:  python sheffield_networks.py
"""
from pathlib import Path

import contextily as cx
import geopandas as gpd
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import osmnx as ox
import pandas as pd
from shapely.geometry import Point

HERE = Path(__file__).parent
OUT = HERE / "output"
OUT.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 120)
ox.settings.use_cache = True
ox.settings.cache_folder = str(HERE / "cache")
ox.settings.log_console = False

CENTRE = (53.3801, -1.4703)          # Peace Gardens, Sheffield (lat, lon)
STATION = (53.3780, -1.4623)         # Sheffield railway station
UNIVERSITY = (53.3817, -1.4886)      # Firth Court, University of Sheffield
WALK_M_PER_MIN = 4.8 * 1000 / 60     # 4.8 km/h
BASEMAP = cx.providers.Esri.WorldGrayCanvas


def fill_small_holes(geom, max_area=40_000):
    """Remove interior rings (city blocks) smaller than max_area m²."""
    from shapely.geometry import MultiPolygon, Polygon
    polys = geom.geoms if isinstance(geom, MultiPolygon) else [geom]
    out = [Polygon(p.exterior, [r for r in p.interiors if Polygon(r).area > max_area]) for p in polys]
    return MultiPolygon(out) if len(out) > 1 else out[0]


def section(title: str) -> None:
    print(f"\n=== {title} ===")


# --- Step 1: download the walking network -------------------------------------
G = ox.graph_from_point(CENTRE, dist=2500, network_type="walk")
Gp = ox.project_graph(G)  # to local UTM so lengths/areas are in metres
section("graph")
print(Gp)
print("crs:", Gp.graph["crs"])
stats = ox.basic_stats(Gp)
for k in ["n", "m", "k_avg", "edge_length_total", "street_length_avg", "streets_per_node_avg",
          "intersection_count", "circuity_avg"]:
    v = stats[k]
    print(f"{k:22s} {v:,.2f}" if isinstance(v, float) else f"{k:22s} {v:,}")

nodes, edges = ox.graph_to_gdfs(Gp)
section("edges GeoDataFrame")
print(edges[["osmid", "name", "highway", "length"]].head().to_string())
print("geometry types:", edges.geom_type.value_counts().to_dict())

fig, ax = ox.plot_graph(Gp, node_size=0, edge_color="#0f766e", edge_linewidth=0.5, bgcolor="white",
                        show=False, close=False, figsize=(8, 8))
ax.set_title("Sheffield walking network within 2.5 km of the Peace Gardens")
fig.savefig(OUT / "01-network.png")
plt.close(fig)

# --- Step 2: shortest path -----------------------------------------------------------
pts = gpd.GeoSeries([Point(STATION[1], STATION[0]), Point(UNIVERSITY[1], UNIVERSITY[0])],
                    crs="EPSG:4326").to_crs(Gp.graph["crs"])
orig, dest = ox.distance.nearest_nodes(Gp, pts.x.tolist(), pts.y.tolist())
route = ox.shortest_path(Gp, orig, dest, weight="length")
route_gdf = ox.routing.route_to_gdf(Gp, route, weight="length")
route_len = route_gdf["length"].sum()
section("shortest path station -> university")
print("nodes in route:", len(route))
print(f"length: {route_len:,.0f} m  |  walking time: {route_len / WALK_M_PER_MIN:.1f} min")
print("streets:", " → ".join(pd.unique(route_gdf["name"].dropna().astype(str)))[:400])

fig, ax = ox.plot_graph_route(Gp, route, route_color="#e8590c", route_linewidth=4, node_size=0,
                              edge_color="#9aa5a1", edge_linewidth=0.4, bgcolor="white",
                              show=False, close=False, figsize=(8, 6))
ax.set_title(f"Shortest walk: station to Firth Court ({route_len / 1000:.2f} km)")
fig.savefig(OUT / "02-shortest-path.png")
plt.close(fig)

# --- Step 3: walking isochrones --------------------------------------------------------
times = [5, 10, 15, 20]
isos = []
for t in sorted(times, reverse=True):
    sub = nx.ego_graph(Gp, orig, radius=t * WALK_M_PER_MIN, distance="length")
    sub_edges = ox.graph_to_gdfs(sub, nodes=False)
    isos.append({"minutes": t, "geometry": fill_small_holes(sub_edges.buffer(30).unary_union.buffer(-10))})
iso = gpd.GeoDataFrame(isos, crs=Gp.graph["crs"])
iso["area_km2"] = (iso.area / 1e6).round(2)
section("isochrones from Sheffield station")
print(iso[["minutes", "area_km2"]].sort_values("minutes").to_string(index=False))
iso.to_crs("EPSG:4326").to_file(OUT / "isochrones.geojson", driver="GeoJSON")

fig, ax = plt.subplots(figsize=(8, 8))
iso.plot(ax=ax, column="minutes", cmap="viridis_r", alpha=0.55, edgecolor="white", linewidth=0.6,
         legend=True, categorical=True, legend_kwds={"title": "Walk (min)", "loc": "lower right"})
pts.iloc[[0]].plot(ax=ax, color="#e8590c", markersize=60, marker="*", zorder=5)
cx.add_basemap(ax, crs=iso.crs, source=BASEMAP, attribution_size=6)
ax.set_title("5–20 minute walking isochrones from Sheffield station")
ax.set_axis_off()
fig.savefig(OUT / "03-isochrones.png")
plt.close(fig)

# --- Step 4: accessibility to supermarkets ------------------------------------------
shops = ox.features_from_point(CENTRE, tags={"shop": "supermarket"}, dist=3000)
shops = shops.to_crs(Gp.graph["crs"])
shops["geometry"] = shops.geometry.representative_point()
section("supermarkets from OSM")
print(len(shops), "supermarkets;", shops["name"].value_counts().head(6).to_dict())

shop_nodes = set(ox.distance.nearest_nodes(Gp, shops.geometry.x.tolist(), shops.geometry.y.tolist()))
dist = nx.multi_source_dijkstra_path_length(Gp, shop_nodes, weight="length")
nodes["walk_min"] = nodes.index.map(dist).astype(float) / WALK_M_PER_MIN
edges["walk_min"] = [np.nanmean([nodes.at[u, "walk_min"], nodes.at[v, "walk_min"]]) for u, v, _ in edges.index]

bins = pd.cut(nodes["walk_min"], [0, 5, 10, 15, 20, np.inf], labels=["0–5", "5–10", "10–15", "15–20", "20+"])
section("share of network nodes by walk time to nearest supermarket")
print((bins.value_counts(normalize=True).sort_index() * 100).round(1).to_string())

fig, ax = plt.subplots(figsize=(8, 8))
edges.plot(ax=ax, column="walk_min", cmap="plasma", vmin=0, vmax=20, linewidth=0.8, legend=True,
           legend_kwds={"label": "Walking minutes to nearest supermarket", "shrink": 0.7})
shops.plot(ax=ax, color="black", markersize=12, marker="s")
ax.set_title("Network accessibility: minutes to the nearest supermarket")
ax.set_axis_off()
fig.savefig(OUT / "04-supermarket-accessibility.png")
plt.close(fig)

# --- Step 5: street orientation ------------------------------------------------------
Gu = ox.bearing.add_edge_bearings(ox.convert.to_undirected(G))
fig, ax = ox.plot_orientation(Gu, title="Sheffield centre street orientation", area=True,
                              color="#0f766e", edgecolor="#0b4f49",
                              title_font={"family": "DejaVu Sans", "size": 14, "weight": "bold"})
fig.savefig(OUT / "05-orientation.png")
plt.close(fig)
section("orientation entropy")
print(round(ox.bearing.orientation_entropy(Gu), 3))

ox.save_graph_geopackage(Gp, OUT / "sheffield_walk.gpkg")