Spatial Weights Matrix Construction for Public Health Surveillance

This guide sits within Spatial Epidemiology Fundamentals & Data Standards and covers the one object every areal cluster method silently depends on: the spatial weights matrix $W$ that encodes which health units are neighbors and how strongly they interact. Get $W$ wrong — the wrong contiguity rule, an unprojected distance band, an unhandled island — and every Moran’s I, Getis-Ord Gi*, and spatial regression coefficient downstream inherits the same defect, usually without any error being raised.

A weights matrix is not a preprocessing detail you configure once and forget. It is a modeling decision with the same evidentiary weight as the choice of statistical test, because it defines the null hypothesis of spatial randomness against which every result is judged. Two analysts who build Queen contiguity versus a 4-nearest-neighbor graph on the same county case rates are answering different questions, and their hotspot maps will differ accordingly. This page treats $W$ as versioned, hashed, reviewable infrastructure — constructed deterministically, validated for connectivity, and bound to every artifact it produces.

Concept & Epidemiological Alignment

Spatial statistics on areal health data — case counts or rates aggregated to counties, census tracts, or ZIP Code Tabulation Areas — are all functions of a weights object $W$, an $n \times n$ matrix whose entry $w_{ij}$ quantifies the spatial relationship between unit $i$ and unit $j$. A positive $w_{ij}$ says “unit $j$ is a neighbor of unit $i$ and its value should influence $i$’s local statistic”; a zero says “these units are not spatially connected.” The diagonal is always zero — a unit is never its own neighbor. Everything a spatial method claims about clustering, autocorrelation, or spillover is a statement about covariance structured through $W$.

The epidemiological question dictates the definition of “neighbor.” Contiguity weights encode the idea that disease processes diffuse across shared administrative boundaries — appropriate when transmission, environmental exposure, or care-seeking crosses adjacent units. Distance-band and kernel weights encode a continuous decay of influence with separation — appropriate when the exposure is a plume, a distance-to-service gradient, or a vector’s flight range rather than a boundary-mediated process. k-nearest-neighbor weights guarantee every unit has exactly $k$ connections regardless of how sparse or dense the geography is — appropriate when unit sizes vary wildly (rural mega-counties beside dense urban tracts) and you need a uniform local sample size. Choosing the family is choosing a hypothesis about how the health signal propagates.

That choice then propagates into every downstream statistic. The same $W$ feeds the permutation inference of Global & Local Moran’s I Implementation, the signed hotspot z-scores of Getis-Ord Gi* Hotspot Detection, and the spatially lagged design matrix of Spatial Lag and Error Regression. Because these methods share the object, they must share the same instance of it: constructing a Queen graph for the global test and a distance band for the local decomposition means the two results are no longer comparable, and a reviewer cannot reconcile them. Build $W$ once, validate it, hash it, and pass that exact object to every consumer.

Weights Family Selection

There is no universally correct weights family; there is a family matched to the geometry of the units and the biology of the exposure. The table below is the decision most areal surveillance work reduces to.

Contiguity is the reflex choice for polygon health data because it needs no distance parameter and respects the boundaries agencies already report on. Queen versus Rook is a genuine fork with measurable consequences — the vertex-only neighbors Queen adds can flip a borderline Moran’s I — and it earns its own focused treatment in Queen vs Rook Contiguity Weights in libpysal. Distance and kernel families move you into parameter territory: you must justify the threshold or bandwidth, and you must be in a projected CRS or the distances are meaningless.

Spatial Data Prerequisites

Two preconditions gate every construction below, and both fail silently rather than loudly.

First, the CRS must be projected before any distance-based weights are built. DistanceBand, KNN, and Kernel all measure separation between centroids, and a separation measured in decimal degrees is not a distance — one degree of longitude is about 111 km at the equator and about 55 km at 60° latitude, so a fixed 50 km band silently reaches different real distances across a study area, and near the poles it degenerates entirely. Enforce a documented projected or equal-area CRS as covered in Coordinate Reference Systems for Public Health before you touch a distance threshold. Contiguity weights are topological and CRS-agnostic in principle, but you should still reproject first so a mixed pipeline never applies a distance rule to unprojected coordinates by accident.

Second, topology must be valid and features stably ordered. Contiguity is derived from shared boundaries, so a duplicated vertex that fractures a shared edge into two near-coincident segments makes a contiguity rule read two genuinely adjacent counties as non-adjacent — an artificially deflated neighbor count and a missed cluster. Repair geometries with make_valid, confirm a single geometry type, and sort on a stable identifier so the row index of $W$ is reproducible. libpysal keys weights by the DataFrame index, so an unsorted input produces a differently-ordered $W$ and non-reproducible output even when the underlying graph is identical.

# Reproject + topology gate before any weights construction.
# Pinned: geopandas==0.14.4, shapely==2.0.4, pyproj==3.6.1
import logging
import geopandas as gpd
from shapely.validation import make_valid

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

def prepare_units(gdf: gpd.GeoDataFrame, projected_crs: str = "EPSG:5070") -> gpd.GeoDataFrame:
    """CONUS Albers by default; pass a local UTM/State Plane for a single region."""
    if gdf.crs is None:
        raise ValueError("Input has no CRS; refusing to guess before distance weights.")
    gdf = gdf.to_crs(projected_crs)                       # metres, not degrees
    gdf["geometry"] = gdf.geometry.apply(make_valid)      # heal fractured boundaries
    assert gdf.is_valid.all(), "Invalid geometry survived make_valid; inspect digitisation."
    assert gdf.geom_type.nunique() == 1, "Mixed geometry types; split before building W."
    gdf = gdf.sort_values("geoid").reset_index(drop=True) # stable index => reproducible W
    log.info("prepared %d units in %s", len(gdf), gdf.crs.to_authority())
    return gdf

Production Implementation: Building and Row-Standardizing W

The construction below builds each of the four families from the same prepared GeoDataFrame, inspects connectivity, and row-standardizes. Row standardization is the step that makes weights comparable across units with different neighbor counts. A raw contiguity matrix gives a corner county with two neighbors the same per-neighbor weight as an interior county with eight, so the interior county’s spatial lag sums eight unit-values and the corner county’s sums two — an artifact of geography, not epidemiology. Row-standardization rescales each row to sum to one:

$$w^{*}_{ij} = \frac{w_{ij}}{\sum_{j} w_{ij}}$$

so that a unit’s spatial lag $\sum_j w^{*}_{ij} x_j$ is a weighted average of its neighbors’ values rather than a sum. After standardization the row for a two-neighbor county assigns $0.5$ to each neighbor and the row for an eight-neighbor county assigns $0.125$ to each, and both lags live on the same scale as $x$ itself. This is the transform Moran’s I and most spatial regressions assume; skip it and the global coefficient is no longer bounded in the familiar way and cross-unit comparisons distort.

# Build all four weights families, standardize, and diagnose connectivity.
# Pinned: libpysal==4.12.1, geopandas==0.14.4, numpy==1.26.4
import logging
import numpy as np
from libpysal.weights import Queen, Rook, KNN, DistanceBand

log = logging.getLogger("weights.build")

def build_weights(gdf, *, knn_k=6, band_metres=50_000):
    """Return a dict of named, row-standardized libpysal W objects."""
    specs = {}

    # --- Contiguity: topological, no distance parameter ---
    specs["queen"] = Queen.from_dataframe(gdf, use_index=True)
    specs["rook"]  = Rook.from_dataframe(gdf, use_index=True)

    # --- k-nearest-neighbor: fixed cardinality, never produces an island ---
    specs["knn"] = KNN.from_dataframe(gdf, k=knn_k)

    # --- Distance band: requires a projected CRS (metres) ---
    specs["dband"] = DistanceBand.from_dataframe(
        gdf, threshold=band_metres, binary=True, silence_warnings=True
    )

    for name, w in specs.items():
        # Row-standardize so each row sums to 1 (transform "R").
        w.transform = "R"
        log.info(
            "W[%s]: n=%d islands=%d mean_neighbors=%.2f min=%d max=%d pct_nonzero=%.4f%%",
            name, w.n, len(w.islands), w.mean_neighbors,
            w.min_neighbors, w.max_neighbors, w.pct_nonzero,
        )
        if w.islands:
            log.warning("W[%s] has %d island(s): %s", name, len(w.islands), w.islands[:10])
    return specs

The logged pct_nonzero and neighbor-count summaries are the first sanity check. A Queen graph on U.S. counties typically averages five to six neighbors; a mean of one or a max_neighbors in the hundreds signals a topology problem (fractured boundaries deflating counts, or a sliver polygon touching everything). Always inspect the full connectivity histogram, not just the mean — a bimodal distribution usually means two datasets were concatenated in different projections.

# Connectivity histogram + island report for a constructed W.
# Pinned: libpysal==4.12.1, numpy==1.26.4
import numpy as np

def connectivity_report(w, name=""):
    card = np.fromiter((len(w.neighbors[i]) for i in w.id_order), dtype=int)
    hist = np.bincount(card)
    report = {
        "name": name,
        "n": int(w.n),
        "islands": list(w.islands),
        "n_islands": len(w.islands),
        "cardinality_histogram": {i: int(c) for i, c in enumerate(hist) if c},
        "mean_neighbors": float(card.mean()),
        "median_neighbors": int(np.median(card)),
    }
    return report

Serialization with a Provenance Hash

A weights matrix that drives a public-health decision must be reconstructable on demand, which means it is versioned and hashed like any other governed artifact. libpysal serializes to two ASCII interchange formats: GAL (adjacency lists, natural for contiguity) and GWT (weighted edge lists, natural for distance/kNN weights that carry a value). Write the file, then hash the specification — the family, parameters, library versions, and input identity — so a reviewer can prove which construction produced which map without re-deriving it.

# Serialize W to GAL/GWT and bind a provenance hash to the spec.
# Pinned: libpysal==4.12.1
import hashlib
import json
import logging
from libpysal.io import open as psopen

log = logging.getLogger("weights.provenance")

def serialize_with_provenance(w, spec: dict, gal_path: str) -> str:
    """Write W to a GAL/GWT file and return the SHA-256 of its full spec."""
    # 1. Persist the graph itself.
    with psopen(gal_path, "w") as out:
        out.write(w)

    # 2. Hash the reconstruction contract, not the floating output.
    provenance = {
        "family": spec["family"],              # queen | rook | knn | dband | kernel
        "params": spec["params"],              # {"k": 6} or {"threshold_m": 50000}
        "transform": w.transform,              # "R" for row-standardized
        "n": int(w.n),
        "n_islands": len(w.islands),
        "input_geoid_sha256": spec["input_geoid_sha256"],  # hash of sorted id list
        "libpysal_version": "4.12.1",
    }
    spec_hash = hashlib.sha256(
        json.dumps(provenance, sort_keys=True).encode()
    ).hexdigest()
    provenance["spec_sha256"] = spec_hash
    log.info("weights provenance=%s", json.dumps(provenance))
    return spec_hash

Store spec_sha256 beside the hotspot map it produced. When a Moran’s I result is challenged months later, the hash is the evidence that the challenged run used a specific, row-standardized Queen graph built from a specific set of county boundaries — not a subtly different graph an analyst rebuilt from a re-downloaded shapefile.

Parameter Selection & Tuning

Contiguity weights have no continuous parameter, only the Queen/Rook rule, which is why they are the low-controversy default. The moment you choose a distance or kNN family you own a tuning decision that must be justified in writing.

• Distance-band threshold $d$. Too small and units fall out of the graph as islands; too large and the matrix densifies toward all-pairs, washing out local structure and blowing up memory. Anchor $d$ to the exposure’s interaction range (a known vector flight distance, a plume radius, a commuting shed), then confirm it against the connectivity report: the smallest defensible $d$ that leaves zero islands is usually the right one. A common diagnostic is to compute the maximum nearest-neighbor distance across all units — any threshold below it guarantees at least one island.
• kNN parameter $k$. $k$ is the local sample size for every unit’s statistic. Small $k$ (3–4) makes local statistics noisy and sensitive to a single neighbor; large $k$ over-smooths and can bridge genuinely separate regions. For county-scale surveillance $k$ between 5 and 8 is typical. Because kNN fixes cardinality, it never produces an island — which is exactly why it is the standard tool for repairing one (see below).
• Kernel bandwidth. A fixed bandwidth treats dense and sparse regions identically; an adaptive bandwidth (set per unit from its $k$-th nearest neighbor) keeps the effective neighbor count stable across uneven density. Adaptive is almost always right for national health data where urban and frontier units coexist.

Whatever the parameter, tune it once, record it in the provenance spec, and never let it drift between the global test and the local decomposition.

Edge Cases & Failure Modes

• Islands / zero-neighbor units. A unit with no neighbors has an empty row, so row-standardization divides by zero and its local statistic is undefined. libpysal surfaces these as w.islands and warns on construction. They are not a nuisance to suppress — they are a decision point. The full detection-and-repair workflow (kNN augmentation, union of contiguity and kNN) is covered in Handling Island Polygons in Spatial Weights.
• Disconnected components. Even with no islands, $W$ can split into two or more components with no edge between them (a mainland and an archipelago each internally connected). Global statistics computed across a disconnected graph mix incomparable populations. Detect components with a connected-components pass over the adjacency and decide deliberately whether to analyze them jointly or separately.
• Self-neighbors and duplicate geometries. Coincident or overlapping polygons (a common artifact of merging multi-source boundary files) can make a unit its own neighbor or inflate a neighbor’s weight. Deduplicate geometries before construction; assert a zero diagonal after.
• Memory at $N > 50\text{k}$. A dense representation of $W$ is $O(n^2)$ and becomes infeasible past tens of thousands of units. Contiguity and kNN graphs are sparse by construction — keep them in libpysal’s sparse form and never materialize the dense matrix. A distance band with a large threshold is the trap: it can densify to near all-pairs and exhaust memory. Cap the threshold, or switch to kNN, which bounds nonzeros at $nk$ regardless of geography.
• Row-standardizing a binary distance band changes its meaning. Row standardization is standard for contiguity and kNN, but for some kernel and distance applications the raw weight is the signal; standardizing it discards the decay magnitude. Decide the transform deliberately and record it in the spec.

Compliance & Audit Controls

The weights matrix is governed infrastructure, and the audit trail for it is small but non-negotiable:

• Version the weights specification, not just the file. The GAL/GWT file is an output; the reproducible contract is the family, parameters, transform, and the hash of the sorted input identifiers. Commit the spec to version control and attach spec_sha256 to every artifact the weights produce.
• Log island reconnections explicitly. If any island was repaired, the audit record must state which units were reconnected and by what rule, because a manual reconnection is an analyst decision that changed the neighbor graph and therefore the results.
• Pin and record libpysal and geopandas versions. Contiguity extraction and index handling have changed across libpysal releases; a matrix is only reproducible against the version that built it.
• Bind the transform to the artifact. “Queen, row-standardized” and “Queen, binary” are different matrices that yield different statistics. Record w.transform in the provenance so a reviewer never has to guess.

Implementation Checklist

• All layers reprojected to a documented projected/equal-area CRS before any distance or kNN weights
• Geometries repaired with make_valid; single geometry type asserted; duplicates removed
• Features sorted on a stable identifier so the row order of $W$ is reproducible
• Weights family chosen from the exposure model and unit geometry, and justified in writing
• Connectivity histogram inspected; mean/min/max neighbor counts sane for the geography
• w.islands checked; any reconnection logged with the reconnected units and rule
• Distance threshold or k recorded in the provenance spec
• Row-standardization (transform = "R") applied deliberately and recorded
• Weights serialized to GAL/GWT and the spec hashed (SHA-256)
• The same W object passed to every downstream statistic, not rebuilt per method
• libpysal and geopandas versions pinned and recorded alongside the spec

Related Topics

• Queen vs Rook Contiguity Weights in libpysal — how the vertex-neighbor rule changes neighbor cardinality and Moran’s I.
• Handling Island Polygons in Spatial Weights — detecting and repairing zero-neighbor units before they break standardization.
• Global & Local Moran’s I Implementation — the autocorrelation test that consumes the row-standardized weights directly.
• Getis-Ord Gi* Hotspot Detection — hotspot z-scores built on the same weights object.
• Coordinate Reference Systems for Public Health — the projection discipline that distance and kNN weights depend on.