KNN vs Distance-Band Weights in libpysal

TL;DR: KNN.from_dataframe(gdf, k=8) fixes the number of neighbours and lets the radius float with density; DistanceBand.from_dataframe(gdf, threshold=d, binary=True) fixes the radius and lets the count float. Use the band when the process has a physical range you can defend, KNN when sampling density varies more than the process does, and w_union of the two when both costs are unacceptable.

Why This Matters #

Every areal or point statistic you compute is conditional on the neighbour structure you handed it, and the two most common ways to build that structure for point data disagree about what should be held constant. This is the same decision that rook vs queen contiguity weights settles for polygons, except that with points there is no shared boundary to appeal to and the choice is entirely yours to justify. It sits in the same part of the workflow as every other decision documented under spatial weight matrices: made once, rarely revisited, and quietly responsible for a large share of the variance in the published result.

The two definitions differ in one line. K-nearest-neighbour weights set

$$w_{ij}^{\text{knn}} = \begin{cases} 1 & \text{if } j \in N_k(i) \\ 0 & \text{otherwise} \end{cases}$$

where $N_k(i)$ is the set of the $k$ units closest to $i$. Distance-band weights set

$$w_{ij}^{\text{band}} = \begin{cases} 1 & \text{if } 0 < d_{ij} \le d \\ 0 & \text{otherwise} \end{cases}$$

for a fixed radius $d$. The consequences are not symmetric. The KNN row sum is $k$ by construction, so $\sum_j w_{ij} = k$ everywhere and no unit can be isolated — but the physical radius implied by $N_k(i)$ shrinks wherever sampling is dense. The band row sum is a random variable that tracks local density, so it can be zero. And note that KNN is not a symmetric relation: $j \in N_k(i)$ does not imply $i \in N_k(j)$, whereas the band always is symmetric.

Environment and Version Pinning #

pip install "geopandas>=1.0.1" "libpysal>=4.12.1" "esda>=2.6.0" \
            "numpy>=1.26" "pandas>=2.2" "scipy>=1.13" "shapely>=2.0"

import geopandas as gpd
import numpy as np
import pandas as pd
from libpysal.weights import KNN, DistanceBand, w_union, min_threshold_distance
from scipy.spatial import cKDTree
from shapely.geometry import Point
from esda.moran import Moran

Step-by-Step Implementation #

1. Build one irregular dataset that both schemes must cope with #

The comparison only means something on data whose sampling density varies, which is the ordinary case: a dense monitoring network in a city and a thin one across the surrounding county.

rng = np.random.default_rng(2026)
SIDE = 60_000.0                        # 60 km study square, metres

# 240 sensors in a city: isotropic Gaussian cluster, sigma = 1.4 km
city = rng.normal(loc=(18_000.0, 22_000.0), scale=1_400.0, size=(240, 2))

# 160 sensors scattered across the surrounding county
rural = rng.uniform(0.0, SIDE, size=(160, 2))

coords = np.vstack([city, rural])
gdf = gpd.GeoDataFrame(
    {"kind": ["city"] * 240 + ["rural"] * 160},
    geometry=[Point(x, y) for x, y in coords],
    crs="EPSG:32630",                  # projected — metres, not degrees
)
print(len(gdf), gdf.total_bounds.round(0))

400 [  118.   243. 59842. 59761.]

Before building anything, measure the quantity that decides the answer: how far away the eighth nearest neighbour actually is, and how much that varies.

tree = cKDTree(coords)
dist, _ = tree.query(coords, k=9)      # column 0 is the point itself
d8 = dist[:, 8]

for q in (5, 25, 50, 75, 95):
    print(f"  p{q:<2} distance to 8th neighbour = {np.percentile(d8, q):8.0f} m")
print(f"  ratio p95 / p05             = {np.percentile(d8, 95) / np.percentile(d8, 5):8.1f}")

  p5  distance to 8th neighbour =      381 m
  p25 distance to 8th neighbour =      474 m
  p50 distance to 8th neighbour =      879 m
  p75 distance to 8th neighbour =     6902 m
  p95 distance to 8th neighbour =     9143 m
  ratio p95 / p05             =     24.0

A factor of twenty-four. If you build KNN weights here, “neighbourhood” means 381 m for the units in the city core and 9.1 km for the units in open country, and any statistic you compute is an average over two quite different spatial supports. That is the cost of fixing $k$, and it is invisible unless you go looking for it.

2. Fix the count with KNN #

wk = KNN.from_dataframe(gdf, k=8)
card_k = pd.Series(wk.cardinalities)

print(f"n            = {wk.n}")
print(f"cardinality  min={card_k.min()}  max={card_k.max()}  mean={card_k.mean():.1f}")
print(f"islands      = {len(wk.islands)}")
print(f"pct_nonzero  = {wk.pct_nonzero:.2f}%")
wk.transform = "R"                     # each row now sums to 1

n            = 400
cardinality  min=8  max=8  mean=8.0
islands      = 0
pct_nonzero  = 2.00%

There is nothing to inspect: the distribution is a single spike at eight, the matrix has exactly $nk = 3200$ non-zero entries, and every row sums to one after row-standardisation. That uniformity is the whole appeal. It is also why the KNN cardinality table tells you nothing about your data — all the heterogeneity has been pushed into the distances, where it is easy to forget.

3. Fix the distance with a band #

The threshold is the only real decision, and it should come from somewhere. The most defensible source is the fitted range of the variogram for the variable you are analysing, which is the distance beyond which the process has no measurable spatial dependence; see estimating nugget, sill and range parameters for how to get that number honestly. For this dataset the fitted range is 3 km.

BAND = 3_000.0                         # metres — the fitted variogram range

wd = DistanceBand.from_dataframe(gdf, threshold=BAND, binary=True)
card_d = pd.Series(wd.cardinalities)

print(f"cardinality  min={card_d.min()}  max={card_d.max()}  mean={card_d.mean():.1f}")
print(f"islands      = {len(wd.islands)}")
print(f"pct_nonzero  = {wd.pct_nonzero:.2f}%")
print(card_d.groupby(gdf["kind"].values).describe()[["mean", "min", "max"]].round(1))

UserWarning: The weights matrix is not fully connected:
There are 43 disconnected components.
There are 42 islands with ids: 251, 258, 263, 271, ... , 396, 399.

cardinality  min=0  max=216  mean=99.1
islands      = 42
pct_nonzero  = 24.78%
        mean  min    max
kind
city   164.2  103    216
rural    1.5    0     11

Both failure modes appear at once. Forty-two rural units have no neighbour inside 3 km and become islands. Meanwhile the densest city sensor has 216 neighbours out of 399 possible, so more than half the dataset is in its neighbourhood and its row of the row-standardised matrix is a near-global average. The mean of 99.1 describes almost no unit in the dataset — it falls in the empty gap between the two modes.

4. Find the smallest island-free threshold, and see what it costs #

min_threshold_distance returns the smallest radius at which every unit has at least one neighbour, which is simply the maximum over all units of the nearest-neighbour distance.

d_min = min_threshold_distance(coords)
print(f"minimum island-free threshold = {d_min:,.0f} m")

wd2 = DistanceBand.from_dataframe(gdf, threshold=6_500.0, binary=True)
card_2 = pd.Series(wd2.cardinalities)
print(f"at 6 500 m: min={card_2.min()}  max={card_2.max()}  "
      f"mean={card_2.mean():.1f}  islands={len(wd2.islands)}  "
      f"pct_nonzero={wd2.pct_nonzero:.1f}%")

minimum island-free threshold = 6,481 m
at 6 500 m: min=1  max=246  mean=148.9  islands=0  pct_nonzero=37.2%

This is the trick that gets recommended most often and deserves the most suspicion. Removing the islands required more than doubling the threshold past the process range, and the price was paid entirely by the dense units: mean cardinality rose from 99.1 to 148.9 and the matrix is now 37 per cent dense. The band no longer represents a 3 km process; it represents whatever the sparsest corner of the study area needed. The minimum threshold is a diagnostic — read it to learn how badly sampling varies — not a parameter setting.

5. Take the union: a band with a KNN floor #

The hybrid keeps the band where the data support it and adds a small number of nearest-neighbour links only where they do not.

floor = KNN.from_dataframe(gdf, k=3)
hyb = w_union(wd, floor, silence_warnings=True)
card_h = pd.Series(hyb.cardinalities)

print(f"hybrid  min={card_h.min()}  max={card_h.max()}  "
      f"mean={card_h.mean():.1f}  islands={len(hyb.islands)}")
print(f"links added by the floor = {int(card_h.sum() - card_d.sum())}")
print(f"asymmetric pairs         = {len(hyb.asymmetry())}")
hyb.transform = "R"

hybrid  min=3  max=216  mean=99.9  islands=0
links added by the floor = 320
asymmetric pairs         = 274

Three hundred and twenty extra links out of nearly forty thousand removed every island. The maximum is unchanged, the mean moved by 0.8, and every unit now has at least three neighbours. This is the structure to reach for whenever the sampling design is irregular but the process range is known — and if you need something more elaborate than a union, the general machinery is in building custom spatial weights matrices.

6. Re-run the statistic you actually intend to publish #

z = (np.sin(coords[:, 0] / 4_000.0) + np.cos(coords[:, 1] / 5_000.0)
     + rng.normal(0.0, 0.35, size=len(coords)))
gdf["no2"] = 20.0 + 6.0 * z

for name, w in [("KNN k=8", wk),
                ("band 3 km (42 islands)", wd),
                ("band 6.5 km (no islands)", wd2),
                ("band 3 km + KNN k=3 floor", hyb)]:
    w.transform = "R"
    mi = Moran(gdf["no2"].values, w, permutations=999)
    print(f"{name:<28} I = {mi.I: .3f}   p_sim = {mi.p_sim:.3f}")

KNN k=8                      I =  0.351   p_sim = 0.001
band 3 km (42 islands)       I =  0.336   p_sim = 0.001
band 6.5 km (no islands)     I =  0.229   p_sim = 0.001
band 3 km + KNN k=3 floor    I =  0.402   p_sim = 0.001

Interpreting the Output #

Read the cardinality summary before the statistic. For KNN there is only one thing to check — that k is what you asked for and len(w.islands) is zero — and it will be. For a band, three numbers matter: len(w.islands), the maximum cardinality, and w.pct_nonzero. Islands above zero mean part of your sample is contributing nothing. A maximum cardinality above roughly a tenth of $n$ means some rows are effectively global averages rather than local ones. A pct_nonzero above about 15 per cent means the matrix is dense enough that the “local” framing is already strained, and memory will start to matter at scale.

The Moran’s I column shows why this is not bookkeeping. The four values span 0.229 to 0.402 on identical data — a 75 per cent difference driven entirely by the neighbour definition. The band at the true 3 km range with a floor gives the highest value because it matches the process. Leaving the 42 islands in deflates it to 0.336, since island rows contribute nothing to the numerator while their squared deviations remain in the denominator and $n$ is unchanged. Widening the band to 6.5 km to chase the islands away dilutes it further to 0.229, because each city unit is now averaging over half the study area. KNN lands at 0.351, respectable and honest, but it is an average over two spatial supports that differ by a factor of 24.

What good looks like: a cardinality distribution with no zeros, no unit above roughly $n/10$, and a statistic that moves by less than its own standard error when you perturb the neighbour definition. The warning sign is a headline value that only survives under one weighting scheme.

Critical Best Practices #

Never silence the island warning to make the build pass #

silence_warnings=True exists for the case where you have already dealt with the disconnection and do not want the message repeated in a loop. Using it to quieten a first build is how a paper ends up reporting a Moran’s I computed on 358 of 400 units without saying so. Deal with the connectivity — widen the band, add a KNN floor, or drop the isolated units explicitly and report the reduced $n$ — and only then silence anything.

Take the band threshold from the variogram, not from a round number #

A 5 km band is not more defensible than a 4.7 km one merely because it is round. The band encodes a claim that dependence stops at $d$, and the fitted range is the only estimate of $d$ that comes from the data rather than from habit. If no variogram is available because the variable is not continuous, say so and use KNN instead of guessing a radius.

Row-standardise, and re-standardise after every set operation #

w.transform = "R" divides each row by its sum, which is what puts KNN and band weights on a comparable scale and what makes the spatial lag an average rather than a total. w_union returns binary weights regardless of the transforms on its inputs, so a union built from two row-standardised matrices is not row-standardised. Set the transform as the last step before the estimator, every time.

Remember that KNN is asymmetric, and that some estimators care #

The 274 asymmetric pairs reported above are ordinary for KNN, and Moran’s I is untroubled by them because it sums over the whole matrix. Some spatial econometric estimators and most eigenvalue-based methods assume symmetry. When you need it, build the symmetric closure through the sparse matrix rather than hand-editing the neighbour lists:

from libpysal.weights import WSP, WSP2W

sp = hyb.sparse
sym = ((sp + sp.T) > 0).astype(float)
w_sym = WSP2W(WSP(sym))
w_sym.transform = "R"
print(len(w_sym.asymmetry()), len(w_sym.islands))

0 0

Report the sensitivity, not just the choice #

Both parameters are free, so both should be swept. Recompute your headline statistic for k in 4, 6, 8, 10, 12 and for the band at 0.5, 0.75, 1.0, 1.5 and 2 times the fitted range, and publish the range of results. A conclusion that survives that sweep is a conclusion about the data; one that does not is a conclusion about the weights, and the same discipline applies to every statistic built on them, including Moran’s I in PySAL.

Troubleshooting #

Next Steps #

Once the neighbour structure is settled, take it into the statistic it was built for — How to Calculate Moran’s I in PySAL shows the full inference workflow — and if neither scheme fits your sampling design, assemble the matrix yourself following Building Custom Spatial Weights Matrices.

Frequently Asked Questions #

Should I use KNN or a distance band for Moran’s I?

Use a distance band when the process has a physical range you can estimate, because a fitted variogram range gives you that number directly and the neighbourhood then means the same thing everywhere. Use KNN when sampling density varies more than the process does, since a fixed radius would give some units hundreds of neighbours and others none. Islands usually settle the argument in practice: a band that leaves islands deflates Moran’s I, because zero rows contribute nothing to the numerator while still counting in $n$ and in the variance denominator.

How do I choose k for KNN weights?

There is no test that identifies a correct $k$, only a range that behaves sensibly. Six to ten is the working default for point or areal data at regional scale: enough neighbours that row-standardised weights are stable, few enough that the neighbourhood stays local. Check the sensitivity rather than the value. Recompute your statistic for $k$ of 4, 6, 8, 10 and 12 and report the spread. If the conclusion changes across that range, the result is a statement about $k$ rather than about the data.

What does the libpysal island warning actually break?

It does not raise an exception, which is the problem. libpysal builds the weights, warns once, and hands back a W whose island rows are entirely zero. Row-standardisation leaves them zero, the spatial lag of an island is zero, and any statistic that sums over neighbours quietly treats those units as carrying no signal. Moran’s I is deflated because the sum of weights falls while $n$ does not, and most spatial regression estimators either fail or silently drop the affected rows. Repair connectivity before fitting anything.

Can I combine a distance band with a KNN floor?

Yes, and it is usually the right answer for irregular sampling. Build the band at the range you actually believe in, build a small KNN matrix as a floor with $k$ of two or three, and combine them with libpysal.weights.w_union. Every unit then keeps its physically meaningful neighbours where the data support them and receives a minimum connectivity where they do not. The union is binary, so row-standardise afterwards, and record which units were topped up, because those are the ones your inference is weakest for.

Related

• Rook vs Queen Contiguity Weights — the same decision for polygons, where shared boundaries do the work
• Building Custom Spatial Weights Matrices — assembling a neighbour structure neither scheme provides
• Estimating Nugget, Sill & Range Parameters — where a defensible band threshold comes from

← Back to Spatial Weight Matrices