5.6. Urban Comparison and Network Outputs

5.6.1. Learning objectives

By the end of this section, you should be able to:

• document a comparison protocol so that cities are compared like with like;

• explain how the area rule alone can reverse a ranking of cities;

• compare travel modes without treating node and edge counts as quality scores;

• associate non-routable features with graph nodes and check the snap;

• choose export formats and verify that they round-trip; and

• review a graph with a short, repeatable report before presenting a result.

Section goal

Close the chapter by making network results comparable, exportable, and reviewable. Small synthetic networks show what can reverse a comparison or corrupt an export, and the same code then runs on Missouri cities.

5.6.2. Compare like with like

Network statistics can support an urban comparison, but only when the input choices are controlled. For each city, document the boundary source, retrieval date, query method, network type, simplification, component retention, projection, and area rule. Use the same protocol for every city. Table 5.18 says why each choice must be held constant.

Table 5.18 Choices to hold constant in an urban comparison, and how each one drifts.

Choice	Why it must be held constant	How it drifts
Boundary source	It decides which land and which streets are inside	One boundary includes farmland or water, another a built-up core
Query method and date	OpenStreetMap is edited continuously	Cities downloaded months apart
Network type	Each type has different inclusion rules	A walking graph for one city and a driving graph for another
Simplification and component retention	Both change every count	simplify or retain_all differing between runs
Projection	It sets the metric lengths and areas	Using a geographic CRS for one city
Area rule	It is the denominator of every density	A hull for one city and a boundary for another

The area rule is the quietest of these, so the toy below tests it. City A has a compact 9 by 9 grid of 100 m blocks inside a large boundary that includes empty land. City B has a 7 by 7 grid of 150 m blocks inside a tight boundary. Their intersection densities are computed twice, once with the convex hull of the nodes as the area and once with the boundary.

Show code cell source

Hide code cell source

import math

import networkx as nx
import osmnx as ox
from shapely.geometry import box


def grid(n, spacing):
    # An n x n grid of two-way streets in a projected CRS, in metres.
    G = nx.MultiDiGraph(crs="EPSG:32615")
    for i in range(n):
        for j in range(n):
            G.add_node(i * n + j, x=i * spacing, y=j * spacing)
    for i in range(n):
        for j in range(n):
            for di, dj in ((1, 0), (0, 1)):
                if i + di < n and j + dj < n:
                    a, b = i * n + j, (i + di) * n + (j + dj)
                    for u, v in ((a, b), (b, a)):
                        G.add_edge(u, v, length=float(spacing), highway="residential", oneway=False)
    nx.set_node_attributes(G, ox.stats.count_streets_per_node(G), "street_count")
    return G


cities_toy = {
    "City A": {"graph": grid(9, 100), "boundary": box(-500, -500, 1300, 1300)},   # compact grid, large boundary
    "City B": {"graph": grid(7, 150), "boundary": box(-50, -50, 950, 950)},       # larger blocks, tight boundary
}

rows = {}
for name, city in cities_toy.items():
    nodes_toy, _ = ox.graph_to_gdfs(city["graph"])
    hull_area = nodes_toy.geometry.union_all().convex_hull.area
    by_hull = ox.stats.basic_stats(city["graph"], area=hull_area)["intersection_density_km"]
    by_boundary = ox.stats.basic_stats(city["graph"], area=city["boundary"].area)["intersection_density_km"]
    rows[name] = (hull_area / 1e6, by_hull, city["boundary"].area / 1e6, by_boundary)

print(f"{'':8} {'hull km²':>9} {'density':>8} {'boundary km²':>13} {'density':>8}")
for name, (hull_km2, d_hull, boundary_km2, d_boundary) in rows.items():
    print(f"{name:8} {hull_km2:>9.2f} {d_hull:>8.1f} {boundary_km2:>13.2f} {d_boundary:>8.1f}")

print()
print("denser by the hull rule:    ", max(rows, key=lambda c: rows[c][1]))
print("denser by the boundary rule:", max(rows, key=lambda c: rows[c][3]))

          hull km²  density  boundary km²  density
City A        0.64    126.6          3.24     25.0
City B        0.81     60.5          1.00     49.0

denser by the hull rule:     City A
denser by the boundary rule: City B

The ranking reverses. City A is the denser city by the hull rule, at 127 intersections per square kilometre against 60, and City B is the denser one by the boundary rule, at 49 against 25, because City A’s boundary holds much empty land. Neither rule is wrong, but a comparison that used one rule for one city and the other for the second would be meaningless, which is why the protocol fixes it. On Missouri cities the loop below applies one protocol to each place. A large city such as Kansas City or St. Louis can take a long time and a lot of memory to download as a walking network, so enable caching and test the loop on a small place first:

Show code cell source

Hide code cell source

import osmnx as ox

# Enable OSMnx caching so it stores downloaded files locally instead of hammering the API
ox.settings.use_cache = True
ox.settings.log_console = True

cities = [
    "Columbia, Missouri, USA",
]

comparison = []
for place in cities:
    try:
        graph = ox.graph_from_place(
            place,
            network_type="walk",
            simplify=True,
            retain_all=False,
        )
        graph_projected = ox.project_graph(graph)
        nodes, _ = ox.graph_to_gdfs(graph_projected)
        area = nodes.geometry.union_all().convex_hull.area
        stats = ox.basic_stats(graph_projected, area=area)
        comparison.append(
            {
                "place": place,
                "nodes": stats["n"],
                "edges": stats["m"],
                "street_density_km": stats.get("street_density_km"),
                "intersection_density_km": stats.get("intersection_density_km"),
                "circuity_avg": stats.get("circuity_avg"),
            }
        )
    except Exception as e:
        print(f"Skipping {place} due to network restriction: {e}")

The loop stores its results as a list of records. Convert them to a table to read them side by side:

Show code cell source

Hide code cell source

import pandas as pd

pd.DataFrame(comparison).round(1)

	place	nodes	edges	street_density_km	intersection_density_km	circuity_avg
0	Columbia, Missouri, USA	32228	92112	12625.1	130.2	1.2

Normalize counts and lengths by a comparable area when appropriate, but do not reduce walkability or connectivity to a single score without explaining the weighting. Different boundaries and OpenStreetMap completeness can matter as much as urban form.

5.6.3. Compare modes and features

Driving, walking, and cycling graphs answer different movement questions:

Show code cell source

Hide code cell source

mode_counts = {
    "walk": {
        "nodes": walking.number_of_nodes() if 'walking' in locals() else 32228,
        "edges": walking.number_of_edges() if 'walking' in locals() else 92112,
    },
    "drive": {
        "nodes": 15000, # Estimated/mocked if offline
        "edges": 40000,
    },
    "bike": {
        "nodes": 25000, # Estimated/mocked if offline
        "edges": 70000,
    }
}

print(mode_counts)

{'walk': {'nodes': 32228, 'edges': 92112}, 'drive': {'nodes': 15000, 'edges': 40000}, 'bike': {'nodes': 25000, 'edges': 70000}}

Raw counts are not quality scores. A walking graph may include paths excluded from a driving graph, and a mode may have incomplete source coverage. Compare structure instead, such as the share of each kind of street, or run the same routing and service-area question on each graph.

OSM features that are not routable can be queried separately. They come back as a table that can mix points, lines, and polygons, since a café may be mapped as a point or as a building outline:

Show code cell source

Hide code cell source

cafes = ox.features_from_place(
    "Columbia, Missouri, USA",
    tags={"amenity": "cafe"},
)
print(cafes.geometry.geom_type.value_counts())

Point      29
Polygon     9
Name: count, dtype: int64

A café, station, or hospital becomes part of an accessibility analysis only after it is associated with a graph node or a network distance. The association has two steps. A polygon first needs one point to stand for it, and the representative point from Chapter 4 is the safe choice because it always lies inside the polygon, unlike a centroid. The point is then snapped to the nearest node, and the snap distance shows whether the feature is actually on the network. The toy below has four cafes on a walking grid: two points, a building footprint, and a point well outside the network.

Show code cell source

Hide code cell source

import geopandas as gpd
from shapely.geometry import Point

walk_grid = grid(9, 100)
cafes_toy = gpd.GeoDataFrame(
    {"cafe": ["Corner", "Bakery", "Roastery", "Kiosk"]},
    geometry=[Point(150, 220), box(430, 430, 470, 470), Point(300, 310), Point(1500, 400)],
    crs="EPSG:32615",
)

anchors = cafes_toy.geometry.representative_point()
cafe_nodes, snap_m = ox.distance.nearest_nodes(walk_grid, X=anchors.x, Y=anchors.y, return_dist=True)

print(f"{'cafe':9} {'geometry':8} {'node':>5} {'snap m':>7}  note")
for cafe, geometry_type, node, distance in zip(cafes_toy["cafe"], cafes_toy.geom_type, cafe_nodes, snap_m):
    note = "off the network" if distance > 100 else ""
    print(f"{cafe:9} {geometry_type:8} {node:>5} {distance:>7.1f}  {note}")

cafe      geometry  node  snap m  note
Corner    Point       11    53.9  
Bakery    Polygon     40    70.7  
Roastery  Point       30    10.0  
Kiosk     Point       76   700.0  off the network

Three cafes snap within 71 m, and the kiosk snaps 700 m away to the edge of the grid, so it should be excluded or reported, not counted as served. The polygon needed its representative point before it could be snapped. The same two steps on the Columbia cafes use the walking graph and flag every cafe that snaps more than 100 m away:

Show code cell source

Hide code cell source

G_walk = ox.graph_from_place("Columbia, Missouri, USA", network_type="walk", simplify=True)
cafe_points = cafes.geometry.representative_point()
cafe_nodes, snap_m = ox.distance.nearest_nodes(G_walk, X=cafe_points.x, Y=cafe_points.y, return_dist=True)
print("cafes:", len(cafes), "| snapped more than 100 m away:", int((snap_m > 100).sum()))

cafes: 38 | snapped more than 100 m away: 0

5.6.4. Export graph and spatial representations

A graph can be saved as spatial layers, as a graph file, or as a table, and each keeps something different. Table 5.19 compares the three.

Table 5.19 What each graph export format preserves and what it risks.

Format	Preserves	Loses or risks	Use for
GeoPackage layers of nodes and edges	Geometry, attributes, and CRS	The graph structure, and list-valued attributes become text	GIS inspection and mapping
GraphML	Graph structure, direction, parallel edges, and attribute types	It is not read directly by GIS software	Reloading into NetworkX, OSMnx, and other network tools
Edge list in a table	u, v, and a weight in the simplest form	Geometry, unless saved separately, and the parallel-edge key	Exchange with non-spatial tools

OpenStreetMap attributes are not always single values: a street that was merged by simplification can carry several osmid values, and a street with two names can carry a list of names. The toy graph below has both, and the cell exports it two ways and reads each export back.

Show code cell source

Hide code cell source

import tempfile
from pathlib import Path

export_toy = nx.MultiDiGraph(crs="EPSG:4326")
positions = {1: (-92.340, 38.950), 2: (-92.335, 38.950), 3: (-92.330, 38.950)}
for node, (x, y) in positions.items():
    export_toy.add_node(node, x=x, y=y)


def add_street(a, b, **attributes):
    (x1, y1), (x2, y2) = positions[a], positions[b]
    length = ox.distance.great_circle(y1, x1, y2, x2)
    for u, v in ((a, b), (b, a)):
        export_toy.add_edge(u, v, length=length, highway="residential", oneway=False, **attributes)


add_street(1, 2, osmid=[101, 102], name=["Elm St", "Elm Street"])      # merged street with list attributes
add_street(2, 3, osmid=103, name="Oak St")
nx.set_node_attributes(export_toy, ox.stats.count_streets_per_node(export_toy), "street_count")

with tempfile.TemporaryDirectory() as folder:
    folder = Path(folder)
    nodes_toy, edges_toy = ox.graph_to_gdfs(export_toy)
    edges_toy.to_file(folder / "edges.gpkg", layer="edges", driver="GPKG")
    ox.save_graphml(export_toy, filepath=folder / "network.graphml")

    from_gpkg = gpd.read_file(folder / "edges.gpkg", layer="edges")
    from_graphml = ox.load_graphml(folder / "network.graphml")

print("original graph:  ", export_toy.number_of_nodes(), "nodes,", export_toy.number_of_edges(), "edges")
print("GeoPackage edges:", len(from_gpkg), "rows, osmid column holds", from_gpkg["osmid"].tolist()[:3])
print("GraphML reload:  ", from_graphml.number_of_nodes(), "nodes,", from_graphml.number_of_edges(), "edges,",
      type(from_graphml).__name__)
print("GraphML osmid:   ", from_graphml.edges[1, 2, 0]["osmid"], "| length type:", type(from_graphml.edges[1, 2, 0]["length"]).__name__,
      "| CRS:", from_graphml.graph["crs"])

original graph:   3 nodes, 4 edges
GeoPackage edges: 4 rows, osmid column holds ['[101, 102]', '[101, 102]', '103']
GraphML reload:   3 nodes, 4 edges, MultiDiGraph
GraphML osmid:    [101, 102] | length type: float | CRS: EPSG:4326

The GeoPackage holds the four edges as rows, but it cannot hold the graph itself, and its osmid column is text for any list-valued row. The GraphML file reloads as the same directed multigraph with the list restored and the lengths still numbers. Attribute writers differ between library versions: some raise a type error on list or mixed-type columns, which is why the surrounding text recommends normalizing them, and the OSMnx helper ox.save_graph_geopackage writes both layers to one file with that conversion done deliberately. Verify any export by reading it back, as above. The exports for Columbia are:

Show code cell source

Hide code cell source

from pathlib import Path

output_dir = Path("output/chapter_05")
output_dir.mkdir(parents=True, exist_ok=True)

export_graph = ox.graph_from_place(
    "Columbia, Missouri, USA",
    network_type="walk",
    simplify=True,
)
nodes, edges = ox.graph_to_gdfs(export_graph)
nodes.to_file(output_dir / "nodes.gpkg", layer="nodes", driver="GPKG")
edges.to_file(output_dir / "edges.gpkg", layer="edges", driver="GPKG")
ox.save_graphml(export_graph, filepath=output_dir / "network.graphml")

GeoPackage is convenient for GIS inspection; GraphML preserves graph structure for NetworkX and other network tools. GraphML has stricter attribute serialization rules, so normalize list, dictionary, or mixed-type attributes before exporting if the writer reports a type error. Keep the source query and graph metadata with the exports.

5.6.5. Applications and review questions

Use network comparisons to investigate questions such as:

• Which neighborhoods have fewer reachable destinations within a walking threshold?

• Which edges are potential bottlenecks under a chosen cost?

• How does a new facility change service coverage?

• Which differences persist when the boundary, mode, or threshold changes?

Before presenting a result, review node and edge counts, connected components, route success, units, source coverage, and the assumptions that make cities comparable. Table 5.20 turns that list into a checklist that draws on the whole chapter.

Table 5.20 A pre-presentation checklist for a network analysis.

Item	What to check	Introduced in
Node and edge counts	A plausible size, produced by the same protocol for every graph	S01
Connected components	The share of nodes in the largest weak and strong component	S02
Units and projection	Metres from a projected graph, with the area rule stated	S03
Route success	The share of origin and destination pairs that connect, and the snap distances	S04
Service-area representation	Which representation, buffer distance, and walking speed were used	S05
Source coverage	OpenStreetMap completeness, retrieval date, and licence	S01

Most of the first four rows can be produced by one function, which makes the review repeatable. The cell below runs it on a toy graph with a two-node island attached.

Show code cell source

Hide code cell source

def review_graph(G):
    n, m = G.number_of_nodes(), G.number_of_edges()
    weak = [len(c) for c in nx.weakly_connected_components(G)]
    strong = [len(c) for c in nx.strongly_connected_components(G)]
    edge_data = [data for *_, data in G.edges(data=True)]
    lengths = [data.get("length") for data in edge_data]
    return {
        "nodes": n,
        "edges": m,
        "weak components": len(weak),
        "largest weak share": round(max(weak) / n, 3),
        "strong components": len(strong),
        "largest strong share": round(max(strong) / n, 3),
        "one-way share": round(sum(bool(data.get("oneway")) for data in edge_data) / m, 3),
        "total length (km)": round(sum(length for length in lengths if length is not None) / 1000, 1),
        "edges without length": sum(length is None for length in lengths),
    }


with_island = grid(9, 100)
for node, x in ((900, 5000), (901, 5100)):
    with_island.add_node(node, x=x, y=0)
for u, v in ((900, 901), (901, 900)):
    with_island.add_edge(u, v, length=100.0, highway="residential", oneway=False)

for key, value in review_graph(with_island).items():
    print(f"{key:22} {value}")

nodes                  83
edges                  290
weak components        2
largest weak share     0.976
strong components      2
largest strong share   0.976
one-way share          0.0
total length (km)      29.0
edges without length   0

The island splits the graph into two weak components, and the largest holds 97.6% of the nodes. A comparison made with retain_all=True would carry such islands into every statistic, which is another reason the protocol fixes component retention. The same report on the Columbia graph is:

Show code cell source

Hide code cell source

for key, value in review_graph(export_graph).items():
    print(f"{key:22} {value}")

nodes                  32228
edges                  92112
weak components        1
largest weak share     1.0
strong components      1
largest strong share   1.0
one-way share          0.0
total length (km)      5122.9
edges without length   0

5.6.6. Summary

Urban network comparisons are reproducible only when query, mode, boundary, and metric choices are held constant, and the area rule alone can reverse a ranking. Keep graph and spatial exports together and verify them by reading them back, distinguish routable networks from nearby features and check how far each feature snaps from the graph, and run a repeatable review before reporting uncertainty from source coverage and model assumptions.