Choosing a Spatial Weights Matrix with libpysal

Pick between queen contiguity, k-nearest neighbours and distance-band weights in libpysal, handle islands and row-standardisation, and see how the choice changes spatial lag features.

The weights matrix defines what “nearby” means, and every spatial lag feature you compute inherits that definition. libpysal offers three practical families: contiguity (Queen, Rook) for polygons that share boundaries, KNN for a fixed neighbour count regardless of geometry, and DistanceBand for a fixed radius. Choosing between them is a modelling decision about the process, not a technical detail — and the wrong choice produces lag features that are subtly uninformative rather than obviously broken.

This page compares the options. For what to do with the resulting lag features, see Spatial Lag and Neighborhood Statistics.

Three definitions of “neighbour” over the same geometry The same nine units are shown three times. Under queen contiguity the centre unit links to all eight surrounding units. Under k nearest neighbours with k of four it links only to its four closest centroids. Under a fixed distance band it links to everything inside a circle, which leaves an outlying unit with no neighbours at all. The same geometry, three neighbour graphs, three different lag features Queen contiguity 8 neighbours · shared edge or corner neighbour count varies by unit KNN, k = 4 exactly 4, always · no islands possible asymmetric: A may pick B, not vice versa DistanceBand island fixed radius · sparse areas get nothing An island’s lag is 0 — which the model reads as a real, extremely low value.

Why This Fails in Geospatial ML Pipelines

Islands are the quiet killer. libpysal emits a warning about disconnected observations and then carries on; the lag for those rows is 0, and zero is a perfectly ordinary value for a standardised variable. So an island parcel — a real island, a unit at the edge of the study area, or a polygon whose neighbour was filtered out upstream — is presented to the model as a place where every neighbouring value happens to be zero. In a dataset with a few percent islands the effect is a diffuse loss of accuracy that no diagnostic points at.

Distance bands fail the same way but worse, because a single radius has to work across a study area whose density varies by orders of magnitude. Set the radius for a city and every rural unit is an island. Set it for the countryside and urban units acquire hundreds of neighbours, which makes the lag an average over half a city and destroys the local signal it was supposed to carry. min_threshold_distance gives the smallest radius that leaves no islands — often far larger than is useful.

The third failure is scale drift between training and inference. A KNN graph built over the training rows is not the graph you get when you build one over an inference batch: k nearest among 50,000 parcels is a different neighbourhood from k nearest among the 300 parcels in today’s tile. Lag features must be computed over a fixed reference geometry — the full parcel layer — not over whatever subset happens to be in the batch, which is the same contract problem described in scaling features consistently between training and inference.

Core Principles

  • Match the graph to the process. Contiguity for administrative or parcel effects, KNN for point observations, distance band only where a physical range genuinely exists.
  • Row-standardise for lag features. w.transform = "r" makes the lag a weighted mean rather than a sum.
  • Never accept islands silently. Count them, and either switch to KNN or mark those rows as missing.
  • Build the graph once, over the reference geometry. Inference reads it; it does not rebuild it.
  • Prefer k in the 4–8 range for KNN. Fewer is noisy, more washes out the local signal.
  • Persist the graph. A .gal/.gwt file, or the sparse matrix, is a versioned artifact.

Production-Ready Code

from __future__ import annotations

import logging
import geopandas as gpd
import numpy as np
from libpysal.weights import KNN, Queen, DistanceBand, W

logger = logging.getLogger(__name__)


def build_weights(gdf: gpd.GeoDataFrame, kind: str = "queen", k: int = 6,
                  threshold: float | None = None, row_standardise: bool = True) -> W:
    """Build a spatial weights matrix with explicit island handling.

    Args:
        gdf: Layer in a PROJECTED CRS — distance-based schemes need metres.
        kind: "queen", "rook", "knn" or "distance".
        k: Neighbour count for KNN.
        threshold: Radius in CRS units for the distance band.

    Raises:
        ValueError: on a geographic CRS, or when a contiguity graph leaves islands.
    """
    if gdf.crs is None or not gdf.crs.is_projected:
        raise ValueError("weights need a projected CRS so distances are in metres")

    if kind == "queen":
        w = Queen.from_dataframe(gdf, use_index=True, silence_warnings=True)
    elif kind == "rook":
        from libpysal.weights import Rook
        w = Rook.from_dataframe(gdf, use_index=True, silence_warnings=True)
    elif kind == "knn":
        w = KNN.from_dataframe(gdf, k=k)
    elif kind == "distance":
        if threshold is None:
            from libpysal.weights import min_threshold_distance
            pts = np.c_[gdf.geometry.centroid.x, gdf.geometry.centroid.y]
            threshold = float(min_threshold_distance(pts))
            logger.info("no threshold given; using island-free minimum %.1f m", threshold)
        w = DistanceBand.from_dataframe(gdf, threshold=threshold, silence_warnings=True)
    else:
        raise ValueError(f"unknown weights kind {kind!r}")

    if w.islands:
        msg = (f"{len(w.islands)} of {w.n} observations have no neighbours under "
               f"'{kind}'. Their spatial lag would be 0, which the model reads as a "
               "real low value. Use kind='knn', widen the threshold, or mask them.")
        if kind in ("queen", "rook"):
            raise ValueError(msg)
        logger.warning(msg)

    if row_standardise:
        w.transform = "r"
    counts = np.array([len(w.neighbors[i]) for i in w.id_order])
    logger.info("%s weights: n=%d, mean neighbours %.1f, min %d, max %d",
                kind, w.n, counts.mean(), counts.min(), counts.max())
    return w


def spatial_lag(w: W, values: np.ndarray) -> np.ndarray:
    """Weighted average of each unit's neighbours. NaN where a unit has none."""
    from libpysal.weights import lag_spatial
    lag = lag_spatial(w, np.asarray(values, dtype="float64"))
    if w.islands:
        island_pos = [w.id_order.index(i) for i in w.islands]
        lag[island_pos] = np.nan          # missing, not zero
    return lag

Step-by-Step Walkthrough

Step 1 — Compare the candidate graphs before committing

import geopandas as gpd

parcels = gpd.read_file("parcels.gpkg").to_crs(32633)

for kind, kwargs in [("queen", {}), ("knn", {"k": 6}), ("distance", {"threshold": 2000})]:
    try:
        w = build_weights(parcels, kind=kind, **kwargs)
    except ValueError as exc:
        print(f"{kind:>9}: rejected — {exc}")

Run this before writing any feature code. The neighbour-count summary tells you immediately whether a distance band is viable for your density range.

Step 2 — Compute lag features and check they are not degenerate

w = build_weights(parcels, kind="knn", k=6)

parcels["ndvi_lag"] = spatial_lag(w, parcels["ndvi_mean"].to_numpy())
parcels["ndvi_anomaly"] = parcels["ndvi_mean"] - parcels["ndvi_lag"]

corr = parcels[["ndvi_mean", "ndvi_lag"]].corr().iloc[0, 1]
print(f"corr(value, lag) = {corr:.3f}")
assert 0.05 < abs(corr) < 0.999, "lag is either uninformative or a copy of the value"

A correlation near 1 means the graph is too tight — often k=1, where the lag is just the nearest neighbour’s value. Near 0 means the graph is too loose to capture local structure.

Step 3 — Persist the graph as an artifact

from libpysal.io import open as psopen

with psopen("artifacts/parcels_knn6.gal", "w") as fh:
    fh.write(w)

Serving reads this file. Rebuilding the graph from an inference batch produces a different neighbourhood for every batch size, which is the lag-feature equivalent of re-fitting a scaler.

Step 4 — Handle islands deliberately at inference

lag = spatial_lag(w, new_values)
missing = np.isnan(lag)
if missing.any():
    logger.warning("%d rows have no neighbours; lag features left as NaN", missing.sum())

Leaving them NaN lets a gradient-boosting model route them down its missing-value branch, which is a far better representation than a fabricated zero.

A zero lag is indistinguishable from a genuinely low neighbourhood A distribution of spatial lag values is drawn, centred well above zero. A spike at exactly zero marks the island units, which sits inside the left tail of the real distribution, so the model treats those rows as unusually low rather than as missing. Distribution of a row-standardised lag feature 0.00.51.0 real neighbourhood means islands forced to 0 The spike is not an outlier the model can learn around — it is a lie about the data. Set them to NaN and let the estimator route them explicitly.

Verification

import numpy as np
from libpysal.weights import lat2W


def test_row_standardised_lag_is_a_weighted_mean():
    """On a 3x3 lattice with all values equal, every lag must equal that value."""
    w = lat2W(3, 3)
    w.transform = "r"
    values = np.full(9, 4.2)
    assert np.allclose(spatial_lag(w, values), 4.2)


def test_islands_become_nan_not_zero():
    w = lat2W(3, 3)
    w.transform = "r"
    w.islands.append(0)                  # simulate a disconnected unit
    lag = spatial_lag(w, np.arange(9, dtype="float64"))
    assert np.isnan(lag[0])


def test_knn_never_produces_islands():
    gdf = make_random_points(n=500, crs=32633)
    w = build_weights(gdf, kind="knn", k=5)
    assert not w.islands


def test_contiguity_with_islands_raises():
    gdf = polygons_with_one_detached_island()
    try:
        build_weights(gdf, kind="queen")
    except ValueError as exc:
        assert "no neighbours" in str(exc)
    else:
        raise AssertionError("expected the island check to fire")

Beyond unit tests, compute the same lag feature under two graphs and correlate them. If queen and KNN lags correlate above about 0.9, the choice does not matter for this dataset and you can pick on robustness grounds. If they correlate weakly, the choice is a real modelling decision and deserves to be evaluated inside the same folds as the model — see spatial cross-validation strategies.

Choosing the scheme from the geometry and the process A decision tree starts by asking whether the units tile the space. If they are contiguous polygons, queen contiguity is chosen. If they are scattered points, k nearest neighbours is chosen. If a physical interaction range is known, a distance band is chosen. All three branches end at the same island check. What are the units? and what is the process? tiling polygons → Queen contiguity scattered points → KNN, k = 4–8 known physical range → DistanceBand then: count islands and decide what they mean

FAQ

Queen or rook contiguity for polygon data?

Queen, in almost every case. Rook requires a shared edge and therefore drops corner-touching neighbours, which on grid-like cadastres yields a sparser, more brittle graph. Rook is preferable only when the process really does move along shared boundaries, such as flow between adjacent river reaches.

What should I do about islands?

Decide explicitly. An island’s lag of zero reads as an extreme low value, not as missing. Either use KNN, where every unit has neighbours by construction, or keep contiguity and set those rows to NaN as the code above does.

Should weights be row-standardised?

For lag features, yes — w.transform = "r" makes the lag a weighted average on the same scale as the variable, independent of neighbour count. Keep binary weights only when the neighbourhood total is what you actually want.

Part of: Spatial Lag and Neighborhood Statistics Part of: Spatial Feature Engineering for Machine Learning