Building Custom Spatial Weights Matrices in Python

TL;DR #

Build a custom libpysal.weights.W with exponential distance decay by passing neighbor and weight dictionaries to W(neighbors, weights, id_order=list(range(n))). Use scipy.spatial.cKDTree.query_ball_point to find candidates within max_dist, compute w = exp(-decay_rate * d) for each pair, exclude self-loops manually, then call w_obj.transform = "r" for row-standardization.

Why This Matters #

Standard contiguity rules — Rook, Queen, and fixed distance bands — assume that spatial interaction depends solely on shared boundaries or Euclidean proximity. Real phenomena are rarely so uniform. Hydrological networks propagate influence along flow direction, transit accessibility follows network travel time, and ecological dispersal decays with resistance surfaces rather than straight-line distance. When the weights topology misrepresents the underlying process, model residuals accumulate unmodeled spatial autocorrelation, biasing coefficient estimates and inflating Type I error rates.

Custom matrix construction, as covered in the spatial weight matrices cluster, decouples neighbor discovery from rigid default rules and lets you inject domain-specific impedance functions directly into the weight topology. This page shows the complete implementation path from a GeoDataFrame to a validated W object ready for spreg regression or ESDA clustering within the broader Core Concepts of Spatial Statistics & Geostatistics workflow.

The diagram below maps the data flow from input geometry to downstream model consumption.

Environment and Version Pinning #

pip install libpysal==4.10.0 geopandas==0.14.3 scipy==1.13.0 numpy==1.26.4

import numpy as np                          # 1.26.4
import geopandas as gpd                     # 0.14.3
from scipy.spatial import cKDTree          # scipy 1.13.0
from libpysal.weights import W             # libpysal 4.10.0

CRS requirement: all distance calculations must use a projected coordinate reference system (meters). Geographic CRS (EPSG:4326) produces degree-based distances that are dimensionally incompatible with Euclidean metrics; cKDTree has no concept of the curvature correction needed for degree inputs.

Step-by-Step Implementation #

Step 1 — Validate the GeoDataFrame and reset the index #

def prepare_geodataframe(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """
    Enforce projected CRS, valid geometries, and zero-based integer index.
    Always call this before any distance-based weight construction.
    """
    if gdf.crs is None:
        raise ValueError("GeoDataFrame has no CRS. Assign one before computing weights.")
    if gdf.crs.is_geographic:
        raise ValueError(
            f"Geographic CRS detected ({gdf.crs.to_epsg()}). "
            "Reproject to a metric CRS, e.g. gdf.to_crs('EPSG:32618')."
        )
    # Repair degenerate geometries (slivers, unclosed rings)
    if not gdf.geometry.is_valid.all():
        gdf = gdf.copy()
        gdf["geometry"] = gdf.geometry.make_valid()
    # Zero-based sequential index is mandatory for W dict-key alignment
    return gdf.reset_index(drop=True)

libpysal.weights.W keys its neighbor and weight dictionaries on integer observation indices. A non-sequential or string-based index causes silent misalignment when you later join the weights back to a regression DataFrame.

Step 2 — Extract centroids and build a cKDTree #

def extract_centroids(gdf: gpd.GeoDataFrame) -> np.ndarray:
    """Return an (n, 2) array of centroid (x, y) coordinates in CRS units."""
    return np.column_stack([
        gdf.geometry.centroid.x,
        gdf.geometry.centroid.y
    ])

centroids = extract_centroids(gdf)
tree = cKDTree(centroids)

cKDTree.query_ball_point returns all points within a given radius in $O(n \log n)$ time, versus the $O(n^2)$ cost of a naive pairwise distance matrix. For datasets above roughly 50,000 observations this difference is the boundary between a sub-second query and a multi-hour computation.

Step 3 — Compute exponential decay weights #

The weight formula used here is:

$$w_{ij} = \exp(-\lambda \cdot d_{ij}), \quad d_{ij} \leq d_{\max}$$

where $\lambda$ is the decay rate and $d_{\max}$ is the hard cutoff beyond which $w_{ij} = 0$.

def compute_decay_weights(
    centroids: np.ndarray,
    tree: cKDTree,
    max_dist: float,
    decay_rate: float
) -> tuple[dict, dict]:
    """
    Build neighbor and weight dicts for libpysal W using exponential decay.

    Parameters
    ----------
    centroids : ndarray, shape (n, 2)
        Centroid coordinates in projected CRS units (meters).
    tree : cKDTree
        Pre-built tree from the same centroids array.
    max_dist : float
        Hard cutoff in the same units as centroids (usually meters).
    decay_rate : float
        Lambda parameter. Recommended starting value: -ln(0.05) / max_dist.

    Returns
    -------
    neighbors : dict  {int: list[int]}
    weights   : dict  {int: list[float]}
    """
    n = len(centroids)
    neighbors: dict[int, list[int]] = {}
    weights: dict[int, list[float]] = {}

    for i in range(n):
        candidate_idx = tree.query_ball_point(centroids[i], r=max_dist)

        nbrs, wts = [], []
        for j in candidate_idx:
            if j == i:
                continue   # diagonal elements must be zero by convention
            d = float(np.linalg.norm(centroids[i] - centroids[j]))
            nbrs.append(j)
            wts.append(float(np.exp(-decay_rate * d)))

        neighbors[i] = nbrs
        weights[i] = wts

    return neighbors, weights

Self-exclusion (if j == i: continue) is mandatory. libpysal does not enforce this automatically; including self-loops inflates diagonal entries and breaks the spatial lag operator $Wy$.

Step 4 — Instantiate libpysal W and apply normalization #

def build_custom_distance_weights(
    gdf: gpd.GeoDataFrame,
    max_dist: float,
    decay_rate: float = 0.5,
    row_standardize: bool = True
) -> W:
    """
    Full pipeline: validate → extract centroids → build tree →
    compute decay weights → instantiate W → normalize.
    """
    gdf = prepare_geodataframe(gdf)
    n_obs = len(gdf)
    centroids = extract_centroids(gdf)
    tree = cKDTree(centroids)

    neighbors, weights = compute_decay_weights(
        centroids, tree, max_dist, decay_rate
    )

    w_obj = W(neighbors, weights, id_order=list(range(n_obs)))

    if row_standardize:
        w_obj.transform = "r"

    return w_obj

Step 5 — Choose the right decay rate #

A practical heuristic: set $\lambda$ so that $w_{ij}$ at max_dist equals your target floor weight $w_{\min}$:

$$\lambda = \frac{-\ln(w_{\min})}{d_{\max}}$$

def recommended_decay_rate(max_dist: float, floor_weight: float = 0.05) -> float:
    """
    Return the lambda that makes exp(-lambda * max_dist) == floor_weight.
    Default floor of 0.05 means edge-of-bandwidth neighbors carry 5% weight.
    """
    return -np.log(floor_weight) / max_dist

# Example: 10 km bandwidth, 5% floor weight
lam = recommended_decay_rate(max_dist=10_000, floor_weight=0.05)
# lam ≈ 0.000300

Step 6 — Validate before use #

def validate_weights(w_obj: W, label: str = "W") -> None:
    """Raise or warn on structural problems that corrupt downstream models."""
    n = w_obj.n
    sparsity = 1.0 - (w_obj.sparse.nnz / (n * n))

    print(f"[{label}] n={n}  islands={len(w_obj.islands)}  "
          f"components={w_obj.n_components}  sparsity={sparsity:.2%}")

    if w_obj.islands:
        raise ValueError(
            f"{len(w_obj.islands)} isolated nodes detected. "
            "Reduce max_dist or apply a k-nearest fallback."
        )
    if w_obj.n_components > 1:
        print(f"  WARNING: {w_obj.n_components} disconnected components. "
              "SAR/SEM models require a fully connected graph.")
    if sparsity < 0.90:
        print(f"  WARNING: matrix density > 10%. "
              "Consider tightening max_dist or switching to KNN.")

validate_weights(w_obj)

Interpreting the Output #

After row-standardization, the spatial lag $Wy$ computes a weighted average of neighbors’ values for each observation. If row sums deviate from 1.0 after .transform = "r", isolated nodes or numerical precision errors are the typical cause — check w_obj.islands first.

Critical Best Practices #

Always project before measuring distance #

cKDTree computes Euclidean distance in the input coordinate space. Degree-based inputs (EPSG:4326) produce distances in degrees — a 0.01-degree difference near the equator is approximately 1.1 km, but at 60° latitude it is approximately 0.55 km. Always call gdf.to_crs("EPSG:XXXXX") before centroid extraction. UTM zones (e.g., EPSG:32618 for UTM 18N) are the most common choice for local analyses; for continental-scale work use an equal-area projection such as EPSG:5070 (Albers Equal Area, contiguous US).

Reset the index before building W #

W keys everything on integers. A DataFrame that retains a non-sequential index from filtering or merging will silently misalign weights with regression variables. gdf.reset_index(drop=True) is a one-liner with no side effects and should be treated as a mandatory precondition.

Use w_obj.sparse for downstream operations #

Never call w_obj.full()[0] on datasets larger than a few thousand observations. The dense $n \times n$ float64 array for $n=10{,}000$ requires 800 MB; for $n=50{,}000$ it requires 20 GB. w_obj.sparse returns a scipy.sparse.csr_matrix that is compatible with spreg, esda, and direct NumPy/SciPy linear algebra without materializing the dense form.

Row-standardize for autoregressive models, skip for flow analysis #

Row-standardization (.transform = "r") normalizes each observation’s outgoing weights to sum to 1.0. This is correct for spatial autoregressive (SAR) and spatial error (SEM) models where the spatial lag must be interpretable as a weighted average. For models that preserve absolute interaction magnitude — migration flow tables, pollutant dispersion with mass conservation — skip row-standardization and consider double-standardization or retain the raw exponential weights.

Implement a k-nearest fallback for sparse datasets #

In datasets with spatially uneven observation density, a fixed max_dist will produce isolated nodes in sparse regions and over-dense neighborhoods in dense regions. A robust production pattern combines both: first attempt the distance-band matrix; if len(w_obj.islands) > 0, fall back to libpysal.weights.KNN.from_dataframe(gdf, k=4) for island observations and patch the two matrices together.

Troubleshooting #

Next Steps #

For a broader understanding of the construction options — contiguity, KNN, kernel, and hybrid topologies — see the parent Spatial Weight Matrices overview. Once your matrix is validated, the most common next task is computing global and local clustering statistics; the how to calculate Moran’s I in PySAL guide shows how to pass a custom W object directly to esda.Moran.

Related #

• Spatial Weight Matrices: Construction & Validation in Python
• How to Calculate Moran’s I in PySAL
• Spatial Autocorrelation Metrics

← Back to Spatial Weight Matrices: Construction & Validation in Python