Spatial Epi Public Health GIS

Queen vs Rook Contiguity Weights in libpysal

Queen and Rook contiguity look interchangeable and are not: on the same polygons they build different neighbor sets, and that difference changes the value and significance of every autocorrelation statistic downstream. This guide, part of Spatial Weights Matrix Construction for Public Health Surveillance, shows exactly where the two rules diverge in libpysal, how to measure the divergence, and how to choose defensibly.

Problem Context & Constraints

Both rules define neighbors by shared boundary. Rook contiguity counts two units as neighbors only if they share a line segment — an edge of nonzero length. Queen contiguity is more permissive: it counts units that share an edge or merely a single vertex — a corner point. The entire difference between the two is the set of corner-touching, diagonal units, and on most real administrative geographies that set is small. It is precisely because the difference is small and easy to overlook that it causes trouble: an analyst swaps the default, the neighbor graph shifts by a few edges, and a borderline Moran’s I crosses its significance threshold without anyone noticing the cause.

The effect is largest on regular or near-gridded layers. Consider a raster-derived surveillance grid or a set of nearly-square sanitary districts. Under Rook, an interior cell has exactly four neighbors; under Queen it has eight, because the four diagonal cells that touch only at a corner now count. Doubling every interior unit’s neighbor count roughly halves each neighbor’s row-standardized weight and pulls in a wider, smoother spatial context. That reshapes the spatial lag, and the spatial lag is the raw material of the global test in Global & Local Moran’s I Implementation and the hotspot z-scores in Getis-Ord Gi* Hotspot Detection.

There is also a subtler failure the Queen rule was designed to prevent. On a perfect grid, if two same-valued high cells sit diagonally adjacent, Rook sees no connection between them at all — a genuine local cluster is invisible because the only shared geometry is a corner. Queen closes that gap. The cost is that Queen can also bridge across a corner where nothing epidemiologically flows, connecting two units that share nothing but a surveyor’s coincidence. Choosing between them is choosing which error you would rather make.

Edge Neighbors versus Vertex Neighbors of a Center Cell A three-by-three grid with the center cell labelled C. The four cells sharing an edge with C — above, below, left, and right, each labelled E — are connected to C by solid lines and are neighbors under both Rook and Queen contiguity. The four diagonal corner cells, each labelled V, touch C only at a single vertex; they are connected by dashed lines and are neighbors only under Queen contiguity. Rook contiguity gives C four neighbors, Queen gives C eight. Neighbors of the center cell C solid = shared edge (E) · dashed = shared vertex only (V) C E E E E V V V V Rook: 4 neighbors (E only) · Queen: 8 neighbors (E + V) The four V cells touch C at a single corner — the entire difference between the two rules lives in those vertex-only neighbors. On gridded layers Queen doubles interior neighbor counts, reshaping every spatial lag

Prerequisites

Fix the environment so the comparison is reproducible and the contiguity extraction is well defined:

  • python 3.11
  • geopandas 0.14.4 (pulls shapely 2.0.x, pyproj 3.6.x)
  • libpysal 4.12.1 — the contiguity builders and W object
  • esda 2.6.0 — Moran’s I for the downstream comparison
  • numpy 1.26.4

Input state assumed below:

  • A polygon GeoDataFrame of areal health units with a stable, sorted index (contiguity keys on the index, so an unsorted frame yields a differently-ordered W).
  • Valid topology. Run make_valid first: a fractured shared edge — one boundary digitized as two near-coincident segments — can make Rook miss a genuinely edge-adjacent pair, muddying the very comparison you are trying to make.
  • A numeric attribute to test (an age-standardized rate, say). Contiguity itself is topological and CRS-agnostic, but reproject first anyway so the same frame can also feed distance-based weights without a coordinate-unit mix-up.

Step-by-Step Solution

Build both graphs from the identical DataFrame, row-standardize each, and compare them on three axes: neighbor cardinality, the exact set of edges that differ, and the effect on a single global statistic computed on the same variable.

# Build Queen and Rook contiguity and compare their neighbor structure.
# Pinned: libpysal==4.12.1, esda==2.6.0, geopandas==0.14.4, numpy==1.26.4
import logging
import numpy as np
import geopandas as gpd
from libpysal.weights import Queen, Rook
from esda.moran import Moran

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("contiguity.compare")

def build_and_compare(gdf: gpd.GeoDataFrame, value_col: str, seed: int = 42):
    # Stable index => identical, reproducible W row order for both rules.
    gdf = gdf.sort_values("geoid").reset_index(drop=True)

    w_queen = Queen.from_dataframe(gdf, use_index=True)
    w_rook  = Rook.from_dataframe(gdf, use_index=True)

    # --- Cardinality comparison (BEFORE standardization, so counts are integers) ---
    card_q = np.fromiter((len(w_queen.neighbors[i]) for i in w_queen.id_order), dtype=int)
    card_r = np.fromiter((len(w_rook.neighbors[i])  for i in w_rook.id_order),  dtype=int)
    log.info("Queen: mean=%.3f min=%d max=%d islands=%d",
             card_q.mean(), card_q.min(), card_q.max(), len(w_queen.islands))
    log.info("Rook:  mean=%.3f min=%d max=%d islands=%d",
             card_r.mean(), card_r.min(), card_r.max(), len(w_rook.islands))

    # --- Which units gained neighbors under Queen? (the vertex-only edges) ---
    extra = {}
    for i in w_queen.id_order:
        gained = set(w_queen.neighbors[i]) - set(w_rook.neighbors[i])
        if gained:
            extra[i] = sorted(gained)
    log.info("units with vertex-only Queen neighbors: %d of %d", len(extra), w_queen.n)

    # --- Effect on a global statistic, same variable, same seed ---
    w_queen.transform = "R"     # row-standardize AFTER counting
    w_rook.transform = "R"
    y = gdf[value_col].to_numpy(float)
    mi_q = Moran(y, w_queen, permutations=999, seed=seed)
    mi_r = Moran(y, w_rook,  permutations=999, seed=seed)
    log.info("Moran's I  Queen=%.4f (p=%.4f)  Rook=%.4f (p=%.4f)",
             mi_q.I, mi_q.p_sim, mi_r.I, mi_r.p_sim)

    return {
        "cardinality": {"queen": card_q, "rook": card_r},
        "vertex_only_edges": extra,
        "moran": {"queen": (mi_q.I, mi_q.p_sim), "rook": (mi_r.I, mi_r.p_sim)},
    }

Three outputs matter. The cardinality summary tells you how much the two rules diverge on this geography — a Queen mean far above the Rook mean signals a near-gridded layer where the choice is consequential; nearly equal means an irregular layer where it barely matters. The vertex-only edge set (extra) is the exact list of connections Queen adds, which you can map to confirm they are real adjacencies and not slivers. The Moran’s I pair, computed on the identical variable with the identical permutation seed, isolates the rule as the only moving part, so any gap in II or its pseudo p-value is attributable to contiguity alone.

Row-standardization is applied after the integer cardinality count, because the counts are what you audit and the standardized weights are what the statistic consumes. For a unit with mim_i neighbors, each row-standardized weight is

wij∗=1mi,w^{*}_{ij} = \frac{1}{m_i},

so moving from Rook (mi=4m_i = 4 interior) to Queen (mi=8m_i = 8 interior) halves the per-neighbor weight and averages over a broader, smoother neighborhood — which is why Queen almost always yields a lower, smoother Moran’s I on gridded data even though it connects more units.

Validation & Edge Cases

1. Cardinality histograms diverge only at the diagonals. Plot or log the bincount of both cardinality arrays. On a clean grid you expect the Rook histogram to concentrate at 4 (interior) with 3 (edge) and 2 (corner) tails, and the Queen histogram to concentrate at 8/5/3. If the Rook histogram shows unexpected 3s where you expect 4s, a fractured boundary is dropping a true edge neighbor — a topology bug, not a contiguity subtlety.

INFO Queen: mean=7.812 min=3 max=8 islands=0
INFO Rook:  mean=3.904 min=2 max=4 islands=0
INFO units with vertex-only Queen neighbors: 361 of 400

2. The Queen-minus-Rook edge diff should be all corners. Every entry in extra should correspond to a pair of units whose boundaries intersect in a point, not a line. Validate it directly rather than trusting it:

# Confirm each Queen-only edge is a true vertex (point) touch, not a missed edge.
# Pinned: geopandas==0.14.4, shapely==2.0.4
from shapely import intersection

def audit_vertex_edges(gdf, extra):
    suspicious = []
    for i, js in extra.items():
        gi = gdf.geometry.iloc[i]
        for j in js:
            shared = intersection(gi, gdf.geometry.iloc[j])
            # a genuine Queen-only edge meets in a point / multipoint, length 0
            if getattr(shared, "length", 0) > 0:
                suspicious.append((i, j))   # shares a line -> Rook should have had it
    return suspicious   # non-empty => topology problem upstream, not a real Queen edge

A non-empty suspicious list means a shared edge was classified as Queen-only, which is impossible under correct topology — it points to invalid geometry or a projection mismatch, and you must fix that before trusting either graph.

3. The statistic flips across the significance threshold. The consequential case is a Moran’s I that is significant under one rule and not the other:

INFO Moran's I  Queen=0.183 (p=0.071)  Rook=0.246 (p=0.038)

Here Rook declares spatial structure at α=0.05\alpha = 0.05 and Queen does not. Neither is “wrong” — they answer different questions. The resolution is to pick the rule before seeing the result, on the merits of the geometry, and record that choice, so the decision cannot be reverse-engineered from the p-value you preferred.

Compliance Notes

The contiguity rule is a modeling decision that must survive review, so log it as one:

  • Record the chosen rule and the libpysal version in the run’s provenance, exactly as the parent guide hashes the full weights spec. “Queen, row-standardized, libpysal 4.12.1” and “Rook, row-standardized, libpysal 4.12.1” are different, non-interchangeable matrices.
  • Justify the choice from geometry, not outcome. State in writing that the layer is near-gridded (Queen, to close diagonal gaps) or irregular with meaningful shared edges (either, with the reason). A choice justified only by the p-value it produced is indefensible.
  • Persist the cardinality summary and the vertex-only edge count. These are the auditable evidence that the graph you used has the neighbor structure you claim, and that its diagonal edges were validated as real touches rather than topology artifacts.