Spatial Networks

📖 Ahead of time…

In this block we cover some of the analytics you can obtain when you consider street networks as spatial graphs rather than as geo-tables.

  • A good example of applying concepts and ideas presented in this block is Boeing (2020) [Boe20b]

  • Boeing (2017) [Boe17] provides a general overview on the osmnx project

  • A brief overview of momepy, the package for urban morphometrics, available in Fleischmann (2019) [Fle19]

💻 Hands-on coding

import pandas
import geopandas
import momepy
import networkx as nx
import contextily
import matplotlib.pyplot as plt

Assuming you have the file locally on the path ../data/:

db = geopandas.read_file("../data/arturo_streets.gpkg")

If you’re online, you can do:

db = geopandas.read_file(
    "http://darribas.org/gds4ae/_downloads/67d5480f98453027d59bf49606a7ad92/arturo_streets.gpkg"
)

To make things easier later, we “explode” the table so it is made up of LINESTRINGS instead of MULTILINESTRINGS:

db_tab = db.explode().reset_index()
/tmp/ipykernel_312/1525410135.py:1: FutureWarning: Currently, index_parts defaults to True, but in the future, it will default to False to be consistent with Pandas. Use `index_parts=True` to keep the current behavior and True/False to silence the warning.
  db_tab = db.explode().reset_index()
../../_images/07-Spatial_networks_9_0.png
db_tab.head()
level_0 level_1 OGC_FID dm_id dist_barri average_quality population_density X Y value geometry
0 0 0 1 1 1606 3.277406 1149.981915 444133.736820 4.482809e+06 NaN LINESTRING (444096.316 4482762.870, 444171.158...
1 1 0 2 2 1606 3.113298 795.151054 444192.038205 4.482878e+06 NaN LINESTRING (444212.943 4482901.091, 444199.698...
2 2 0 3 3 1606 3.143822 478.869601 444134.537507 4.482885e+06 NaN LINESTRING (444097.968 4482915.826, 444171.158...
3 3 0 4 4 1603 5.740980 2749.417163 445612.690578 4.479336e+06 NaN LINESTRING (445608.884 4479346.815, 445616.438...
4 4 0 5 5 1603 5.468163 2653.489421 445606.319326 4.479354e+06 NaN LINESTRING (445603.738 4479362.060, 445608.884...

Analysing street geo-tables

Length

length = db_tab.to_crs(
    epsg=32630 # Expressed in metres
).geometry.length
length.head()
0    118.699481
1     62.210799
2     95.164472
3     23.503065
4     16.090295
dtype: float64
ax = db_tab.assign(
    length=length
).plot(
    "length", 
    scheme="fisherjenkssampled", 
    k=9, 
    legend=True, 
    linewidth=0.5,
    figsize=(12, 12),
    cmap="magma"
)
contextily.add_basemap(
    ax, 
    crs=db_tab.crs, 
    source=contextily.providers.CartoDB.PositronNoLabels,
    alpha=0.5
)
ax.set_title("Street segment length");
../../_images/07-Spatial_networks_14_0.png

Challenge

Create a quantile choropleth of length for the section of the network with dist_barri values starting by 01.

Bonus tip: you can create the subset of the network using the isin method:

ids = ['0101', '0102', '0103', '0104', '0105', '0106']
subnet = db_tab[db_tab['dist_barri'].isin(ids)]

Linearity

linearity = momepy.Linearity(db_tab).series
linearity.head()
0    1.000000
1    0.999999
2    1.000000
3    1.000000
4    1.000000
dtype: float64
ax = db_tab.assign(
    linearity=linearity
).plot(
    "linearity", 
    scheme="fisherjenkssampled", 
    k=9, 
    legend=True, 
    linewidth=0.5,
    figsize=(12, 12),
    cmap="magma"
)
contextily.add_basemap(
    ax, 
    crs=db_tab.crs, 
    source=contextily.providers.CartoDB.PositronNoLabels,
    alpha=0.5
)
ax.set_title("Street segment linearity");
../../_images/07-Spatial_networks_18_0.png

Challenge

Create a choropleth of linearity for the subnet table you have created above

Streets as spatial graphs

From geo-table to spatial graph:

db_graph = momepy.gdf_to_nx(db_tab)

db_graph
<networkx.classes.multigraph.MultiGraph at 0x7fd9843d4b80>

Now db_graph is a different animal than db that emphasizes connections rather than attributes.

db_graph.is_directed()
False
db_graph.is_multigraph()
True

The (first and last) coordinates of each street segment become the ID for each segment in the graph:

print(db_tab.loc[0, "geometry"])
LINESTRING (444096.3161762458 4482762.870216271, 444171.158127317 4482855.001910598)
l = db_tab.loc[0, "geometry"]
l.coords
<shapely.coords.CoordinateSequence at 0x7fd974ab4940>
node0a, node0b = edge0 = list(
    db_tab.loc[0, "geometry"].coords
)
edge0
[(444096.3161762458, 4482762.870216271),
 (444171.15812731703, 4482855.001910598)]

We can use those to extract adjacencies to each node:

db_graph[node0a]
AdjacencyView({(444171.15812731703, 4482855.001910598): {0: {'level_0': 0, 'level_1': 0, 'OGC_FID': '1', 'dm_id': '1', 'dist_barri': '1606', 'average_quality': 3.277406, 'population_density': 1149.981915, 'X': 444133.736820226, 'Y': 4482808.89166328, 'value': nan, 'geometry': <shapely.geometry.linestring.LineString object at 0x7fd97668bbe0>, 'mm_len': 118.69948078964639}}, (444083.8275243509, 4482747.422611062): {538: {'level_0': 538, 'level_1': 0, 'OGC_FID': '539', 'dm_id': '539', 'dist_barri': '1606', 'average_quality': 3.225291, 'population_density': 2833.605516, 'X': 444090.105664431, 'Y': 4482755.13506047, 'value': nan, 'geometry': <shapely.geometry.linestring.LineString object at 0x7fd985c4ad30>, 'mm_len': 19.864413729824115}}})

We can access edge information for each pair of nodes with a concatenated dict query:

db_graph[node0a][node0b]
AtlasView({0: {'level_0': 0, 'level_1': 0, 'OGC_FID': '1', 'dm_id': '1', 'dist_barri': '1606', 'average_quality': 3.277406, 'population_density': 1149.981915, 'X': 444133.736820226, 'Y': 4482808.89166328, 'value': nan, 'geometry': <shapely.geometry.linestring.LineString object at 0x7fd97668bbe0>, 'mm_len': 118.69948078964639}})
db_graph[node0a][node0b][0]
{'level_0': 0,
 'level_1': 0,
 'OGC_FID': '1',
 'dm_id': '1',
 'dist_barri': '1606',
 'average_quality': 3.277406,
 'population_density': 1149.981915,
 'X': 444133.736820226,
 'Y': 4482808.89166328,
 'value': nan,
 'geometry': <shapely.geometry.linestring.LineString at 0x7fd97668bbe0>,
 'mm_len': 118.69948078964639}
db_graph[node0a][node0b][0]["geometry"]
../../_images/07-Spatial_networks_35_0.svg

If we need all the node IDs:

list(
    db_graph.nodes
)[:5] # Limit to the first five elements
[(444096.3161762458, 4482762.870216271),
 (444171.15812731703, 4482855.001910598),
 (444212.942998509, 4482901.090971609),
 (444097.96831143444, 4482915.825653204),
 (445608.8837672261, 4479346.814511424)]

And same for edges:

list(
    db_graph.edges
)[:5] # Limit to the first five elements
[((444096.3161762458, 4482762.870216271),
  (444171.15812731703, 4482855.001910598),
  0),
 ((444096.3161762458, 4482762.870216271),
  (444083.8275243509, 4482747.422611062),
  538),
 ((444171.15812731703, 4482855.001910598),
  (444212.942998509, 4482901.090971609),
  1),
 ((444171.15812731703, 4482855.001910598),
  (444097.96831143444, 4482915.825653204),
  2),
 ((444212.942998509, 4482901.090971609),
  (444254.9705938099, 4482866.143285849),
  10886)]

Or:

db_graph.edges[node0a, node0b, 0]
{'level_0': 0,
 'level_1': 0,
 'OGC_FID': '1',
 'dm_id': '1',
 'dist_barri': '1606',
 'average_quality': 3.277406,
 'population_density': 1149.981915,
 'X': 444133.736820226,
 'Y': 4482808.89166328,
 'value': nan,
 'geometry': <shapely.geometry.linestring.LineString at 0x7fd97668bbe0>,
 'mm_len': 118.69948078964639}

If you want fast access to adjacencies:

db_graph.adj[node0a]
AdjacencyView({(444171.15812731703, 4482855.001910598): {0: {'level_0': 0, 'level_1': 0, 'OGC_FID': '1', 'dm_id': '1', 'dist_barri': '1606', 'average_quality': 3.277406, 'population_density': 1149.981915, 'X': 444133.736820226, 'Y': 4482808.89166328, 'value': nan, 'geometry': <shapely.geometry.linestring.LineString object at 0x7fd97668bbe0>, 'mm_len': 118.69948078964639}}, (444083.8275243509, 4482747.422611062): {538: {'level_0': 538, 'level_1': 0, 'OGC_FID': '539', 'dm_id': '539', 'dist_barri': '1606', 'average_quality': 3.225291, 'population_density': 2833.605516, 'X': 444090.105664431, 'Y': 4482755.13506047, 'value': nan, 'geometry': <shapely.geometry.linestring.LineString object at 0x7fd985c4ad30>, 'mm_len': 19.864413729824115}}})

Challenge

Create the graph version of subnet and consider the street segment indexed in the table as 53271. Check the adjacencies on both ends of the segment using db_graph and subnet. Are they the same in both graphs? Why?

Analysing graphs

There are many ways to extract information and descriptives from a graph. In this section we will explore a few that can tell us important information about the position of a node or edge in the network and about the broader characteristics of sections of the graph.

Degree

Degree tells us the number of neighbors of every edge, that is how many other nodes it is directly connected to.

degree = list(db_graph.degree)
degree[:5]
[((444096.3161762458, 4482762.870216271), 2),
 ((444171.15812731703, 4482855.001910598), 3),
 ((444212.942998509, 4482901.090971609), 3),
 ((444097.96831143444, 4482915.825653204), 3),
 ((445608.8837672261, 4479346.814511424), 2)]

Node centrality

Fraction of nodes a node is connected to:

nc = pandas.Series(
    nx.degree_centrality(db_graph)
)
nc.head()
444096.316176  4.482763e+06    0.00004
444171.158127  4.482855e+06    0.00006
444212.942999  4.482901e+06    0.00006
444097.968311  4.482916e+06    0.00006
445608.883767  4.479347e+06    0.00004
dtype: float64
nc.plot.hist(bins=100, figsize=(6, 3));
../../_images/07-Spatial_networks_53_0.png

Tip

Other variations of centrality measures are available in networkx. The are computationally demanding but relatively straightforward to calculate using the library. For a few of those, you can check:

Challenge

Create a histogram of degree for db_graph. Now replicate the figure for the case of subnet. What can you learn about the two graphs by doing this exercise?

Meshedness

The messedness of a graph captures the degree of node edge density as compared to that of nodes. Higher meshedness is related to denser, more inter-connected grids.

%time meshd = momepy.meshedness(db_graph, distance=500)
CPU times: user 1min 9s, sys: 353 ms, total: 1min 10s
Wall time: 1min 9s
meshd.nodes[node0a]
{'meshedness': 0.058823529411764705}
pandas.Series(
    {i: meshd.nodes[i]["meshedness"] for i in meshd.nodes}
).plot.hist(bins=100, figsize=(9, 4));
../../_images/07-Spatial_networks_60_0.png

Challenge

Replicate the computation of meshedness for sub_graph using a threshold of 250m and 500m. How do the distributions of both compare with each other?

Betweenness centrality

How often do shortest-path routes pass through a given node?

This is computationally very demanding, so we will work on a subset of the full graph:

ids = ['0101', '0102', '0103', '0104', '0105', '0106']
subnet = db_tab[db_tab['dist_barri'].isin(ids)]
sub_graph = momepy.gdf_to_nx(subnet)
node_sub = subnet.loc[53271, 'geometry'].coords[0]

Calculating it is trivial with momepy:

%%time
betweenness = momepy.betweenness_centrality(sub_graph)
CPU times: user 18.6 s, sys: 0 ns, total: 18.6 s
Wall time: 18.6 s

As with meshedness, we obtain another graph in return with the information attached to it:

betweenness.nodes[node_sub]
{'betweenness': 0.0011679130514459708}

Attaching information to street segments

The trick here is to be able to transfer back the information stored as graphs into geo-tables so we can apply everything we already now about manipulating and mapping data in that structure. With momepy, we can bring a graph back into a geo-table:

nodes = momepy.nx_to_gdf(
    meshd, points=True, lines=False
)
nodes.head()
meshedness nodeID geometry
0 0.058824 0 POINT (444096.316 4482762.870)
1 0.092308 1 POINT (444171.158 4482855.002)
2 0.101449 2 POINT (444212.943 4482901.091)
3 0.065574 3 POINT (444097.968 4482915.826)
4 0.000000 4 POINT (445608.884 4479346.815)
ax = nodes.plot(
    "meshedness", 
    scheme="fisherjenkssampled",
    markersize=0.1,
    legend=True, 
    figsize=(12, 12)
)
contextily.add_basemap(
    ax, 
    crs=nodes.crs,
    source=contextily.providers.CartoDB.DarkMatterNoLabels
)
ax.set_title("Meshedness");
../../_images/07-Spatial_networks_72_0.png

With other measures index on node IDs, we can use joining machinery in pandas:

nc.head()
444096.316176  4.482763e+06    0.00004
444171.158127  4.482855e+06    0.00006
444212.942999  4.482901e+06    0.00006
444097.968311  4.482916e+06    0.00006
445608.883767  4.479347e+06    0.00004
dtype: float64
degree_tab = pandas.DataFrame(
    degree, columns=["id", "degree"]
)
degree_tab.index = pandas.MultiIndex.from_tuples(
    degree_tab["id"]
)
degree_tab = degree_tab["degree"]
degree_tab.head()
444096.316176  4.482763e+06    2
444171.158127  4.482855e+06    3
444212.942999  4.482901e+06    3
444097.968311  4.482916e+06    3
445608.883767  4.479347e+06    2
Name: degree, dtype: int64
net_stats = pandas.DataFrame(
    {"degree": degree_tab, "centrality": nc},
)
net_stats.index.names = ["x", "y"]
net_stats.head()
degree centrality
x y
444096.316176 4.482763e+06 2 0.00004
444171.158127 4.482855e+06 3 0.00006
444212.942999 4.482901e+06 3 0.00006
444097.968311 4.482916e+06 3 0.00006
445608.883767 4.479347e+06 2 0.00004
net_stats_geo = nodes.assign(
    x=nodes.geometry.x
).assign(
    y=nodes.geometry.y
).set_index(
    ["x", "y"]
).join(net_stats)

net_stats_geo.head()
meshedness nodeID geometry degree centrality
x y
444096.316176 4.482763e+06 0.058824 0 POINT (444096.316 4482762.870) 2 0.00004
444171.158127 4.482855e+06 0.092308 1 POINT (444171.158 4482855.002) 3 0.00006
444212.942999 4.482901e+06 0.101449 2 POINT (444212.943 4482901.091) 3 0.00006
444097.968311 4.482916e+06 0.065574 3 POINT (444097.968 4482915.826) 3 0.00006
445608.883767 4.479347e+06 0.000000 4 POINT (445608.884 4479346.815) 2 0.00004
f, axs = plt.subplots(1, 2, figsize=(18, 9))
vars_to_plot = ["degree", "centrality"]
for i in range(2):
    net_stats_geo.plot(
        vars_to_plot[i], 
        scheme="fisherjenkssampled",
        markersize=0.2,
        legend=True, 
        ax=axs[i]
    )
    contextily.add_basemap(
        axs[i], 
        crs=nodes.crs,
        source=contextily.providers.CartoDB.DarkMatterNoLabels
    )
    axs[i].set_title(f"Node {vars_to_plot[i]}")
../../_images/07-Spatial_networks_78_0.png

Challenge

Create choropleths for node and betweenness centrality for sub_graph. How do they compare?

🐾 Next steps

If you found the content in this block useful, the following resources represent some suggestions on where to go next:

  • The NetworkX tutorial is a great place to get a better grasp of the data structures we use to represent (spatial) graphs

  • Parts of the block benefit from the section on urban networks in Geoff Boeing’s excellent course on Urban Data Science

  • If you are interested in urban morphometric analysis (the study of the shape of different elements making up cities), the momepy library is an excellent reference to absorb, including its user guide