Row-Standardising and Kernel Weights in libpysal

TL;DR: w.transform rewrites a libpysal weights object in place. Use w.transform = "R" for Moran’s I and spatial lag models, so that w.sparse @ y is a neighbourhood mean on the same scale as y; use "B" for Getis-Ord Gi*, which compares neighbourhood sums. For distance decay inside a band, build Kernel.from_dataframe(pts, bandwidth=1800, function="triangular").

Why This Matters #

Choosing neighbours is only half of building a weights matrix. Once you have decided who counts as a neighbour — the subject of KNN vs distance-band weights in libpysal — you still have to decide how much each of them counts. That second decision lives in a single attribute, w.transform, and it is the one most likely to be left at whatever the last function call set it to. Two analysts with identical neighbour lists can report Moran’s I of 0.59 and 0.67 on the same data, and both are computing the statistic correctly; they have simply standardised differently.

The transform is not a cosmetic choice. It changes what the spatial lag is, which changes the scale the statistic sits on, which changes whether a published threshold means anything. This page walks the four transformations libpysal offers, shows the same weights object under each, and then moves to kernel weights, where the decay function replaces the transform as the thing that decides how much a distant neighbour is worth. Everything here sits inside Spatial Weight Matrices, which is in turn part of Core Concepts of Spatial Statistics & Geostatistics.

Environment and Version Pinning #

pip install "libpysal==4.15.0" "esda==2.10.0" "geopandas==1.1.4" \
            "numpy==2.5.1" "scipy==1.18.0" "shapely==2.1.2" "mgwr==2.2.1"

import warnings

import geopandas as gpd
import numpy as np
from shapely.geometry import box

from libpysal.weights import Queen, Kernel, fill_diagonal
from esda.moran import Moran
from esda.getisord import G_Local

Step-by-Step Implementation #

1. Build a lattice and its Queen weights #

A regular lattice is the right test bed because the cardinalities are known in advance: corner cells have three Queen neighbours, edge cells five, interior cells eight. Any change in a statistic is then attributable to the transform rather than to irregular geometry.

CELL = 1000.0                                   # 1 km cells, metric CRS
cells = [box(x * CELL, y * CELL, (x + 1) * CELL, (y + 1) * CELL)
         for y in range(8) for x in range(8)]
gdf = gpd.GeoDataFrame(geometry=cells, crs="EPSG:32630")

rng = np.random.default_rng(11)
cx = gdf.geometry.centroid.x.values / 1000.0
cy = gdf.geometry.centroid.y.values / 1000.0
gdf["no2"] = 22.0 + 1.8 * (cx + 0.5 * cy) + rng.normal(0, 3.0, len(gdf))
y = gdf["no2"].values

w = Queen.from_dataframe(gdf, use_index=True)
print(f"n = {w.n}, links = {int(w.s0)}, cardinality 0 = {w.cardinalities[0]}, "
      f"cardinality 27 = {w.cardinalities[27]}")
print(f"transform as built: {w.transform!r}")

n = 64, links = 420, cardinality 0 = 3, cardinality 27 = 8
transform as built: 'O'

Note the transform at construction: 'O', for original, not 'B'. For contiguity weights the original values happen to be ones, so 'O' and 'B' coincide numerically — but for k-nearest or kernel weights they do not, and code that assumes a freshly built object is binary will be wrong.

2. Look at one row under all four transforms #

w.transform is a property with a setter. Assigning to it rebuilds every weight from the stored original values and resets the cached scalars such as w.s0, the sum of all weights.

for t in ("B", "R", "D", "V"):
    ww = Queen.from_dataframe(gdf, use_index=True)   # a fresh object each time
    ww.transform = t
    r0 = np.round(ww.weights[0], 6)                  # corner cell, 3 neighbours
    r27 = np.round(ww.weights[27], 6)                # interior cell, 8 neighbours
    print(f"{t}  row 0  {r0}  sum={r0.sum():.4f}")
    print(f"   row 27 {r27[:4]} ... sum={r27.sum():.4f}   s0={ww.s0:.4f}")

B  row 0  [1. 1. 1.]  sum=3.0000
   row 27 [1. 1. 1. 1.] ... sum=8.0000   s0=420.0000
R  row 0  [0.333333 0.333333 0.333333]  sum=1.0000
   row 27 [0.125 0.125 0.125 0.125] ... sum=1.0000   s0=64.0000
D  row 0  [0.002381 0.002381 0.002381]  sum=0.0071
   row 27 [0.002381 0.002381 0.002381 0.002381] ... sum=0.0190   s0=1.0000
V  row 0  [0.227503 0.227503 0.227503]  sum=0.6825
   row 27 [0.139317 0.139317 0.139317 0.139317] ... sum=1.1145   s0=64.0000

Written out, B sets every neighbour weight to 1, so $S_0$ is the number of directed links. R divides each row by its own row sum, so every row sums to 1 and $S_0 = n$. D divides every weight by the global sum, so the whole matrix sums to 1 while the relative weights are untouched. V divides each row by its Euclidean norm $q_i = \sqrt{\sum_j w_{ij}^2}$ and then rescales globally so that $S_0 = n$; the consequence is that $\sum_j w_{ij}^2$ is identical for every unit, which is what stabilises the variance of the lag. Row sums under V grow as $\sqrt{k_i}$: $0.6825 = 0.394\sqrt{3}$ and $1.1145 = 0.394\sqrt{8}$.

3. Watch the spatial lag and Moran’s I move #

The spatial lag is w.sparse @ y. Print it alongside Moran and the difference between the transforms stops being abstract.

for t in ("B", "R", "D", "V"):
    ww = Queen.from_dataframe(gdf, use_index=True)
    ww.transform = t
    lag = ww.sparse @ y
    mi = Moran(y, ww, transformation=t, permutations=999)
    print(f"{t}  s0={ww.s0:9.4f}  I={mi.I:.4f}  E[I]={mi.EI:.4f}  "
          f"z={mi.z_norm:6.3f}  lag[0]={lag[0]:8.3f}  mean(lag)={lag.mean():7.3f}")
print(f"mean(y) = {y.mean():.3f}")

B  s0= 420.0000  I=0.5910  E[I]=-0.0159  z= 9.362  lag[0]=  76.228  mean(lag)=214.044
R  s0=  64.0000  I=0.6672  E[I]=-0.0159  z=10.162  lag[0]=  25.409  mean(lag)= 32.663
D  s0=   1.0000  I=0.5910  E[I]=-0.0159  z= 9.362  lag[0]=   0.181  mean(lag)=  0.510
V  s0=  64.0000  I=0.6262  E[I]=-0.0159  z= 9.803  lag[0]=  17.342  mean(lag)= 32.639

Moran’s I is

$$I = \frac{n}{S_0} \cdot \frac{\mathbf{z}^{\top}\mathbf{W}\mathbf{z}}{\mathbf{z}^{\top}\mathbf{z}}, \qquad \mathbf{z} = \mathbf{y} - \bar{y}.$$

Because $S_0$ appears in the denominator, multiplying every weight by a constant $c$ multiplies both $S_0$ and the numerator by $c$ and leaves $I$ untouched — which is exactly why B and D give the identical 0.5910. D is nothing but B divided by 420. R and V are row-wise rescalings, not global ones, so they genuinely change the answer. Under R, $S_0 = n$ and the leading factor disappears entirely, leaving $I$ as the slope of the ordinary least-squares regression of the lag on $\mathbf{z}$ — the line drawn on every Moran scatterplot. That equality is the practical reason R is the convention for calculating Moran’s I in PySAL, and the same argument applies to the spatial lag term $\rho \mathbf{W}\mathbf{y}$ in a lag model, where $\rho$ is interpretable only when $\mathbf{W}\mathbf{y}$ is a neighbourhood average.

Look at mean(lag) in each row. Under R it is 32.663 against a data mean of 32.638: the lag lives on the same scale as the variable. Under B it is 214.044, six and a half times larger, because it is a sum over a variable number of neighbours. Under D it is 0.510, three orders of magnitude too small.

4. Why Getis-Ord Gi* wants binary weights instead #

Gi* is a different animal. It compares the total of a neighbourhood against the total of the whole map:

$$G_i^{*} = \frac{\sum_j w_{ij} x_j}{\sum_j x_j}, \qquad \mathbb{E}[G_i^{*}] = \frac{W_i}{n}, \quad W_i = \sum_j w_{ij}.$$

esda standardises this with the classic Getis-Ord variance, which substitutes $\sum_j w_{ij}^2 = W_i$ — an identity that holds only when the weights are zero or one. Row-standardise and $W_i = 1$ for every unit, so the expected value and the variance become the same for all $n$ units and the statistic loses every trace of how many neighbours a unit has.

for t in ("B", "R"):
    wg = Queen.from_dataframe(gdf, use_index=True)
    wg.transform = t
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")        # esda notes how it sets the self-weight
        g = G_Local(y, wg, transform=t, star=True, permutations=0)
    print(f"{t}  G*[0]={g.Gs[0]:.5f}  E[G*][0]={g.EGs[0]:.5f}  z[0]={g.Zs[0]:7.4f}  "
          f"z range=[{g.Zs.min():.3f}, {g.Zs.max():.3f}]  "
          f"|z|>1.96: {int((abs(g.Zs) > 1.96).sum())}")

B  G*[0]=0.04772  E[G*][0]=0.06250  z[0]=-3.1392  z range=[-4.152, 3.959]  |z|>1.96: 31
R  G*[0]=0.01193  E[G*][0]=0.01562  z[0]=-1.5318  z range=[-1.626, 1.516]  |z|>1.96: 0

Thirty-one significant units become zero. Nothing about the data or the neighbour lists changed. This is the single most consequential transform mistake in practice, and it is silent: no exception, no warning about significance, just a map with no hot spots on it. Always set transform="B" when running Getis-Ord Gi* hot spot analysis in Python.

5. Kernel weights: decay inside the band #

A distance band gives every neighbour inside the radius the same weight and everyone outside nothing, which is a strange model of a smooth process. A kernel keeps the band but grades the weight by $z_{ij} = d_{ij} / b_i$, the distance as a fraction of the bandwidth.

pts = gdf.copy()
pts["geometry"] = gdf.geometry.centroid

for fn in ("triangular", "uniform", "quadratic", "gaussian"):
    kw = Kernel.from_dataframe(pts, bandwidth=1800.0, function=fn,
                               silence_warnings=True)
    print(f"{fn:10s} k={kw.cardinalities[0]}  w[0]={np.round(kw.weights[0], 4)}"
          f"  s0={kw.s0:.3f}")

triangular k=4  w[0]=[1.     0.4444 0.2143 0.4444]  s0=205.563
uniform    k=4  w[0]=[0.5 0.5 0.5 0.5]  s0=242.000
quadratic  k=4  w[0]=[0.75   0.5185 0.287  0.5185]  s0=220.407
gaussian   k=4  w[0]=[0.3989 0.3419 0.293  0.3419]  s0=159.544

The corner cell has four entries: itself at distance 0, two rook neighbours at 1000 m ($z = 0.556$) and one diagonal neighbour at 1414 m ($z = 0.786$). libpysal implements $K(z) = 1 - z$ for triangular, a flat $1/2$ for uniform, $\tfrac{3}{4}(1 - z^2)$ for quadratic, $\tfrac{15}{16}(1 - z^2)^2$ for quartic and $(2\pi)^{-1/2}e^{-z^2/2}$ for gaussian, all evaluated only for $|z| \le 1$. The neighbour set here is identical to Queen contiguity; what differs is that the diagonal neighbour is worth 0.214 against 0.444 for the rook neighbours under a triangular kernel, but 0.293 against 0.342 under a gaussian one. Uniform is the degenerate case: it is a binary distance band wearing a kernel’s clothes.

6. Fixed against adaptive bandwidth, and the diagonal #

fixed=True gives every unit the same bandwidth in map units; fixed=False gives every unit the distance to its k-th nearest neighbour. When bandwidth is omitted, libpysal derives one from k — the maximum k-th-neighbour distance over all units for a fixed kernel, which is a worst-case radius.

ka = Kernel.from_dataframe(pts, fixed=False, k=8, function="triangular",
                           silence_warnings=True)
kf = Kernel.from_dataframe(pts, fixed=True, k=8, function="triangular",
                           silence_warnings=True)
print(f"adaptive: bw[0]={ka.bandwidth[0][0]:.1f}  bw[27]={ka.bandwidth[27][0]:.1f}  "
      f"cards {min(ka.cardinalities.values())}-{max(ka.cardinalities.values())}")
print(f"fixed   : bw[0]={kf.bandwidth[0][0]:.1f}  bw[27]={kf.bandwidth[27][0]:.1f}  "
      f"cards {min(kf.cardinalities.values())}-{max(kf.cardinalities.values())}")
print("adaptive w[27] =", np.round(ka.weights[27], 4))

adaptive: bw[0]=2828.4  bw[27]=1414.2  cards 9-9
fixed   : bw[0]=2828.4  bw[27]=2828.4  cards 9-25
adaptive w[27] = [1.     0.2929 0.2929 0.2929 0.2929 0.     0.     0.    0.    ]

The adaptive kernel gives every unit exactly nine entries; the fixed one gives the corner cell nine and an interior cell twenty-five. Note also that k=8 counts eight neighbours besides the unit itself, and that the k-th neighbour sits exactly at the bandwidth, so a triangular kernel gives it a weight of zero. Four of cell 27’s nine entries are therefore worth nothing.

The diagonal is the other thing to get right. Every kernel includes each unit as its own neighbour at distance 0, so the self-weight is $K(0)$: 1 for triangular, 0.5 for uniform, 0.399 for gaussian. diagonal=True overrides that with 1.0.

kg = Kernel.from_dataframe(pts, bandwidth=1800.0, function="gaussian",
                           silence_warnings=True)
kd = Kernel.from_dataframe(pts, bandwidth=1800.0, function="gaussian",
                           diagonal=True, silence_warnings=True)
print("gaussian    w[0] =", np.round(kg.weights[0], 4))
print("diagonal    w[0] =", np.round(kd.weights[0], 4))
kd.transform = "R"
print("diagonal, R w[0] =", np.round(kd.weights[0], 4),
      " sum =", round(sum(kd.weights[0]), 4))

gaussian    w[0] = [0.3989 0.3419 0.293  0.3419]
diagonal    w[0] = [1.     0.3419 0.293  0.3419]
diagonal, R w[0] = [0.5059 0.173  0.1482 0.173 ]  sum = 1.0

Interpreting the Output #

Read w.s0 first: it tells you which transform is actually in force without printing a single weight. s0 == n means R or V, s0 == 1 means D, and s0 equal to the link count means B or an untouched binary object. Then check mean(w.sparse @ y) against y.mean(). If they agree to a couple of decimal places, the lag is a neighbourhood average and any statistic built on it will be on its conventional scale. If they differ by a factor of five or a factor of a thousand, you are about to interpret a sum or a fraction as if it were a mean.

Good looks like this: w.transform == "R" for anything with a lag in it, "B" for Getis-Ord, Moran’s I between $-1$ and $1$ with $E[I] = -1/(n-1)$, and a kernel object whose row weights fall monotonically with distance. The warning signs are Moran’s I above 1 (almost always a non-row-standardised matrix), a Gi* map with no significant units where an obvious gradient exists (row-standardised weights), and a kernel row whose last few weights are exactly zero (an adaptive bandwidth landing on the k-th neighbour).

Critical Best Practices #

w.transform mutates, and esda re-transforms behind you #

Assigning to w.transform rewrites the object you hold, not a copy. Worse, esda.moran.Moran takes transformation="r" by default and applies it to the object you passed in, so a weights matrix you carefully set to binary comes back row-standardised.

wm = Queen.from_dataframe(gdf, use_index=True)
wm.transform = "B"
print("before Moran:", wm.transform, round(wm.s0, 3))
_ = Moran(y, wm)                       # default transformation="r"
print("after  Moran:", wm.transform, round(wm.s0, 3))

before Moran: B 420.0
after  Moran: R 64.0

G_Local does the same with its own transform argument. In a pipeline that computes several statistics, build a fresh weights object per statistic or reset the transform immediately before each call — never assume the object survived the last one.

Transformations chain from the original, except that D does not #

libpysal derives each transform from the stored 'O' weights, so w.transform = "R" followed by w.transform = "B" gives clean binary weights. The 'D' branch is the exception: it reads w.s0 before it restores the originals, so the constant it divides by is the current transform’s total rather than the original one.

wd = Queen.from_dataframe(gdf, use_index=True)
wd.transform = "R"
wd.transform = "D"
print("R then D  s0 =", round(wd.s0, 4))    # should be 1.0
wd2 = Queen.from_dataframe(gdf, use_index=True)
wd2.transform = "D"
print("fresh D   s0 =", round(wd2.s0, 4))

R then D  s0 = 6.5625
fresh D   s0 = 1.0

If you need D, apply it to a freshly constructed object.

B and D are the same test, so do not report them as a robustness check #

Because Moran’s I divides by $S_0$, any uniform rescaling of the matrix cancels. Running the statistic under B and again under D and finding 0.5910 both times demonstrates nothing except that the arithmetic is right. A genuine sensitivity check varies the neighbour definition — contiguity against k-nearest against a distance band — or contrasts a row-wise transform (R, V) with a global one.

Zero the kernel diagonal before computing Moran’s I #

A kernel weights object always includes each unit as its own neighbour, so w.sparse @ y contains $y_i$ itself and the statistic is partly measuring a variable against itself.

kt = Kernel.from_dataframe(pts, bandwidth=1800.0, function="triangular",
                           silence_warnings=True)
print("self-weight kept :", round(Moran(y, kt).I, 4))
kz = fill_diagonal(
    Kernel.from_dataframe(pts, bandwidth=1800.0, function="triangular",
                          silence_warnings=True), 0.0)
print("zero diagonal    :", round(Moran(y, kz).I, 4))

self-weight kept : 0.7743
zero diagonal    : 0.6686

That is sixteen points of inflation from a term that should not be there. The zero-diagonal figure of 0.6686 is within a whisker of the row-standardised Queen result of 0.6672, which is the honest comparison.

Do not hand pre-transformed kernels to GWR #

Geographically weighted regression uses the same distance-decay machinery, but mgwr builds and row-standardises its own kernel at every regression point and selects the bandwidth by cross-validation or AICc. Pass coordinates, not a W:

from mgwr.sel_bw import Sel_BW
from mgwr.gwr import GWR

coords = np.column_stack([gdf.geometry.centroid.x, gdf.geometry.centroid.y])
Y = gdf["no2"].values.reshape(-1, 1)
X = traffic.reshape(-1, 1)

bw = Sel_BW(coords, Y, X, fixed=False, kernel="bisquare").search(criterion="AICc")
res = GWR(coords, Y, X, bw, fixed=False, kernel="bisquare").fit()

With fixed=False, bw is a neighbour count, not a distance — the same adaptive idea as Kernel(fixed=False, k=...), and it is worth reporting alongside the coefficients.

Troubleshooting #

Next Steps #

Transformation only matters once the neighbour set is settled, so if you are still choosing between a fixed radius and a fixed count, work through KNN vs distance-band weights in libpysal first, then return here to decide the weighting. With both settled, take the row-standardised matrix into How to Calculate Moran’s I in PySAL and the binary one into Getis-Ord Gi* Hot Spot Analysis in Python.

Frequently Asked Questions #

Why does Moran’s I change under row standardisation but not under double standardisation?

Moran’s I carries the factor $n$ divided by $S_0$, the sum of all weights. Multiplying every weight by the same constant multiplies both the numerator and $S_0$ by that constant, so the ratio is unchanged. Double standardisation divides every weight by $S_0$, a single global constant, so it is exactly that kind of uniform rescaling and returns the binary answer. Row standardisation divides each row by its own row sum, a different constant per unit, which genuinely reweights the evidence and moves the statistic.

Should I ever use the V variance-stabilising transform?

Rarely, and only deliberately. V divides each row by its Euclidean norm and then rescales globally so that $S_0$ equals $n$. The effect is that the sum of squared weights is identical for every unit, which equalises the sampling variance of the spatial lag when cardinalities differ sharply. That is useful when a few units have many more neighbours than the rest. The cost is that row sums are no longer 1, so the lag is not a neighbourhood mean and most published results are not comparable.

Does a kernel weights object need row-standardising?

For Moran’s I and spatial lag models, yes, and you should zero the diagonal first. A libpysal Kernel object includes each unit as its own neighbour, so the raw lag contains the unit’s own value and inflates the statistic. On the lattice used here, a triangular kernel gives Moran’s I of 0.7743 with the self-weight kept and 0.6686 once the diagonal is zeroed with libpysal.weights.fill_diagonal. For Geographically Weighted Regression the kernel is handled internally and you should not pre-transform anything.

Should I pick a fixed or an adaptive bandwidth?

Use a fixed bandwidth when the process has a physical range in metres and the observations are roughly evenly spaced, because every unit then sees the same geographic footprint. Use an adaptive bandwidth when density varies, because a fixed radius gives dense areas many neighbours and sparse areas almost none. On the example lattice, a fixed bandwidth chosen from k equal to 8 produced cardinalities from 9 to 25, while the adaptive version gave exactly 9 everywhere.

Related

• KNN vs Distance-Band Weights in libpysal — choosing the neighbour set before you choose the weighting
• How to Calculate Moran’s I in PySAL — the statistic that assumes row-standardised weights
• Getis-Ord Gi* Hot Spot Analysis in Python — the statistic that assumes binary weights

← Back to Spatial Weight Matrices