[[7]{.chapter-number}  [The Spatial `Graph`]{.chapter-title}]{#sec-graphs .quarto-section-identifier}

7 The Spatial `Graph`

In formal spatial models, we need a way of defining relationships between observations. In the spatial analysis and econometrics literature, the data structure used to relate observations to one another is canonically called a spatial weights matrix, or \(W\), and is used to specify the hypothesized interaction structure that may exist between spatial units (Getis, 2009). These relationships between units form a special kind of network graph embedded in geographical space, and using the libpysal.graph module, we can explore these relationships in depth.

Code

%load_ext jupyter_black
%load_ext watermark
%watermark -v -a "author: eli knaap" -d -u -p libpysal

Author: author: eli knaap

Last updated: 2025-11-23

Python implementation: CPython
Python version       : 3.12.12
IPython version      : 9.7.0

libpysal: 4.13.0

Note

The Graph class is only available in versions of libpysal >=4.9. With versions prior to that, you need to use the W class, which works a little differently than demonstrated here.

Code

import contextily as ctx
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import networkx as nx
import osmnx as ox
import pandarm as pdna

from geosnap import DataStore
from geosnap import io as gio
from libpysal.graph import Graph

datasets = DataStore()

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

The graph will use the geodataframe index to identify unique observations, so when we have something meaningful (like a FIPS code) that serves as an “id variable”, then we should set that as the index.

Code

dc = gio.get_acs(datasets, state_fips="11", years=2021)
dc = dc.to_crs(dc.estimate_utm_crs())
dc = dc.set_index("geoid")
dc.plot()
plt.show()

/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/geosnap/io/util.py:273: UserWarning: Unable to find local adjustment year for 2021. Attempting from online data
  warn(
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/geosnap/io/constructors.py:218: UserWarning: Currency columns unavailable at this resolution; not adjusting for inflation
  warn(

The most common way of creating a spatial interaction graph is by “building” one from an existing geodataframe. The libpysal.graph.Graph class has a several methods available to generate these relationships, all of which use the Graph.build_* convention.

7.1 Contiguity Graphs

A contiguity graph assumes observations have a relationship when they share an edge and/or vertex. These are conceptually simple relationships that make sense for polygon-based data, and they encode “neighborhoods” when observations share sides (edges), corners (nodes), or a combination of the two. When observation data are encoded as points, it’s technically possible to use contiguity relationships if we rely on certain assumptions, however, strictly speaking it’s not possible for points to be be contiguous with one another.

Since many datasets in urban studies are based on polygons (e.g. administrative units), and a great deal of early research in spatial analysis was based on polygon lattice data (think counties in the U.S.), contiguity relationships are one of the most common forms of spatial Graphs used in practice.

Note

To use “contiguity” relationships using point-based observation data, it is common to use a Vornoi tesselation, which results in a set of implied polygons that encode each point (based on creating equidistance boundaries). When point data are used to create contiguity Graphs, libpysal uses this technique under the hood.

7.1.1 Rook

Code

g_rook = Graph.build_contiguity(dc)
g_rook.n

The number of observations in our graph (the .n attribute) should match the number of rows in the dataframe from which it was constructed.

Code

dc.shape[0]

One way to interact with the graph structure directly is to use the adjacency attribute which returns the neighbor relationships as a multi-indexed pandas series.

Code

g_rook.adjacency.head()

focal         neighbor    
110010001011  110010001021    1
              110010001022    1
              110010001023    1
              110010041003    1
              110010055032    1
Name: weight, dtype: int64

Code

g_rook.adjacency.reset_index()

	focal	neighbor	weight
0	110010001011	110010001021	1
1	110010001011	110010001022	1
2	110010001011	110010001023	1
3	110010001011	110010041003	1
4	110010001011	110010055032	1
...	...	...	...
2943	110019800001	110010108002	1
2944	110019800001	110010108003	1
2945	110019800001	110010108004	1
2946	110019800001	110010108005	1
2947	110019800001	110010110021	1

2948 rows × 3 columns

“Slicing” into a Graph object returns a pandas series where the index is the neighbor’s id/key and the value is the “spatial weight” that encodes the strength of the connection between observations.

Code

g_rook["110019800001"].head()

neighbor
110010001023    1
110010047021    1
110010056022    1
110010056023    1
110010058011    1
Name: weight, dtype: int64

Note this is shorthand for slicing into the adjacency list itself

Code

g_rook.adjacency.loc["110019800001"].head()

neighbor
110010001023    1
110010047021    1
110010056022    1
110010056023    1
110010058011    1
Name: weight, dtype: int64

If you prefer, the neighbors and weights values are also encoded as dictionaries on the Graph, available under the corresponding attributes.

Code

g_rook.neighbors["110019800001"][:5]

('110010001023',
 '110010047021',
 '110010056022',
 '110010056023',
 '110010058011')

Code

g_rook.weights["110019800001"][:5]

(1, 1, 1, 1, 1)

The cardinalities attribute stores the number of connections/relationships for each observation. In graph terminology, that would be the number of edges for each node, and in simple terms, this is the number of neighbors for each unit of analysis

Code

g_rook.cardinalities.head()

focal
110010001011    5
110010001021    8
110010001022    5
110010001023    8
110010002011    3
Name: cardinalities, dtype: int64

Code

g_rook.cardinalities.describe()

count    571.000000
mean       5.162872
std        1.919363
min        2.000000
25%        4.000000
50%        5.000000
75%        6.000000
max       30.000000
Name: cardinalities, dtype: float64

We can get a sense for what the neighbor distribution looks like by plotting a histogram of the cardinalities

Code

g_rook.cardinalities.hist()
plt.show()

The attribute pct_nonzero stores the share of entries in the \(n \times n\) connectivity matrix that are non-zero. If every observation were connected to every other observation, this attribute would equal 1

Code

g_rook.pct_nonzero

0.9041807625421343

In this graph, less than 1% of the entries in that matrix are non-zero, showing that this is a very sparse connectivity graph indeed. Note, to compute the pct_nonzero measure ourselves, we could alternatively rely on properties of the graph:

Code

(g_rook.n_edges / g_rook.n_nodes**2) * 100

0.9041807625421342

To see the full \(n \times n\) matrix representation of the graph, the sparse (scipy) representation is available under the sparse attribute

Code

g_rook.sparse

<Compressed Sparse Row sparse array of dtype 'float64'
    with 2948 stored elements and shape (571, 571)>

Code

g_rook.sparse.todense()

array([[0., 1., 1., ..., 0., 0., 0.],
       [1., 0., 1., ..., 0., 0., 0.],
       [1., 1., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 1., 0.],
       [0., 0., 0., ..., 1., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], shape=(571, 571))

To see the dense version, just convert from sparse to dense using scipy conventions. Both rows and columns of the matrix representation are ordered as found in the geodataframe from which the Graph was constructed (or in the order given, when using a _from_* method). Currently this order is stored under the unique_ids attribute.

Code

g_rook.unique_ids

Index(['110010001011', '110010001021', '110010001022', '110010001023',
       '110010002011', '110010002012', '110010002021', '110010002022',
       '110010002023', '110010002024',
       ...
       '110010109002', '110010110011', '110010110012', '110010110013',
       '110010110021', '110010110022', '110010111001', '110010111002',
       '110010111003', '110019800001'],
      dtype='object', name='focal', length=571)

Code

g_rook.unique_ids.equals(dc.index)

True

Summing across rows of the matrix is another way to count neighbors for each observation

Code

g_rook.sparse.todense().sum(axis=1)[:20].astype(int)

g_rook.cardinalities.values[:20]

array([5, 8, 5, 8, 3, 6, 6, 6, 4, 6, 5, 3, 9, 6, 6, 7, 8, 3, 6, 7])

Classic spatial connectivity graphs are typically very sparse. That is, we generally do not consider every observation to have a relationship with every other observation. Rather, we tend to encode relationships such that observations are influenced directly by other observations that are nearby in physical space. But because nearby observations also interact with their neighbors, it is possible for an influence from one observation to influence another observation far away by traversing the graph.

One useful way to understand how the Graph encodes spatial relationships is to use the plot method, which embeds the graph in planar space showing how observations are connected to one another. In the default plot, each observation (or its centroid, if the observations are polygons) is shown as a dot (or node in graph terminology), and a line (or edge) is drawn between nodes that are connected. Since we are using a rook contiguity rule, the edge, here, indicates that the observations share a common border.

Code

ax = g_rook.plot(dc, node_kws=dict(alpha=0.4), edge_kws=dict(alpha=0.4), figsize=(5, 6))
ctx.add_basemap(ax, source=ctx.providers.CartoDB.Positron, crs=dc.crs)
ax.axis("off")
plt.show()

The explore method creates an interactive webmap version of the same plot, though the arguments are a bit different.

Code

m = dc.explore(tiles="CartoDB Positron", tooltip=["geoid"])
g_rook.explore(
    dc,
    m=m,
)

Make this Notebook Trusted to load map: File -> Trust Notebook