Spatial Autocorrelation & Moran's I in Python

msg = "ready to learn spatial statistics"
print("You said:", msg)

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

import json
STATES = "https://cdn.jsdelivr.net/gh/PublicaMundi/MappingAPI@master/data/geojson/us-states.json"
gj = json.load(open_url(STATES))
print("AREAS  states:", len(gj["features"]), "polygons")
print("first feature props:", gj["features"][0]["properties"])

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(4.8, 3.8))
sc = ax.scatter(topo.x, topo.y, c=topo.z, cmap="viridis", s=80, edgecolor="k", linewidth=.4)
ax.set(title="Elevation soundings", xlabel="x", ylabel="y"); ax.set_aspect("equal")
plt.colorbar(sc, ax=ax, label="elevation")

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

y    = np.log(dens)
cent = np.array([state_centroid(f) for f in feats])
print("y ranges", round(float(y.min()),2), "to", round(float(y.max()),2))
print("centroids:", cent.shape)

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

from scipy.spatial import cKDTree
import numpy as np

def knn_W(coords, k=6, row_standardize=True):
    n = len(coords); tree = cKDTree(coords)
    _, idx = tree.query(coords, k=k+1)        # k+1: the nearest point is itself
    W = np.zeros((n, n))
    for i in range(n):
        for j in idx[i, 1:]:                  # skip self (column 0)
            W[i, j] = 1.0
    if row_standardize:
        W = W / W.sum(1, keepdims=True)
    return W

W = knn_W(cent, k=6)
print("W is", W.shape, "| every row sums to 1:", np.allclose(W.sum(1), 1))

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

ylag = W @ y
print("corr(value, neighbour average):", round(float(np.corrcoef(y, ylag)[0, 1]), 3))

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

import numpy as np
def morans_I(y, W):
    z  = y - y.mean()
    S0 = W.sum()
    return (len(y)/S0) * np.sum(W * np.outer(z, z)) / np.sum(z**2)

I = morans_I(y, W)
print("Moran's I (log density) =", round(float(I), 3))

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

rng = np.random.default_rng(42)
ref = np.array([morans_I(rng.permutation(y), W) for _ in range(999)])
p   = (np.sum(ref >= I) + 1) / (999 + 1)
print(f"pseudo p-value = {p:.3f}")
print("expected if random =", round(-1/(len(y)-1), 3))

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.

import matplotlib.pyplot as plt
zy   = (y - y.mean())/y.std()
zlag = W @ zy
fig, ax = plt.subplots(figsize=(4.8, 4.8))
ax.axhline(0, color="#aaa", lw=.8); ax.axvline(0, color="#aaa", lw=.8)
ax.scatter(zy, zlag, s=32, color="#14233D")
slope = np.polyfit(zy, zlag, 1)[0]
xs = np.linspace(zy.min(), zy.max(), 2)
ax.plot(xs, slope*xs, color="#D6552F", lw=2, label=f"slope = Moran's I = {slope:.3f}")
ax.set(xlabel="value (standardised)", ylabel="spatial lag", title="Moran scatterplot"); ax.legend()

In the interactive version this line is blanked out and you write it yourself, with a hint, a solution, and an automatic check.