City2Graph: Urban Data as Heterogeneous Graphs

Chapter 1: Types, From Zero

Chapter 0 ended with a promise that a graph can hold relations. This chapter is about the moment that promise starts to strain, which happens sooner than you expect — on your second data source.

Here is the scenario, concretely. You have built a graph of tessellation cells around buildings, connected when they touch. It works. Now you load the street network and want the cells to know which street they front onto, because “which street” is how a shop gets customers. You have two obvious moves.

Move one: pretend everything is the same kind of thing

Dump the street segments into the same node set as the cells. Now you have one node list of, say, 240 cells plus 310 segments = 550 nodes, and one edge list mixing cell–cell adjacency with cell–street frontage.

This is a homogeneous graph: one node type, one edge type. Every message-passing rule in Chapter 8 will work on it immediately. And it is broken in two independent ways.

The first break is about features. A cell has a footprint area in square metres, a land-use share, a building height. A street segment has a carriageway width, a speed limit, a road class. To put both in one node feature matrix you must give them the same number of columns in the same order — so column 0 is “area in m² if you are a cell, speed limit in km/h if you are a street.” The model averages over neighbours; it will average 340 m² with 30 km/h and produce 185 of nothing.

The second break is about rules. In a homogeneous graph there is exactly one weight matrix W applied to every neighbour. So the rule the model learns for “how much do I take from a cell that touches me” is forced to be the same rule as “how much do I take from the street I face.” Those are different physical relations with different meanings and different reliability. You have not simplified the model; you have constrained it to be wrong.

Move two: keep two graphs

Build the cell graph and the street graph separately, run two models, concatenate the embeddings. This preserves the features but severs the thing you wanted: information cannot flow from a street to a cell during learning, because there is no edge to carry it. You have two monologues where you wanted a conversation.

The canonical edge type is a triple

Here is the piece of notation you must be fluent in before anything else makes sense. An edge type is written as an ordered triple:

For example ("place", "faced_to", "movement"). Read it aloud as “a place is faced_to a movement.” Three separate facts live in that one tuple, and forgetting any of them causes a specific bug:

Why the source and target type when the relation name already implies them? Because is_nearby is a perfectly sensible relation between many pairs of types — a shop near a stop, a stop near a segment, a cell near a cell. The relation name alone is ambiguous; the triple never is. This is why libraries in this space, including PyTorch Geometric and city2graph, key everything by the full triple.

Counting the machinery a type system costs you

Types are not free, and knowing the price is how you decide when to pay it. Suppose a graph with T node types and R edge types, where a homogeneous layer would need one weight matrix W of size d × d.

A relational layer of the RGCN kind needs one such matrix per edge type, plus a self-transform per node type. With d = 64, T = 3 and R = 5:

Eight times the parameters for one layer — which sounds alarming until you notice the absolute number is 32,768, roughly the size of a rounding error in any modern model, and that the alternative was a model structurally unable to distinguish a wall from a road. Types are cheap in parameters and expensive in bookkeeping. The bookkeeping is what the library removes.

There is a real cost, though, and it is not parameters — it is data per type. Each edge type gets its own weight matrix, so each edge type needs enough edges to estimate it. A relation with 40 edges in the whole city will learn noise. When we get to the honest-limits chapter, “how many edges does my rarest relation have?” will be the first diagnostic on the list.

How this is stored: node-keyed and edge-keyed dictionaries

Concretely, a heterogeneous graph is two dictionaries. In city2graph, before any tensor exists, they are dictionaries of GeoDataFrames — the table-plus-geometry object from GeoPandas:

python
# nodes: one GeoDataFrame per node type, keyed by type name
nodes = {
    "place":    place_gdf,      # 240 rows, index = place_id, geometry = polygon
    "movement": movement_gdf,   # 310 rows, index = movement_id, geometry = linestring
}

# edges: one GeoDataFrame per edge TRIPLE, keyed by the tuple
edges = {
    ("place",    "touched_to",   "place"):    p2p_gdf,   # MultiIndex (place_id, place_id)
    ("movement", "connected_to", "movement"): m2m_gdf,   # MultiIndex (movement_id, movement_id)
    ("place",    "faced_to",     "movement"): p2m_gdf,   # MultiIndex (place_id, movement_id)
}

Three details in that snippet are load-bearing and are the source of most first-day errors.

One: node identity is the index, not a column. city2graph takes node identifiers from the GeoDataFrame index. If your ids are sitting in a column called place_id, you must set_index("place_id") first or the library will use the default 0..n−1 integer index and your edges will point at the wrong rows — silently, because integers are valid ids.

Two: edges carry a MultiIndex. The edge GeoDataFrame’s index has two levels: source id then target id. That is how a row knows which pair it describes. Everything else in the row — length, travel time, angle — becomes a potential edge feature.

Three: the same relation name can appear in several triples and they remain distinct objects, because the key is the whole tuple.

Undirected is a lie that has to be told carefully

“Cell A touches cell B” is symmetric; you would not write it twice. But message passing needs a direction: a message goes from somewhere to somewhere. So an undirected relation stored as one row per pair must be expanded into two directed edges before any GNN sees it.

city2graph makes this explicit rather than magical. The gdf_to_pyg converter takes directed=False by default, and in that mode every edge (u, v) is symmetrized by also adding (v, u), with edge attributes duplicated and self-loops not duplicated. If you had already stored both directions yourself, that is an error the library refuses rather than silently doubling your graph — it rejects already-bidirectional inputs in undirected mode because they cannot be safely deduplicated on the way back out.

Cross-type relations are subtler. If a place faces a movement, the reverse fact is “a movement is faced by a place” — and that reverse edge cannot live in the same triple, because its source type is now movement. So it needs a triple of its own. The library’s reverse_edge_types="auto" default generates ("movement", "rev_faced_to", "place") for you. Set it to None and you get strict mode, where any undirected cross-type edge raises instead. We will count the exact tensors this produces in Chapter 7.

A worked type census

Let us fix a running example that the next several chapters extend. A small study area in a city, with:

and edge types:

Four edge types, three node types. Note immediately what the census tells you: the rarest relation is (stop, is_nearby, movement) at 92 rows, and it is the one whose weight matrix will be estimated from the least evidence. Note also that place and movement are dense enough that a mistake in either will dominate the loss. Reading a type census this way, before training anything, catches more problems than any amount of hyperparameter search.

Chapter 2: Morphology Graphs

Morphology is the study of urban form: the shape and arrangement of the built fabric. It is the oldest of the four graph families in this lesson and the one with the most surprising construction, so we build it slowly.

Start from the naive design and let it fail. Nodes are buildings, edges connect buildings that touch. In a dense terraced street this works — row houses share walls, so the graph is a chain. Move one street over to detached suburban housing and every building touches nothing. Your graph has 300 nodes and 0 edges. The representation collapsed because it was measuring walls, and the thing that actually connects a suburb is the space between the buildings.

The move: tessellate the space, not the objects

The standard fix in urban morphometrics is to stop making buildings the nodes and instead partition the ground into cells, one per building, where each cell is the region of space closer to that building than to any other. This is a morphological tessellation, and it is a Voronoi diagram taken with respect to building footprints rather than points.

Think about what you gain in one sentence: the tessellation cells always tile the plane, so they always touch their neighbours, whether or not the buildings do. The detached suburb now has a connected graph, and the cell area is itself a meaningful feature — it is roughly “how much room this building has”, the quantity a planner would call plot size.

There is a refinement that matters enormously in practice. A raw tessellation would let a cell spill across a motorway, because Euclidean distance does not know about barriers. So the tessellation is enclosed: the street network (and any other barrier you supply) is used first to cut the study area into enclosures, and cells are generated within each enclosure independently. Two houses either side of a dual carriageway now get cells that stop at the road, which is the correct answer — they are neighbours on a map and strangers in life.

python
# the tessellation step is available on its own, and it is what
# morphological_graph() calls internally
from city2graph.utils import create_tessellation

cells = create_tessellation(buildings_gdf)     # barrier-free
# with barriers, enclosures are derived from the street geometry first

Place and movement: two node types, and why

Now the second design decision, and it is the one that makes the whole graph heterogeneous. Urban morphology has two complementary things in it: the space you occupy and the space you move through. city2graph names them directly:

And three relations between them, which is exactly what morphological_graph(buildings_gdf, segments_gdf) returns:

Three relations, and each one is a genuinely different physical channel. Two cells that touch share a wall or a fence — that channel carries similarity of building age, of plot size, of who your neighbour is. Two segments that connect carry traffic. A cell facing a street is how a building gets access, footfall, noise and deliveries. A model that had to use one weight matrix for all three would be learning their average, which describes nothing.

The function returns them as the two dictionaries from Chapter 1:

python
import city2graph as c2g

nodes, edges = c2g.morphological_graph(
    buildings_gdf, segments_gdf, center_point, distance=500
)

nodes["place"]                                        # tessellation cells
nodes["movement"]                                     # street segments
edges[("place", "touched_to", "place")]           # cell adjacency
edges[("movement", "connected_to", "movement")]  # street topology
edges[("place", "faced_to", "movement")]        # frontage interface

Worked example: queen against rook on a 3 × 3 block

The contiguity rule is a parameter (contiguity="queen" by default, or "rook"), and people pick one without thinking. Compute the difference by hand once and you will never pick carelessly again.

Take nine tessellation cells arranged in a perfect 3 × 3 grid, numbered

Rook contiguity connects cells that share an edge. Count horizontally: each row has 2 adjacent pairs (1–2, 2–3), and there are 3 rows, giving 2 × 3 = 6. Count vertically: each column has 2 adjacent pairs, and there are 3 columns, giving another 6. Total

Queen contiguity also connects cells that share only a corner. Every 2 × 2 block of cells contributes exactly 2 diagonal pairs, and a 3 × 3 grid contains 2 × 2 = 4 such blocks:

Check both with the handshake identity — the sum of all degrees must be twice the edge count. Under rook, the 4 corner cells have degree 2, the 4 edge-midpoint cells degree 3, and the centre degree 4:

Under queen the corners have degree 3, the edge midpoints 5, and the centre 8:

So queen gives you 67% more edges (20 vs 12) on the same geometry, and the extra ones are exactly the corner-touchers. Which is correct? It depends on what the edge is supposed to mean. If “touched_to” stands for shared wall, shared drainage, shared fence disputes, rook is the honest rule. If it stands for you can see each other and your children play together, queen is. And if your cells came from a tessellation of real building footprints, exact corner touches are numerically fragile — a coordinate rounded in the seventh decimal place flips a queen edge on or off. That fragility is a real argument for rook on messy data.

The distance parameter is a network distance, and that matters

Real study areas are cut out of a bigger city, and how you cut matters. The distance argument to morphological_graph is not a radius on a map. It is a shortest-path distance along the street network from your center_point: segments beyond that path distance are removed, and tessellation cells are kept only if their own distance via those segments is within the budget.

This distinction is worth an example. A house 200 m away as the crow flies, on the far side of a railway cutting with the nearest bridge 900 m north, is at a network distance of roughly 200 + 2 × 900 = 2,000 m by the time you walk there and back along the cutting. A Euclidean 500 m filter includes it; a network 500 m filter correctly excludes it. If your question involves people walking, only one of those answers is about the real world.

Two related parameters deserve a sentence each because their names do not fully give them away:

Reading the returned cells

Two output options change what you can do downstream. keep_buildings=True preserves building information in the tessellation output, so a cell knows which building it wraps — you need this whenever your node features are building attributes (height, age, use). keep_segments=True (the default) keeps the original LineString of each street segment in a column named segment_geometry on the movement nodes, which you need for drawing and for length-based features, because the node geometry itself may have been simplified.

There is also a plural sibling, morphological_graphs, for the common case of building the same graph at several radii — 400 m, 800 m, 1,600 m — in one shared pass instead of re-tessellating three times. Tessellation is the expensive step, so this is not a small saving.

What morphology gives a model that features do not

Finish this chapter with the payoff, because it is easy to build a morphology graph and forget why. Suppose you want to predict which cells contain retail. Table features (area, height, distance to centre) get you some way. What the graph adds is the ability for the model to learn a pattern of arrangement: retail sits on cells that face a well-connected segment, are small, and have small neighbours — a high-street signature that is invisible in any single row and obvious two hops out.

That is the concrete meaning of “the structure is the signal.” The cell’s own features say what it is. Its neighbours’ features say what kind of place it is in. Only the second one distinguishes a small building on a high street from an identical small building behind a warehouse.

Chapter 3: Network Graphs

Streets and transit are already networks, which sounds like it should make this the easy chapter. It is the chapter with the most traps, because “already a network” hides a question nobody asks out loud: which network? There are two, they are duals of each other, and choosing wrong quietly ruins your analysis.

Primal and dual: the same streets, two graphs

In the primal graph, junctions are nodes and street segments are edges. This is the road map as you would draw it. It is the right object when you want to route: shortest path from A to B, travel time, isochrones. Lengths live on edges, which is where routing algorithms expect them.

In the dual graph, the roles swap: segments become nodes and two segment-nodes are joined when the segments met at a junction in the primal. This is the right object when the street itself is the thing you are predicting about — how busy is this street, what land use is on it, how central is it in the network. Street-level attributes (width, class, speed limit) are node features here, not edge features, which is exactly what a GNN wants.

Notice that Chapter 2 already committed to the dual: the movement nodes are street segments, and connected_to means “these two segments meet.” That was not an accident. A morphology graph is about form, and form is a property of streets, not of the abstract points where they cross.

city2graph gives you both directions explicitly:

python
import city2graph as c2g

# LineStrings -> primal graph (junction nodes + segment edges)
junction_nodes, segment_edges = c2g.segments_to_graph(segments_gdf)

# primal -> dual (segment nodes + adjacency edges)
dual_nodes, dual_edges = c2g.dual_graph((junction_nodes, segment_edges))

In the dual, a node’s geometry is the centroid of the original segment and an edge’s geometry is a LineString joining two such centroids. Pass keep_original_geom=True and the original LineString is preserved in an original_geometry column, which you will want for mapping.

Worked example: how many dual edges does a junction create?

Here is an arithmetic trap that catches people the first time they build a dual graph and wonder why it is so much denser than the primal.

In the primal, a junction where k segments meet is a single node of degree k. In the dual, that same junction becomes a clique: every one of those k segments is now adjacent to every other, giving

For an ordinary four-way crossroads, k = 4 and that is 4 × 3 / 2 = 6 dual edges from one primal node. For a five-way junction, 5 × 4 / 2 = 10. The growth is quadratic in junction degree, so a city with a handful of complex intersections gets a dual graph noticeably denser than a naive count would predict. Budget for it before you build a country-scale dual and run out of memory.

Now the modelling consequence, which is more interesting than the memory one. That clique says all four arms of the crossroads are mutually adjacent — including a pair you cannot legally turn between. The dual graph encodes meeting, not permitted movement. If turn restrictions matter to your question, they have to be added as edge attributes or by deleting edges; nothing about the construction knows about them.

Transit: a GTFS feed is a timetable, not a graph

Public transport arrives as GTFS — the General Transit Feed Specification — which is a zip of CSV files: stops.txt, routes.txt, trips.txt, stop_times.txt, calendar.txt, sometimes frequencies.txt. It is a schedule. It has no edges in it at all. What it has is, for each trip, an ordered list of stops with arrival and departure times, and separately a calendar saying which days that trip runs.

So building a transit graph means aggregating a timetable into relations, and every choice in that aggregation is a modelling choice. city2graph loads the feed into DuckDB and then summarizes it:

python
import city2graph as c2g

gtfs = c2g.load_gtfs("itm_london_gtfs.zip")   # -> in-memory DuckDB connection

nodes, edges = c2g.travel_summary_graph(
    gtfs,
    calendar_start="20250601",
    calendar_end="20250601",
    start_time="07:00:00",
    end_time="10:00:00",
)

load_gtfs reads every .txt in the archive as a table, registers a time-parsing function, and materializes stops.geometry as points from stop_lon/stop_lat when coordinates are present. travel_summary_graph then produces stop nodes and stop-to-stop edges with two attributes that carry the whole meaning of the graph:

Worked example: what “service-weighted average” actually computes

Suppose stops X and Y are consecutive on two different routes on the same Tuesday. Route 12 is a frequent bus that runs the leg 40 times and takes on average 180 seconds. Route 45 is an express that runs it 10 times and, because of a different stopping pattern, takes 300 seconds.

The naive average of the two routes is (180 + 300)/2 = 240 seconds. The service-weighted average asks instead: if you sampled a random traversal, how long would it take?

and the edge’s frequency is 40 + 10 = 50 traversals. The 36-second gap between 240 and 204 is not rounding; it is the difference between averaging routes and averaging service, and over a whole network it systematically biases every travel-time estimate towards the timetable of the rarest service. Weighting by frequency is the correct default and it is what the library does.

Frequency-based services, and the arithmetic that expands them

Some feeds do not list every trip. Instead frequencies.txt says “this pattern runs every 10 minutes from 07:00 to 08:00”, which is one row standing for many trips. With use_frequencies=True (the default), headway-based services are expanded into effective trip counts. The documented arithmetic is exactly what you would do by hand:

and if your calendar window covers five active service days, that same row contributes 5 × 6 = 30 traversals to frequency. Turn the flag off and that row contributes 1, understating the service by a factor of 30. If your betweenness centrality ranking looks bizarre, this flag is the first thing to check.

Directed or not, and why transit is not symmetric

travel_summary_graph takes directed=False by default, merging opposite-direction edges. That is usually right for accessibility work. But be aware of what it hides: one-way loops, services that run inbound in the morning and outbound in the evening, and legs where the return trip takes materially longer because it climbs a hill or crosses the busier side of a gyratory. If your question is “can people get to work and home again”, keep directed=True and check both directions.

Shared mobility, briefly

The other transportation input is GBFS — the General Bikeshare Feed Specification — loaded with c2g.load_gbfs(path), which produces DuckDB tables with station or vehicle point geometry. GBFS is a snapshot of a live system rather than a schedule, so it does not aggregate into a timetable graph the way GTFS does; you use it as a point layer, and connect it to everything else with the proximity graphs of Chapter 4.

Putting a stop into the morphology graph

This is where the type system starts to pay off. You now have transit stops (points) and street segments (lines) from two entirely different sources with no shared identifiers. There is no join key. What there is, is geography — and geography is enough:

python
# add a stop layer to the morphology graph by spatial proximity
layers = {"place": nodes["place"],
          "movement": nodes["movement"],
          "stop": stop_nodes}

bridge_n, bridge_e = c2g.bridge_nodes(layers, proximity_method="knn", k=1)
# -> new edge types named ("stop", "is_nearby", "movement"), etc.

bridge_nodes builds directed proximity edges between every ordered pair of node layers, generating a new is_nearby relation for each. You can restrict which layers act as sources or targets with source_node_types and target_node_types, which matters because the default “every ordered pair” grows quadratically in the number of layers: three layers give 3 × 2 = 6 directed layer pairs, five layers give 5 × 4 = 20.

Chapter 4: Proximity Graphs

Morphology gave us edges from touching. Networks gave us edges from meeting. This chapter handles the case that covers most urban datasets: you have a pile of points and no relations at all. Shops, schools, bus stops, bike docks, air-quality sensors, crime reports. Nothing in the file says which of them are connected, because in reality none of them are — you have to invent the edges.

That word is doing real work. When you build a proximity graph you are not discovering structure, you are asserting a hypothesis: things this close influence each other. Different rules encode different hypotheses, they give wildly different graphs on the same points, and the choice is a modelling decision you should be able to defend. So let us build every rule by hand on the same five points.

The five points we will use for everything

Coordinates in metres, in a projected coordinate system (about which there is a warning at the end of this chapter that will save you a week):

All ten pairwise Euclidean distances, computed the only way there is — Pythagoras. Two worked in full so the rest are checkable:

Rule 1: k-nearest neighbours — each node gets exactly k friends

knn_graph(gdf, k=5) connects every node to its k closest others. Take k = 2 and read the table row by row:

Five nodes × 2 = 10 directed edges, guaranteed, always. Now collapse them to unordered pairs and something instructive appears:

Four of those are reciprocated — {A,B}, {B,C}, {C,E}, {D,E} appear in both directions. Two are one-way: A lists D but D does not list A (D prefers B and E); D lists B but B does not list D (B prefers A and C). Check the arithmetic closes:

The virtue of KNN is that degree is bounded and uniform: every node has exactly k outgoing edges, so dense downtowns and empty suburbs get equal representation. The vice is the same fact seen from the other side: in a sparse area, the “nearest” five may be kilometres off and meaningless, and KNN will connect them anyway with total confidence.

Rule 2: fixed radius — a real distance, and a quadratic dial

fixed_radius_graph(gdf, radius=...) connects every pair within a threshold. This is often called a Gilbert graph. It fixes the KNN vice: an isolated point stays isolated, which is the honest answer.

With radius 50 m on our five points, the surviving pairs are those under 50: A–B (31.62), B–C (31.62), B–D (41.23), C–E (46.10). Four edges. Push the radius to 55 and you add A–D (53.85), B–E (53.15) and D–E (50.25) — seven edges, nearly double, for a 10% change in the parameter.

That sensitivity is not an accident, and the general law is worth committing to memory. For points at density λ per square kilometre, the expected number of neighbours within radius r is the density times the area of the disc:

Take a moderately busy commercial district at λ = 400 points per km² and r = 300 m = 0.3 km:

Double the radius to 600 m and the disc area quadruples:

With 2,000 points that is 2,000 × 452 ≈ 904,000 directed edges from a single parameter nudge. So: k is a linear dial and r is a quadratic one. If you are going to sweep a parameter, sweep the one whose blow-up you can predict.

Rule 3 and 4: Gabriel and relative-neighbourhood — parameter-free rules

Both KNN and fixed radius need a number you chose. The geometric graphs need nothing, because they define adjacency by a test rather than a threshold. Two of them are in the library and both are one line to state:

Work the RNG test once on our points, on the pair A–C (60.00 m). Is there a point closer to A than 60 and closer to C than 60? B is at 31.62 from A and 31.62 from C — both under 60. So A–C is not an RNG edge: B lies in the lens between them, and the relation A–C is better explained as A–B–C. Now test A–B (31.62): we need a point under 31.62 from A and under 31.62 from B, and the smallest other distances are A–D at 53.85 and B–C at 31.62 — nothing qualifies, so A–B is an RNG edge.

The family nests, which is a useful fact to keep: the minimum spanning tree is contained in the RNG, which is contained in the Gabriel graph, which is contained in the Delaunay triangulation. city2graph ships all four — euclidean_minimum_spanning_tree, relative_neighborhood_graph, gabriel_graph, delaunay_graph — so you can walk up that ladder and watch the edge count grow while never once picking a threshold.

Rule 5: Waxman — edges as coin flips

The last rule is probabilistic, and it exists because real interaction does not have a hard cutoff. People mostly shop nearby and occasionally travel far. The Waxman model gives every pair a connection probability that decays exponentially with distance:

where β (between 0 and 1) scales the overall likelihood and r₀ sets how fast it decays — larger r₀ means long-range links stay plausible.

Compute it by hand with β = 0.5 and r₀ = 100 m. For A–B at 31.62 m:

and for the far pair A–E at 83.22 m:

The nearest pair is only 1.67 times more likely to connect than the farthest one, which tells you immediately that r₀ = 100 is large relative to this point set — the decay is barely biting. Sum the probability over all ten pairs and you get the expected edge count:

About three edges out of ten possible pairs — and a different three every time you run it unless you pass seed. That is the defining property of a Waxman graph and the reason it is a modelling tool rather than a data-derived one: it is a random draw from a generative model of interaction, useful for null hypotheses (“is my real network more clustered than distance alone predicts?”) and for simulation, not for representing observed structure.

Contiguity: the polygon version of proximity

Everything above assumed points. When your units are polygons — census tracts, wards, output areas, tessellation cells — the natural relation is shared boundary, and that is contiguity_graph(gdf, contiguity="queen"), using the same queen/rook distinction we computed in Chapter 2. Under the hood it uses libpysal’s spatial weights, which is the mature implementation of this in the Python geospatial stack.

Edge weights come from the distance between polygon centroids under the chosen metric, and the metric changes the geometry of the edge as well as its number:

That third option is the one that turns a toy into a tool. “These two wards are adjacent” is a weak statement when the shared boundary is a river with no bridge; “these two wards are adjacent and 1,900 m apart along the walking network” is a strong one. Every proximity function in the module accepts distance_metric="network" with a network_gdf, so network distance is available throughout, not just here.

Containment: polygons and the points inside them

The last constructor in this family answers “which zone is this shop in?” and returns a heterogeneous graph directly:

python
# zones (polygons) linked to the POIs they contain
nodes, edges = c2g.group_nodes(zones_gdf, poi_gdf)          # predicate="covered_by" by default

The default predicate is "covered_by" rather than "within", and the difference is one sentence with real consequences: covered_by includes points that lie exactly on the boundary, within excludes them. Bus stops are placed on kerbs, kerbs are digitised as boundaries, and a surprising fraction of your points will sit precisely on a polygon edge. With the stricter predicate they silently belong to no zone at all and vanish from every zonal aggregate you compute.

The projection warning, with the number attached

Every distance in this chapter assumed metres. If your GeoDataFrame is in EPSG:4326 — plain latitude and longitude, which is what nearly every downloaded file gives you — then your coordinates are degrees, and a degree is not a distance.

A degree of latitude is about 111,320 m everywhere. A degree of longitude shrinks with the cosine of latitude. At Liverpool, latitude 53.4°:

so one unit east–west is 66 km while one unit north–south is 111 km. The ratio is

Run KNN on those coordinates and every distance is stretched north–south by 68% relative to east–west. You do not get an error. You do not get an obviously silly map. You get the wrong five nearest neighbours, systematically biased along one axis, on every node, and nothing downstream will ever tell you. Reproject to a metric CRS — for Britain that is EPSG:27700, and UTM zones or a local national grid elsewhere — before you compute a single distance.

Chapter 5: Mobility Graphs

The three families so far built edges from geometry — things touch, things meet, things are close. This one is different and, for many questions, better: the edges are observed. Somebody counted how many people went from here to there, and that count is the edge.

This is origin–destination data, and it is everywhere once you look: census migration tables (how many people moved from each district to each other district), commuting matrices from travel surveys, bike-share trip records, mobile-phone-derived flow estimates, transit smart-card taps. The 2021 England and Wales census, used in the library’s own tutorials, publishes migration flows between all middle-layer super output areas — a genuinely large OD matrix.

Two input shapes for the same information

OD data arrives in one of two forms and od_matrix_to_graph accepts both via matrix_type:

Adjacency (matrix_type="adjacency"): a square table whose rows are origins and columns are destinations. A square pandas DataFrame indexed by zone id, or a square NumPy array whose ordering matches zones_gdf.

Edge list (matrix_type="edgelist", the default): one row per flow, with an origin column, a destination column, and one or more numeric weight columns.

The adjacency form is convenient to read and quadratic to store: n zones need n² cells whether or not there is any flow. For the 7,264 middle-layer areas of England and Wales that is 7,264² = 52,765,696 cells, the overwhelming majority of them zero. The edge list only stores what happened, which is why real published OD data is almost always distributed that way.

Worked example: a four-zone commuting matrix, end to end

Four zones, daily commuting counts. Rows are origins, columns destinations:

Start by counting what is there. Four zones give 4 × 4 = 16 cells, of which the 4 diagonal entries are self-flows, leaving 4 × 3 = 12 possible directed edges. The diagonal is dropped by default (include_self_loops=False) because “people who live and work in zone 3” is usually a node attribute, not an edge — a self-loop in a GNN just adds the node to its own neighbourhood, which most architectures already do.

Step 1: directed, with a threshold

Set threshold=20. In directed mode the rule is keep the flow when weight ≥ threshold. Walk the twelve values: 120 ✓, 45 ✓, 8 ✗, 95 ✓, 210 ✓, 30 ✓, 40 ✓, 180 ✓, 12 ✗, 5 ✗, 25 ✓, 15 ✗. Four fail, so

Step 2: undirected, where the order of operations changes the answer

Now set directed=False. The documented rule is precise and worth reading twice: for each unordered pair the weight becomes the sum of both directions, and when a threshold is provided in undirected mode it is applied after this summation. So sum first:

Five undirected edges out of six possible pairs. Now look hard at the last row, because it contains the entire lesson of this chapter. The pair {Z3, Z4} has 12 one way and 15 the other — both below the threshold of 20, so in directed mode both were deleted and Z4 lost a connection entirely. In undirected mode they sum to 27 and the edge survives.

Step 3: the call

python
import city2graph as c2g

nodes, edges = c2g.od_matrix_to_graph(
    od_df, zones_gdf,
    matrix_type="edgelist",
    source_col="origin",
    target_col="destination",
    weight_cols=["flow"],
    zone_id_col="zone_id",
    threshold=20,
    directed=False,
)

The nodes GeoDataFrame comes back indexed by the zone identifier; the edges GeoDataFrame carries a MultiIndex on (source_id, target_id) — the same contract as every other constructor in the library, which is the point of having a library.

Where the edge geometry comes from, and why it is a polite fiction

With compute_edge_geometry=True (the default) each edge gets a LineString drawn between the two zone centroids. This is what produces those striking migration-flow maps — a spray of straight lines across a country.

It is worth being clear-eyed that nobody travelled along that line. It is not a route; it is a visual representation of a pair. The centroid itself is already an abstraction — the centre of a crescent-shaped ward may lie in the park it wraps around, or in the sea. For plotting, fine. For any computation involving the length of that line, not fine: if you want travel distance between zones, compute it on the network from Chapter 3, do not measure the pretty straight line.

Multiple weights, and the column you must name

Real OD tables carry several counts at once: total flow, plus a breakdown by mode or purpose. Pass them all in weight_cols and they are all preserved as edge attributes, which is exactly what you want for a GNN — edge features are a first-class input.

But thresholding needs a single number to compare, so when there is more than one weight column you must designate threshold_col. Forget it and the library cannot know whether “at least 20” means twenty commuters in total or twenty cyclists. The general rule that this special case illustrates: when several columns could serve a role, the API makes you name one rather than guessing, and every time it makes you name one it has saved you a silent wrong answer.

Reading a mobility graph: what its structure means

Once built, an OD graph answers questions no proximity graph can. Two examples, both computable directly from the object:

Weighted degree is functional size. Sum the flow on all edges into a zone and you have the number of people who arrive there daily — a direct measure of how much of a destination a place is, independent of its population or its area. Employment centres light up; dormitory suburbs do not.

Asymmetry is function. Compare in-flow to out-flow. A zone with 5,000 in and 400 out is a job centre. Reverse those and it is a commuter suburb. Roughly equal, and it is mixed-use. That single ratio, which is trivial to compute on a directed OD graph and impossible to see in any table of zone attributes, is one of the cleanest functional classifications in urban analysis — and it is exactly the kind of signal the land-use models in Chapter 9 exploit.

Chapter 6: Metapaths

You now have four graph families and, if you have been stacking them, one heterogeneous graph with several node types. Here is the question that motivates this chapter, and it is the question urban analysis actually asks:

“Which neighbourhoods are connected to each other by public transport?”

Look at your graph and notice there is no relation that answers it. You have (area, is_nearby, stop) and (stop, connects, stop). Nowhere is there an (area, ..., area) edge meaning “reachable by transit.” The relation you want is not stored; it is implied, three hops deep, and it is a composition of relations you do have:

A named, typed sequence of edge types like that is a metapath. Materializing it — walking every such path and writing down which areas end up joined — produces a genuinely new edge type, and that is what add_metapaths does.

Why not just use a deeper GNN?

This is the right objection and it deserves a real answer, because a three-layer message-passing network does propagate information three hops.

The difference is that a three-layer GNN propagates along every three-hop path, indiscriminately. It mixes area → stop → stop → area with area → stop → area → stop and with area → segment → segment → area, all in the same hidden state, all through the same activation. The typed path you cared about is in there somewhere, diluted by all the others.

A metapath instead names the composition and gives it its own edge type, which means:

That last row is not a footnote. Stacking layers to reach three hops is exactly the recipe for over-smoothing, where every node’s representation converges towards the same average. A metapath gets you the long-range relation in a one-layer model, because the distance was collapsed at construction time.

Composition is matrix multiplication — let us do it by hand

The mechanism is completely concrete. Represent each typed relation as a binary matrix and multiply. The library’s own description is exactly this: the operation multiplies typed adjacency tables to connect terminal node pairs, and can aggregate numeric edge attributes along the way.

Toy example. Three areas a1, a2, a3 and three transit stops s1, s2, s3.

A is the (area, near, stop) matrix, 3 × 3, with rows = areas and columns = stops:

S is the (stop, connects, stop) matrix. The transit line runs s1 – s2 – s3, so s1 and s2 connect, s2 and s3 connect, s1 and s3 do not:

The metapath matrix is M = A · S · A^T. Do it in two steps.

Step one: A · S — “which stops can I reach in one transit hop from my area?”

Each row of A picks out rows of S and adds them.

Read a3’s result aloud: because a3 sits by two stops on the line, one transit hop from a3 reaches all three stops. Already this is a fact that no single stored relation contained.

Step two: multiply by A^T — “which areas sit by those stops?”

A^T is A with rows and columns swapped, so its rows are stops and its columns are areas:

So the composed relation is

Reading the answer, entry by entry

Every number in M is a path count, not a yes/no. Verify a few by walking them:

The diagonal entry of 2 is the classic surprise. a3 is connected to itself by two distinct transit paths, because it borders two stops on the same line. That is a real fact — it is a measure of how well-served a3 is — but as an edge it is a self-loop, and self-loops are normally dropped when materializing a metapath because a GNN already includes a node in its own update.

Off the diagonal there are 6 non-zero entries, and because S is symmetric and A appears on both sides, M is symmetric too — so those 6 directed entries describe 3 undirected relations: {a1,a2}, {a1,a3}, {a2,a3}. All three areas are now mutually connected by transit, from an input in which no area was connected to any other area at all.

The real API

python
import city2graph as c2g

nodes, edges = c2g.add_metapaths(
    (nodes, edges),
    sequence=[
        ("area", "is_nearby", "stop"),
        ("stop", "connects",  "stop"),
        ("stop", "is_nearby", "area"),
    ],
    new_relation_name="transit_accessible",
    edge_attr="travel_time_sec",
    edge_attr_agg="sum",
)
# edges now contains ("area", "transit_accessible", "area")

Four things about that call are worth stating explicitly.

The sequence must chain. Each step’s target type has to be the next step’s source type, and the path must contain at least two steps. A mismatch is a specification error, not a silent empty result — which is the behaviour you want, because an empty metapath looks exactly like a metapath over data with no connections.

Name your relation. If new_relation_name is left as None the edges are named metapath_0. Two metapaths later you will have metapath_0 and no memory of which composition it was. Name it after the thing it means.

Attributes aggregate along the path. edge_attr names numeric edge columns to accumulate, with edge_attr_agg being "sum" (the default) or "mean". Summing travel_time_sec along area → stop → stop → area gives total journey time including both walking legs, which is the number a passenger experiences. Choose "mean" and you get an average per hop, which is almost never what anyone wants — but it is there because for some attributes (a road-quality score, say) an average is the sensible summary.

Direction matters. With directed=False (the default) both edge directions are accepted where available in the input graph. Set it True for genuinely one-way compositions, of which one-way streets and directional commuting are the obvious urban cases.

There is also add_metapaths_by_weight for the case where you want the composition governed by an edge weight budget rather than a fixed hop sequence — the natural fit for “compose until cumulative travel time exceeds 20 minutes.”

The cost, which is quadratic and will surprise you

Metapaths densify. That is their purpose, and it is also their danger, so budget it before you run it.

Suppose each area is near p stops, and each stop connects to q other stops, and each stop has p areas near it. A single area reaches roughly

With a modest p = 3 and q = 4 that is 3 × 4 × 3 = 36 new edges per area. Over 5,000 areas that is 180,000 metapath edges from an input of 5,000 × 3 = 15,000 nearby-edges and a few thousand transit links. A twelvefold increase.

Now be less modest. In a dense city centre an area might be near p = 10 stops, each connecting to q = 8 others:

This is the mechanism by which a metapath on a city-scale graph turns into an out-of-memory error, and note the shape of it — the blow-up is worst exactly where the city is densest, so it will not show up in your test run on a quiet suburb. Check p and q on your densest tile before you run the whole city.

Which metapaths are worth materializing?

A short list of the compositions that recur in urban work, so you have a starting vocabulary:

Notice that the first three all end where they began — area to area, place to place. That is the common shape, because the point of a metapath is usually to create a relation among the units you are actually predicting about, expressed through the infrastructure that connects them.

Chapter 7: Graph → Tensors

Everything so far has been GeoDataFrames: tables with geometry, human-readable indices, string identifiers. A neural network cannot consume any of that. It wants dense float tensors with integer indices and no strings anywhere.

The conversion is a single function, gdf_to_pyg, and it is worth understanding in detail because it is where geospatial thinking hands over to tensor thinking, and every mistake you make here is silent.

What PyTorch Geometric wants, stated plainly

A PyG graph is a small number of tensors. For a homogeneous graph, a Data object holds:

The one to stare at is edge_index, because its layout is unlike anything in a GeoDataFrame. It is 2 rows by E columns, not E rows by 2 columns. Column j is one edge: edge_index[0, j] is the source, edge_index[1, j] is the target. This transposed layout exists so that gathering all source features is a single contiguous index operation, x[edge_index[0]], which is the innermost operation of every message-passing layer ever written.

And those are positional integers, not your identifiers. Node “E01006512” becomes row 0. Node “place_00317” becomes row 1. The mapping from your ids to those positions is created during conversion and stored in metadata, which is how pyg_to_gdf can put your identifiers back afterwards.

A HeteroData object is the same idea, indexed by type: data['place'].x, data['place', 'faced_to', 'movement'].edge_index, and so on. Each node type gets its own independent 0-based numbering, which is the subtlety that trips people up — place 0 and movement 0 are different nodes and there is nothing in the tensors to remind you.

The call, with every argument that matters

python
import city2graph as c2g

data = c2g.gdf_to_pyg(
    nodes, edges,
    node_feature_cols={"place":    ["area", "height", "far"],
                       "movement": ["length", "width"]},
    node_label_cols={"place": ["land_use"]},
    edge_feature_cols={("place", "faced_to", "movement"): ["frontage_m"]},
    device="cpu",
    directed=False,
    reverse_edge_types="auto",
    keep_geom=True,
)

The function detects homogeneous versus heterogeneous from the input structure: pass single GeoDataFrames and you get Data, pass dictionaries and you get HeteroData. There is no mode flag, which is convenient right up until you accidentally pass a single frame and wonder where your types went.

Note that different node types may have different numbers of features — place has 3 here and movement has 2 — which is the concrete realization of the promise made back in Chapter 1. There is no shared column space and no padding.

Worked example: counting every number in the object

Take the running study area from Chapter 1 and convert it. Inputs:

Node tensors

At float32, four bytes each:

Edge tensors, after symmetrization

Apply the Chapter 1 sanity rule. The two same-type undirected relations double inside their own triple:

The cross-type relation keeps its 390 columns and gains a new triple ("movement", "rev_faced_to", "place") with another 390:

So the object has four edge types where your input dictionary had three, and the total edge column count is

Each stored as two int64 values:

Total graph: about 96 KiB. A whole neighbourhood, fully typed, fits in a fraction of the space one photograph of it would occupy — which is a useful thing to know when someone tells you graph methods do not scale. It is the metapaths and the dense proximity radii that cost memory, not the representation.

The four gotchas, in the order you will hit them

One: your identifiers must be the index. Said in Chapter 1, repeated here because this is where it bites. Node ids come from the GeoDataFrame index; edge endpoints come from the edge frame’s two-level MultiIndex. If either is a column instead, conversion either raises or — worse — succeeds against the default integer index.

Two: features must be numeric. A column of land-use strings cannot become a float tensor. One-hot or ordinally encode before conversion. If node_feature_cols is omitted the converter resolves feature columns itself, which is convenient and occasionally picks up a column you did not intend as a feature — being explicit costs one line and removes a class of mystery.

Three: keep_geom changes what round-trips. With keep_geom=True the original geometries are serialized into metadata and come back exactly. With False they are reconstructed from node positions, so your curved street segments return as straight lines between centroids. For training, false is leaner. For anything you will map afterwards, true.

Four: parallel edges need permission. Two bus routes between the same pair of stops are two rows with the same (source, target) index. In default mode that is rejected. Pass multigraph=True to promote the two-level index to a keyed (source, target, key) contract so parallel rows survive.

The round trip, and why it is the feature that matters

The whole design of this library rests on one property: conversions go both ways without loss.

python
# GeoDataFrames -> PyG -> back to GeoDataFrames
data = c2g.gdf_to_pyg(nodes, edges)
nodes2, edges2 = c2g.pyg_to_gdf(data)

# and to the other two graph libraries
G      = c2g.gdf_to_nx(nodes, edges)     # NetworkX: algorithms, exploration
G_rx   = c2g.nx_to_rx(G)                 # rustworkx: the same algorithms, much faster
G_pyg  = c2g.pyg_to_nx(data)             # PyG -> NetworkX

Think about what the round trip actually buys. You train a model, get a 16-dimensional embedding per node, and now need to look at it. Without round-tripping you are manually joining a tensor row index back to a zone code through whatever mapping you remembered to save. With it, you call pyg_to_gdf, attach the embedding as a column, and plot a map. The gap between “the model trained” and “I can see what it learned” is where most geospatial ML projects die, and closing it is not a convenience feature.

There is also validate_pyg(data), which returns the graph metadata and checks structural consistency — worth calling once in any pipeline you intend to run unattended.

Chapter 8: A Hetero-GNN, By Hand

The graph is built and converted. This chapter runs a model over it — slowly enough that you can do the arithmetic yourself, because the point is not to memorize an architecture but to see exactly what a type-aware layer computes and what a type-blind one would compute instead.

Message passing in three lines

Every graph neural network layer is the same three steps, applied to every node simultaneously:

In a heterogeneous layer, steps 1 and 2 happen once per edge type, each with its own weight matrix, and step 3 combines the per-relation results. That is the entire difference, and the worked example below makes its consequences visible.

The toy graph

Three place nodes and two movement nodes, each with a two-dimensional feature vector:

Edges: p1 touches p2, p2 touches p3. p1 and p2 both face m1; p3 faces m2. After symmetrization and reverse-edge generation, the messages arriving at places come along two relations: touched_to from other places, and rev_faced_to from movements.

Weights, chosen to be hand-computable rather than realistic:

Aggregation is the mean over each relation; the update sums the self term and the two relation terms, then applies a ReLU.

Node p1

Self term: W_self·[1.0, 0.0] = [1.0, 0.0].

Touched neighbours of p1: just p2. Mean = [0.0, 1.0]. Transform:

Faced-by neighbours of p1: just m1. Mean = [2.0, 0.0]. Transform:

Node p2

Self: [0.0, 1.0]. Touched neighbours: p1 and p3, so the mean is

Faced-by: m1 only, giving [1.0, 0.0] as before.

Node p3

Self: [0.5, 0.5]. Touched: p2 only → W_touch·[0.0, 1.0] = [0.0, 0.5]. Faced-by: m2 only:

Every value is non-negative so the ReLU changes nothing. Three updated place vectors: [2.0, 0.5], [1.375, 1.125], [0.5, 1.5].

Now do it type-blind, and compare

Collapse the types. One weight matrix W = W_touch for every neighbour regardless of what it is. Node p1 now has neighbours {p2, m1} in one undifferentiated set:

Compare to the typed answer of [2.0, 0.5]. Two things went wrong, and they are different failures.

The averaging is dimensionally incoherent. The mean [1.0, 0.5] blends a place vector with a movement vector as if their coordinates measured the same quantity. In the real graph, place feature 0 might be plot area in square metres and movement feature 0 carriageway width in metres. Their average is a number with no units and no meaning, and every subsequent layer builds on it.

The single weight matrix cannot separate the channels. Even with sensible features, the typed model applied a different transform to the street signal than to the neighbour signal — it downweighted the second dimension of street messages by a factor of four. The blind model has one knob for both. Whatever it learns for “how much of my neighbour’s land use should I absorb” is forcibly also what it learns for “how much of my street’s width should I absorb.”

Semantic attention: letting the model weigh the relations

Fixing one weight matrix per relation is the RGCN answer. The HAN answer — heterogeneous graph attention, which is what the city2graph case study uses — goes further: run the graph once per metapath, producing a separate embedding of each node under each relation, then learn how much each relation matters.

Suppose a node has two metapath embeddings: one from walk-based accessibility and one from multimodal accessibility. The semantic attention layer projects each, applies a nonlinearity, and dots it against a learned query vector to get one scalar per metapath. Say those scalars come out as

Softmax them:

and the final embedding is the weighted sum. With z_walk = [1.0, 0.0] and z_multi = [0.0, 1.0]:

Those two β values are the most useful numbers the whole model produces, and not because of the embedding. They are a read-out: this model, on this city, found walking accessibility 1.65 times as informative as multimodal accessibility for whatever it was asked to reconstruct. That is a sentence you can put in front of a planner. Very little in deep learning translates that directly into a domain claim.

What do you train it on when there are no labels?

Here is the practical problem with urban GNNs: nobody has labelled 240 tessellation cells with the answer. Land use is patchy, sold-price data is sparse, and “which places are similar” has no ground truth at all.

The standard escape is a graph autoencoder, which is self-supervised: train the model to reconstruct the graph itself. Encode every node to a low-dimensional vector, then ask the decoder to predict, from a pair of vectors, whether that pair is an edge. The loss is over real edges (should score high) and sampled non-edges (should score low), so the supervision is the graph structure you already have.

The case study accompanying the library does exactly this, with two models it names directly: GATGAE — a two-layer graph-attention encoder with a DistMult structure decoder for the homogeneous contiguity graph — and HANGAE — a two-layer HAN encoder with semantic attention across metapaths and DistMult per relation. Its published configuration is worth reading as a sanity check on scale: 23 input features, hidden dimension 8 with 4 attention heads (so 32 concatenated), output embedding dimension 16, dropout 0.6, learning rate 0.005, up to 5,000 epochs with early-stopping patience 100, negative sampling ratio 1.0, and 8 clusters for the downstream K-Means. The whole thing ran on an Apple M2 CPU with 16 GB of RAM, no CUDA.

DistMult, in one paragraph, because it is the decoder

The decoder scores a candidate edge (u, r, v) from the two embeddings and a learned diagonal vector per relation:

Work one with three-dimensional embeddings. Let z_u = [0.6, 0.2, 0.9], z_v = [0.5, 0.8, 0.1] and the relation vector d_r = [1.0, 0.5, 2.0]:

Because the relation vector is diagonal, DistMult is symmetric in u and v — swapping them gives the same score. For undirected urban relations (contiguity, accessibility) that is exactly the right inductive bias and it is cheap. For a genuinely asymmetric relation such as one-way commuting flow it is the wrong model, and you would need something like a translational or bilinear decoder instead. Knowing which of your relations are symmetric is therefore a modelling prerequisite, not a detail.

The pipeline, end to end

python
import torch
import city2graph as c2g
from torch_geometric.nn import HeteroConv, SAGEConv

# 1. build -> 2. compose -> 3. convert
nodes, edges = c2g.morphological_graph(buildings_gdf, segments_gdf, distance=500)
nodes, edges = c2g.add_metapaths((nodes, edges), sequence=walk_path,
                                 new_relation_name="walk_accessible")
data = c2g.gdf_to_pyg(nodes, edges, node_feature_cols=feat_cols)

# 4. one heterogeneous layer: a separate conv per edge type, summed per node
conv = HeteroConv({
    et: SAGEConv((-1, -1), 16) for et in data.edge_types
}, aggr="sum")

x_dict = conv(data.x_dict, data.edge_index_dict)
z = x_dict["place"]                       # [240, 16] embeddings

# 5. put the embeddings back on the map
nodes["place"]["emb_0"] = z[:, 0].detach().numpy()
nodes["place"].plot(column="emb_0", legend=True)

The (-1, -1) in SAGEConv is PyG’s lazy-initialization idiom: it means “infer the source and target input widths from the first batch”, which is precisely what you need when different node types have different feature counts. Without it you would have to hand-specify the input dimension for every one of the edge types.

Step 5 is the one people leave out and should not. An embedding you have not mapped is an embedding you have not checked. Plot the first few dimensions over the study area: if they look like noise, or if they perfectly trace your data-collection boundary, the model learned an artefact rather than the city.

Chapter 9: Real Uses, Honest Limits

Time to be specific about what this machinery is actually good for, and then equally specific about where it misleads. The second half of this chapter is the more valuable one.

Use 1: accessibility, computed properly

“How many people live within a 15-minute walk of a pharmacy?” is the single most common question in applied urban analysis, and the standard answer is wrong in a way you can quantify.

The standard answer draws a circle. Walking speed is about 4.8 km/h, so 15 minutes is

and a circle of radius 1.2 km has area

But you cannot walk through buildings. The ratio of network distance to straight-line distance is called circuity, and in ordinary street grids it runs around 1.3 to 1.4. Take 1.35: covering 1,200 m of pavement gets you only

The buffer claims nearly twice the catchment that exists. And that error is not uniform — it is worst exactly where circuity is worst, which is cul-de-sac suburbs, areas severed by railways, and anywhere a river runs. So a national buffer-based accessibility study systematically flatters precisely the places that are hardest to get around, which is the opposite of what the study was commissioned to find.

The graph version does not have this problem, because it walks the actual network:

python
from city2graph.utils import create_isochrone, filter_graph_by_distance

# everything reachable within 1,200 m along the network from a set of origins
reachable = filter_graph_by_distance(
    graph, center_points, threshold=1200, edge_attr="length"
)

# the reachable polygon, three rings out of one distance computation
iso = create_isochrone(
    graph,
    center_point=center_points,
    threshold=[400, 800, 1200],
    edge_attr="length",
    method="concave_hull_alpha",
)
# -> GeoDataFrame, one row per ring, columns: threshold, geometry

The two signatures, and why the keyword names are load-bearing

Those two calls look interchangeable and they are not. Written out in full, the parameters are:

python
filter_graph_by_distance(
    graph,             # 1st positional
    center_point,      # 2nd positional
    threshold,         # 3rd, no default
    edge_attr="length",
    node_id_col=None,
)

create_isochrone(
    graph=None,        # 1st positional
    nodes=None,        # 2nd, not center_point
    edges=None,        # 3rd positional
    center_point=None,
    threshold=None,    # singular, takes a list
    edge_attr=None,
    cut_edge_types=None,
    method="concave_hull_knn",
    **kwargs,          # absorbs any keyword typo
)

Read the second one twice. The second positional slot is nodes, not center_point. And the threshold parameter is threshold, singular, even when what you hand it is a list. Those two facts produce the two commonest ways to mis-call this pair, and they fail very differently.

The general lesson is bigger than this library. Any function ending in **kwargs so that it can forward per-method options has, in exchange, given up the ability to tell you that you misspelled a parameter. In that design a typo is not a syntax error. It is a silent no-op, and the complaint surfaces somewhere else entirely, phrased as though you had forgotten something.

What actually comes back

filter_graph_by_distance returns the same type it was handed. Give it a NetworkX graph and you get a node-induced subgraph; give it an edges GeoDataFrame and you get a GeoDataFrame of the surviving edges, nodes not included. Multiple origins are unioned rather than intersected: a node survives if it is within the threshold of any centre, which is what you want for “near a pharmacy” and emphatically not what you want for “near a pharmacy and a school”.

create_isochrone changes shape according to what you pass to threshold. A scalar gives a one-row GeoDataFrame holding a single geometry. A sequence gives one row per threshold with columns threshold and geometry — and it builds them from one Dijkstra run whose cutoff is the largest threshold, then slices that single distance table three times. Three rings cost one traversal, not three, so once you have paid for the 1,200 m ring the 400 and 800 are nearly free. Ask for them. An empty list, incidentally, raises ValueError: threshold sequence must not be empty.

Worked example: the library’s own three-node graph

The docstring ships an example small enough to check by hand, which is the best kind. Three nodes on a line at 0, 10 and 20 metres, two edges of length 10, one origin at Point(1, 0), threshold 12:

python
G.add_node(0, pos=(0, 0))
G.add_node(1, pos=(10, 0))
G.add_node(2, pos=(20, 0))
G.add_edge(0, 1, length=10)
G.add_edge(1, 2, length=10)

filter_graph_by_distance(G, Point(1, 0), threshold=12)  # nodes -> [0, 1]

Trace it in two steps. First the origin is snapped to the nearest node by straight-line distance, using a KD-tree over the node positions:

Then one Dijkstra from node 0, keeping every node whose path cost is ≤ 12:

which is the [0, 1] the docstring prints. Notice what that already tells you: node 2 is only 19 m from the origin in a straight line, a circuity of barely 1.05, and it is still cut. A network threshold is a threshold on paths, never on proximity.

Now nudge the origin twenty centimetres, from Point(4.9, 0) to Point(5.1, 0), and watch the snap flip:

Twenty centimetres of origin, the same graph and the same threshold, and the answer went from two thirds of the network to all of it. Snapping is not an implementation detail. It is a modelling decision you are making whether or not you noticed making it, and its blast radius is the distance between junctions.

The isochrone is conservative, and the buffer was optimistic

One more consequence hides in the phrase node-induced. The reachable set is a set of nodes whose path cost is ≤ the threshold, and an edge survives only when both of its endpoints do. Nothing is clipped mid-block. So the ring stops at the last junction inside your budget, not at the point along the street where the budget actually ran out.

Put numbers on it, with a 1,200 m budget. On a dense grid with junctions every 100 m the last reachable junction sits at exactly 1,200 m and you lose nothing. On a suburban layout with junctions every 350 m the reachable ones are at

so the ring stops at 1,050 m and 150 m of perfectly walkable pavement is discarded. That is 12.5% of the radius, and because area goes as the square of the radius the shortfall in area is

Two knobs that change what the polygon means

edge_attr names the edge column that gets summed along a path. Pass "length" and the threshold is in metres, so what you have drawn is really an iso-distance. Pass "travel_time" and the threshold is in minutes and you have a true isochrone. The function does not change; the only thing that changed is which column you told it to add up. Name it explicitly every time, because although the docstring advertises a default of "travel_time", the traversal falls back to "length" when edge_attr is None — a discrepancy you will never notice on a graph carrying both columns, and will notice very sharply on one that does not.

method decides how a cloud of reachable nodes becomes a polygon, and the four options are not stylistic:

And cut_edge_types drops whole relation types before the walk begins, which is how you carve a walk-only isochrone out of a multimodal graph: pass cut_edge_types=[("bus_stop", "is_next_to", "bus_stop")] and the bus network stops carrying anyone, so the polygon that comes back is what your legs alone can reach.

Use 2: land use and function, inferred from structure

Official land-use data is coarse, stale, and legally rather than functionally defined — a building zoned “commercial” in 1998 may now be flats. What a heterogeneous GNN offers is a functional classification derived from how a place is connected rather than how it is labelled.

This is exactly the case study published with the library: clustering urban functions in Liverpool. The setup, from the paper abstract, is worth stating in full because it is the clearest available template for this kind of work. Graph autoencoder models were run over census units with three relation types — spatial contiguity, walk-based accessibility, and multimodal accessibility — and the result reported is that compared with the homogeneous model, the heterogeneous models identified spatially coherent clusters that aligned more closely with defined accessibility patterns.

Read that claim carefully, because it is honest in a way results sections often are not. It does not say accuracy went up — there is no ground truth to be accurate against. It says the clusters were spatially coherent and matched accessibility structure. For an unsupervised method on an unlabelled problem, that is the right kind of claim, and it is the kind you should expect to be able to make about your own work.

The inputs were all open: ONS Output Area boundaries and population-weighted centroids for December 2021, Overture Maps places, land use and transportation layers from the 2025-12-17.0 release, and the UK Department for Transport’s Bus Open Data GTFS timetables for the North West.

Use 3: everything network science already knew

Not every use needs a neural network. Once the graph exists, the entire toolbox of network analysis applies for free, because gdf_to_nx hands you a NetworkX object:

And when the graph gets big enough that NetworkX crawls, nx_to_rx moves it to rustworkx for the same algorithms at a different speed, then rx_to_nx brings it back. That is the round-trip design paying off again.

Limit 1: the modifiable areal unit problem

This is the deepest limitation and it has been known in geography for a century. Results change when you change the zones, and there is no correct zoning.

Concretely: run your model on census output areas (roughly 125 households each) and you get one clustering. Run the identical model on the larger super output areas and you get a different one, with different boundaries and sometimes different conclusions. Neither is wrong. The units are administrative artefacts, drawn for enumeration convenience, and any statistic computed on them inherits their arbitrariness.

Tessellation cells are a partial escape, since they are derived from the buildings rather than from administrative history — but only partial, because the tessellation itself has parameters (enclosures, barriers, buffers) that change the cells. There is no zoning-free view of a city.

Limit 2: every edge is an assumption you made

The proximity chapter should have made you uneasy, and it should have. Setting k=5 in a KNN graph asserts that five is the number of neighbours that matter. Choosing a 500 m radius asserts a catchment. Choosing queen over rook asserts that corner-touching is meaningful.

None of those are measurements. The model then learns from the structure you asserted and produces a result that looks like a finding about the city but is partly a finding about your parameters. The mitigation is the same as for MAUP: sweep the parameter, report the sensitivity. If your conclusion holds for k in 3 to 8, say so. If it only holds at k = 5, you have discovered something about k.

Limit 3: correlation, in a domain that badly wants causation

Suppose the model learns that cells near well-connected streets have higher retail density. True, useful, and utterly silent on the question everyone in the room is about to ask: if we build a new connection here, will retail follow?

The graph cannot answer that. It is an observational snapshot of a city whose form and function co-evolved over centuries; the streets are where they are partly because of the retail. Predicting from structure is legitimate. Predicting the effect of changing the structure is a causal question requiring a causal design — a natural experiment, a difference-in-differences around an actual intervention, an instrument. The temptation to slide from one to the other in a slide deck is enormous and must be resisted, because planning decisions get made on those slides.

Limit 4: the data has politics

OpenStreetMap and Overture Maps are volunteer-derived and unevenly complete: wealthy, well-mapped districts have more detail than poor ones, and “fewer amenities in the data” is indistinguishable from “fewer amenities.” GTFS covers scheduled services, so informal transport — which is the majority of transport in much of the world — is absent, and a model trained on GTFS will conclude those neighbourhoods are unreachable. Mobile-phone mobility data under-represents anyone without a smartphone.

An accessibility score computed on incomplete data reads as a statement about a neighbourhood when it is partly a statement about who bothered to map it. If that score then informs where investment goes, the data gap becomes a funding gap. Check coverage per area before you compare areas.

Limit 5: the sparse-relation trap, and other engineering realities

The last row deserves a sentence because it silently corrupts more analyses than any of the others. If you clip your study area and then build the graph, every node near the boundary lost half its neighbours, so betweenness collapses at the edges and your “most peripheral areas” are simply the ones you cut. Build with context, then clip.

Chapter 10: Using The Library

Everything up to here was the reasoning. This chapter is the reference you will come back to: install, get data, build, convert, and the map of the whole API in one table.

Install

city2graph supports Python 3.12 to 3.14. The base install gives you graph construction and spatial network analysis with GeoPandas, NetworkX and rustworkx, and deliberately does not pull in PyTorch:

bash
pip install city2graph

Add the GNN half only when you need it. The cpu extra is the right default for development and small-scale work:

bash
pip install "city2graph[cpu]"          # PyTorch + PyTorch Geometric, CPU
pip install "city2graph[cu130]"        # with CUDA 13.0; also cu126, cu128

Two details from the installation docs that will save you a debugging hour. The supported CUDA extras are cu126, cu128 and cu130; the cpu, cu126 and cu130 extras use PyTorch 2.13 or newer, while cu128 pins to PyTorch 2.11 because PyTorch stopped publishing CUDA 12.8 wheels after that release.

conda-forge also carries the package (conda install -c conda-forge city2graph), but you must add PyTorch and PyTorch Geometric yourself, and the project’s own docs recommend pip or uv for the smoothest path when you need them — conda is no longer officially supported by PyTorch or PyG.

Given how many geospatial dependencies are involved (GeoPandas, shapely, pyproj, libpysal, momepy, DuckDB), pin your environment. The case study repository does this with uv and a committed uv.lock, plus a .python-version file pinning 3.12.8. Copy that pattern; a geospatial stack that resolves differently next month is not a reproducible result.

Getting data with no GIS work at all

The data module fetches Overture Maps directly, either for a bounding box or for a named place resolved through Nominatim:

python
import city2graph as c2g

# by place name (geocoded to a boundary polygon)
data = c2g.load_overture_data(
    place_name="Liverpool, UK",
    types=["building", "segment", "connector", "place", "land_use"],
    output_dir="data/raw/overture",
)
buildings = data["building"]
segments  = data["segment"]

# or just the administrative boundary on its own
boundary = c2g.get_boundaries("Liverpool, UK")

The available Overture types include building, segment (the LineStrings that make up the transportation network), connector (the points joining them), place (points of interest), land_use, water, infrastructure and the division layers. area and place_name are mutually exclusive — pass one.

Raw Overture segments are not immediately usable as a street network, because a single segment record can span several connectors. process_overture_segments is the cleanup step: it splits segments at their connectors and can generate the barrier geometries that the enclosed tessellation of Chapter 2 wants. Skipping it is the most common reason a first morphological graph comes out with implausibly long streets and a nearly empty connected_to relation.

One reproducibility note: pass an explicit release (for example the case study’s 2025-12-17.0). Overture keeps only the most recent monthly releases available, so a pipeline that silently takes “latest” will not rebuild the same graph in six months.

The full pipeline, one screen

python
import city2graph as c2g

# ---- 1. data ----------------------------------------------------------
ov = c2g.load_overture_data(place_name="Liverpool, UK",
                            types=["building", "segment", "connector"])
buildings = ov["building"].to_crs(27700)      # METRIC CRS. do this first.
segments  = c2g.process_overture_segments(ov["segment"]).to_crs(27700)

# ---- 2. morphology ----------------------------------------------------
nodes, edges = c2g.morphological_graph(
    buildings, segments,
    center_point=city_centre, distance=1000,
    contiguity="queen", keep_buildings=True,
)

# ---- 3. transit, then bridge it in ------------------------------------
gtfs = c2g.load_gtfs("north_west_gtfs.zip")
stop_nodes, stop_edges = c2g.travel_summary_graph(
    gtfs, calendar_start="20250603", calendar_end="20250603",
    start_time="07:00:00", end_time="10:00:00",
)
nodes["stop"] = stop_nodes.to_crs(27700)
edges[("stop", "connects", "stop")] = stop_edges.to_crs(27700)

_, near = c2g.bridge_nodes(nodes, proximity_method="knn", k=2,
                           source_node_types=["place"], target_node_types=["stop"])
edges.update(near)

# ---- 4. compose a metapath --------------------------------------------
nodes, edges = c2g.add_metapaths(
    (nodes, edges),
    sequence=[("place", "is_nearby", "stop"),
              ("stop",  "connects",  "stop"),
              ("stop",  "is_nearby", "place")],
    new_relation_name="transit_accessible",
    edge_attr="travel_time_sec", edge_attr_agg="sum",
)

# ---- 5. tensors -------------------------------------------------------
data = c2g.gdf_to_pyg(nodes, edges, node_feature_cols=feat_cols)
c2g.validate_pyg(data)

# ---- 6. train, then come back to the map ------------------------------
# ... your model ...
nodes_out, edges_out = c2g.pyg_to_gdf(data)

Step 1 contains the single most important line in the whole script, and it is the to_crs(27700). Do it before anything measures anything.

The API map

Cheat sheet: the decisions, with defaults worth keeping

Where to go next

The library’s own tutorials each cover one family from this lesson end to end on real data: morphological graphs from Overture Maps and OpenStreetMap; GTFS to public transit graphs, including betweenness ranking and walk-plus-transit isochrones; OD matrices to mobility graphs, up to the full 2021 census migration flows between all middle-layer areas of England and Wales; spatial proximity graphs over Tokyo points of interest under Euclidean, Manhattan and network distance; and metapath construction for heterogeneous GNNs. There is also a FOSS4G 2026 workshop repository, GeoAI in Practice, that walks from open data to a graph-autoencoder clustering pipeline.

On this site, the natural next steps are CS224W: Heterogeneous Graphs for RGCN, HAN and HGT in full detail, Graph Neural Networks for the message-passing foundations, and CS224W: GNN Design Space for how the aggregation and update choices in this chapter interact.

References

1. Sato, Y., Pietrostefani, E., Mahabir, R., & Arribas-Bel, D. “City2Graph: A Python library for Heterogeneous Graph Neural Networks and spatial analysis in urban systems.” Computers, Environment and Urban Systems, 130, 102492, 2026. doi:10.1016/j.compenvurbsys.2026.102492
2. city2graph source repository and README. github.com/c2g-dev/city2graph (BSD 3-Clause)
3. city2graph documentation — overview, installation, tutorials and full Python API reference. city2graph.net
4. Sato, Y. “Case Study Data for City2Graph: Clustering Urban Functions in Liverpool.” Zenodo, 2026. doi:10.5281/zenodo.18396285
5. city2graph case study repository — GATGAE, HANGAE and the K-Means baseline, with the experiment configuration quoted in Chapter 8. github.com/c2g-dev/city2graph-case-study
6. Archived software releases. doi:10.5281/zenodo.15858845
7. Overture Maps Foundation — buildings, segments, connectors, places and land use. overturemaps.org (© OpenStreetMap contributors)
8. PyTorch Geometric — Data, HeteroData and the heterogeneous convolution wrappers. pytorch-geometric.readthedocs.io