Chapter 10: Spatial Relationships and Weights¶
Part II — Geographic Data Science
Learning Objectives¶
By the end of this chapter you will be able to:
- Compute distance metrics and contiguity-based weights
- Use libpysal to build spatial neighborhood structures
- Prepare neighborhoods for simulation models
# Standard imports — add chapter-specific imports below
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
Distance and Proximity¶
Every technique so far in Part II has asked what's at a location. This chapter asks a different question: given two locations, are they neighbors? The answer turns out to depend entirely on which of two philosophies you pick.
A fixed number of neighbors (k-nearest). Every point gets exactly k neighbors — the k closest, whatever their actual distance. In a dense cluster of points those neighbors sit close by; in a sparse area, "closest" might still mean tens of kilometers away. libpysal's KNN builds exactly this:
np.random.seed(42)
n_points = 30
coords = np.random.rand(n_points, 2) * 10
points_gdf = gpd.GeoDataFrame(
geometry=gpd.points_from_xy(coords[:, 0], coords[:, 1])
)
from libpysal.weights import KNN
w_knn = KNN.from_dataframe(points_gdf, k=5)
print("Neighbors of point 0:", w_knn.neighbors[0])
Neighbors of point 0: [np.int64(12), np.int64(27), np.int64(16), np.int64(22), np.int64(3)]
A fixed distance (distance band). Every point's neighbors are whoever falls within a given radius — no fixed count, so a point in a dense cluster picks up many neighbors and an isolated point might have none at all:
from libpysal.weights import DistanceBand
w_band = DistanceBand.from_dataframe(points_gdf, threshold=2.0)
print("Neighbors of point 0:", w_band.neighbors[0])
Neighbors of point 0: [12]
/home/sergio/develop/github/lambdageo/ebooks/geospatial-modeling-python/.venv/lib/python3.12/site-packages/libpysal/weights/util.py:819: UserWarning: The weights matrix is not fully connected: There are 5 disconnected components. There is 1 island with id: 6. w = W(neighbors, weights, ids, **kwargs) /home/sergio/develop/github/lambdageo/ebooks/geospatial-modeling-python/.venv/lib/python3.12/site-packages/libpysal/weights/distance.py:893: UserWarning: The weights matrix is not fully connected: There are 5 disconnected components. There is 1 island with id: 6. W.__init__(
Both w_knn and w_band are libpysal W objects — the library's central data structure for a spatial neighborhood, and the object every technique in this chapter eventually produces one way or another. A W object knows how to draw itself directly onto a GeoDataFrame's plot, which is the fastest way to build intuition for how differently these two philosophies carve up the same 30 points:
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
points_gdf.plot(ax=axes[0], color="grey")
w_knn.plot(points_gdf, ax=axes[0], edge_kws=dict(color="orangered", linewidth=0.8))
axes[0].set_title("KNN (k=5)")
points_gdf.plot(ax=axes[1], color="grey")
w_band.plot(points_gdf, ax=axes[1], edge_kws=dict(color="orangered", linewidth=0.8))
axes[1].set_title("DistanceBand (2.0)")
plt.show()
Contiguity Weights: Queen and Rook¶
Distance-based neighbors make sense for points. For polygons — the São Luís census tracts that anchor this section, continuing Part II's Maranhão thread — a more natural question is usually contiguity: does this polygon touch that one?
Two definitions of "touch" are standard, borrowed from how the chess pieces move:
- Rook contiguity: polygons are neighbors only if they share an edge (a line segment of positive length) — not just a corner point.
- Queen contiguity: polygons are neighbors if they share an edge or even a single vertex — a strictly looser, larger definition than Rook.
Every Queen neighbor pair is automatically a Rook-or-shares-a-corner pair; the two definitions only disagree at the corners, where four polygons can meet at a single point without sharing any edge with each other.
url = (
"https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/"
"malhas_de_setores_censitarios__divisoes_intramunicipais/censo_2022/setores/shp/UF/MA_setores_CD2022.zip"
)
setores_ma = gpd.read_file(url)
sao_luis = setores_ma.query("NM_MUN == 'São Luís'").reset_index(drop=True)
print(len(sao_luis), "census tracts in São Luís")
1736 census tracts in São Luís
from libpysal.weights import Queen, Rook
w_queen = Queen.from_dataframe(sao_luis, use_index=False)
w_rook = Rook.from_dataframe(sao_luis, use_index=False)
print(f"Queen: {w_queen.mean_neighbors:.1f} average neighbors per tract")
print(f"Rook: {w_rook.mean_neighbors:.1f} average neighbors per tract")
Queen: 6.0 average neighbors per tract Rook: 5.6 average neighbors per tract
fig, axes = plt.subplots(1, 2, figsize=(12, 6))
sao_luis_web = sao_luis.to_crs(epsg=3857) # Web Mercator, for a cleaner basemap-style plot
for ax, w, title in zip(axes, [w_queen, w_rook], ["Queen", "Rook"]):
sao_luis_web.plot(ax=ax, color="whitesmoke", edgecolor="lightgray")
w.plot(
sao_luis_web,
edge_kws=dict(linewidth=0.8, color="orangered"),
node_kws=dict(marker=".", s=3),
ax=ax,
)
ax.set_axis_off()
ax.set_title(f"{title} contiguity — São Luís census tracts")
plt.show()
Queen's edges visibly outnumber Rook's — every corner-touching pair Rook misses, Queen catches. For most spatial analysis this difference is a minor technicality; for a cellular automaton in Part III, it's the entire rule set: a fire model using Rook neighbors treats a fire that only touches a corner as non-adjacent, while a Queen-neighbor model lets it spread diagonally. Chapter 18 picks this exact distinction back up.
libpysal Weights Objects¶
A W object is more than the .neighbors dictionary glimpsed above — it's a sparse representation of what would otherwise be an enormous, mostly-empty matrix (most polygons don't neighbor most other polygons). Three attributes cover what you'll reach for most often:
print("Neighbors of tract 0:", w_queen.neighbors[0])
print("Weights of tract 0: ", w_queen.weights[0]) # all 1.0 — unweighted, so far
print("Cardinality (neighbor count) of tract 0:", w_queen.cardinalities[0])
Neighbors of tract 0: [1, 2, 3, 4, 6, 7] Weights of tract 0: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] Cardinality (neighbor count) of tract 0: 6
Every weight is currently 1.0 — a tract with 6 neighbors contributes 6 to any neighbor-based sum, one with 3 neighbors contributes only 3, simply because it has fewer of them. Row-standardization fixes that by dividing each tract's weights by its own neighbor count, so every tract's weights sum to exactly 1.0 regardless of how many neighbors it has — turning a neighbor sum into a neighbor average. This matters enormously for Chapter 12's spatial autocorrelation statistics, which are built directly on top of row-standardized weights:
w_queen.transform = "r" # 'r' for row-standardized
print("Weights of tract 0 after row-standardization:", w_queen.weights[0])
print("Sum of tract 0's weights:", sum(w_queen.weights[0]))
# The full, dense N×N matrix representation — useful for small W, expensive for large ones
W_matrix, ids = w_queen.full()
print("Full matrix shape:", W_matrix.shape)
Weights of tract 0 after row-standardization: [0.16666666666666666, 0.16666666666666666, 0.16666666666666666, 0.16666666666666666, 0.16666666666666666, 0.16666666666666666] Sum of tract 0's weights: 1.0 Full matrix shape: (1736, 1736)
That last line is worth pausing on: full() materializes every one of São Luís's N × N tract pairs, the overwhelming majority of them 0.0. For a few hundred tracts that's a manageable few-hundred-thousand-cell array; for tens of thousands of polygons it stops being one. The sparse .neighbors/.weights representation libpysal uses internally exists precisely to avoid ever needing that dense matrix in memory — call .full() only when you genuinely need it, such as feeding a W into a library that doesn't understand libpysal's native format.
Caching Neighborhoods for Performance¶
Computing Queen or Rook contiguity means checking, for every pair of polygons, whether their boundaries touch — work that grows roughly with the square of the polygon count. For a few hundred census tracts that finishes in a moment; for a state-wide or national tract layer, recomputing it from scratch on every run of an analysis becomes the slowest step by far, even though the underlying geometry never changed.
The fix is the same one any expensive, input-determined computation reaches for: compute once, cache the result, and only recompute when the input actually changes. libpysal W objects support this directly through pickle:
import pickle
from pathlib import Path
cache_path = Path("sao_luis_queen.pkl")
if cache_path.exists():
with open(cache_path, "rb") as f:
w_queen_cached = pickle.load(f)
print("Loaded cached weights —", w_queen_cached.n, "tracts")
else:
w_queen_cached = Queen.from_dataframe(sao_luis, use_index=False)
with open(cache_path, "wb") as f:
pickle.dump(w_queen_cached, f)
print("Computed and cached weights —", w_queen_cached.n, "tracts")
Loaded cached weights — 1736 tracts
The pattern generalizes well beyond libpysal specifically: any pipeline that derives an expensive result from a stable input benefits from keying a cache off exactly what determines that result — here, implicitly, the tract geometries themselves. A production geospatial pipeline processing many derived layers typically makes this explicit, hashing the inputs and parameters into a cache key rather than relying on a single hardcoded filename the way the cell above does for simplicity.
From Weights to Simulation Neighbors¶
Every neighborhood built in this chapter — KNN, distance band, Queen, Rook — answers the same underlying question a cellular automaton needs answered before it can take a single simulation step: for this cell, who influences it next? DisSModel's CellularAutomaton.create_neighborhood(), which Chapter 18 introduces in full, is doing exactly the computation this chapter just walked through by hand, for whichever substrate the model runs on:
- On a raster grid (Chapter 8), neighbors come from array slicing —
grid[i-1:i+2, j-1:j+2]for a Moore (queen-equivalent) neighborhood — because a regular grid's adjacency is implicit in its row/column indices. - On a vector
GeoDataFrame(Chapter 7), neighbors come from exactly thelibpysalWobjects built in this chapter — adjacency has to be computed explicitly, because nothing about a polygon's row position in a table implies anything about its neighbors in space.
The caching concern from the previous section is not incidental to this bridge — it's the whole reason create_neighborhood() exists as a separate, cacheable step rather than being recomputed inside every single simulation tick. A model that runs for a thousand time steps needs its neighborhood structure exactly once, not a thousand times.
Exercises¶
- k versus radius. Rerun the Distance and Proximity KNN example with
k=3andk=10. Does every point still end up with a symmetric relationship — if A is a KNN-neighbor of B, is B necessarily a KNN-neighbor of A? (Hint: think about a point on the edge of the cloud versus one in the dense center.) - Count the Queen-only edges. Using
w_queenandw_rook(before row-standardization), find at least one pair of tracts that are Queen neighbors but not Rook neighbors. (Hint: comparew_queen.neighbors[i]andw_rook.neighbors[i]for a few values ofi.) - Row-standardize DistanceBand. Apply
.transform = "r"tow_bandfrom the Distance and Proximity section and inspectw_band.weights[0]before and after. Does a point with many neighbors end up with smaller or larger individual weights than a point with few? - Cache invalidation. The caching example above never checks whether
sao_luisitself has changed since the cache was written. Sketch, in a few lines of pseudocode, how you would extendcache_pathto include something derived from the input data itself, so that a genuinely new dataset doesn't silently load a stale, unrelated cachedW.
# Your code here
Summary¶
Key concepts introduced¶
- Two philosophies of "neighbor": fixed-count (
KNN) and fixed-distance (DistanceBand), for point data - Contiguity for polygons: Queen (shares an edge or a vertex) versus Rook (shares an edge only)
- The
libpysalWobject —.neighbors,.weights,.cardinalities— and row-standardization (.transform = "r") as the step that turns a neighbor sum into a neighbor average .full()and why the sparseWrepresentation exists: a dense matrix grows with the square of the polygon count- Caching an expensive, geometry-derived
Wobject rather than recomputing it on every run - The direct line from every neighborhood built here to
CellularAutomaton.create_neighborhood()in Part III — the same adjacency question, asked once per model run instead of once per chapter section
Chapter 11 takes a brief detour to a different kind of spatial data — individual point locations, before any aggregation into polygons at all — before Chapter 12 puts these weights to their first real analytical use: measuring whether a variable's values cluster in space more than chance alone would predict.
Further Reading¶
- libpysal documentation, Spatial Weights: https://pysal.org/libpysal/notebooks/weights.html
- Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer Academic — the classical reference for contiguity-based spatial weights
- libpysal documentation, KNN and DistanceBand weights: https://pysal.org/libpysal/generated/libpysal.weights.KNN.html