critwire · Learn · Learn's data

Global Moran's I s4

Prerequisites: Spatial weights (s3)

One number for "do similar values sit together?", and how to test it: global Moran's I, its expectation, variance and z-score in plain SQL, the same statistic in PySAL, a pair-by-pair check that both use the same neighbours, and a permutation test. Everything runs in your own PostgreSQL with PostGIS and the Python environment, on the lab tables (the setup page), with the Tasmanian devil's records.

Concept

One number for "do similar values sit together?", and how to test it:

Why the band is twice the edge: touching hexagons' centroids are 1.73 edges apart and the next ring's 3, so a 2-edge band holds exactly the six touching cells

Worked example

The cells, and a check of their geometry.

Points can't be the units for Moran's I: a record has no value to compare with its neighbours', and many records share one spot. So the records are counted into equal-area hexagons over Tasmania's main island (EPSG:3577, Australian Albers), at two cell sizes, and each cell's value is y = ln(1 + n). Its centroid is kept as two plain numbers, cx and cy, so SQL and PySAL measure distances from the same coordinates. The neighbour rule is a view: two cells of one grid are neighbours when their centroids are within twice the edge length. Unlike most blocks, this one keeps its work: run it without a ROLLBACK, and its tables stay until the lab data is loaded again. It drops them first, so running it again rebuilds them, and every download that reads them starts with it.

Checked 2026-11-02 on PostgreSQL 16.15, PostGIS 3.4.2, GEOS 3.12.1, GDAL 3.8.4; Python 3.12.3 with requirements.txt sha256 d5bbe28024cc.

s4-cells · PostgreSQL with PostGIS 

DROP VIEW IF EXISTS s4_neighbours;
DROP TABLE IF EXISTS s4_cells;
CREATE TABLE s4_cells AS
WITH island AS (
  SELECT (part).geom
  FROM (SELECT ST_Dump(ST_Transform(geom, 3577)) AS part
        FROM admin1_shapes WHERE country_code = 'AU' AND name = 'Tasmania') AS parts
  ORDER BY ST_Area((part).geom) DESC
  LIMIT 1
), devils AS (
  SELECT ST_Transform(o.geom, 3577) AS geom
  FROM occurrences AS o JOIN species AS s ON s.id = o.species_id
  WHERE s.scientific_name = 'Sarcophilus harrisii'
), grid AS (
  SELECT e.edge_m, h.i, h.j, h.geom
  FROM island, (VALUES (5000), (10000)) AS e (edge_m), ST_HexagonGrid(e.edge_m, island.geom) AS h
  WHERE ST_Intersects(h.geom, island.geom)
)
SELECT row_number() OVER (ORDER BY g.edge_m, g.i, g.j)::integer AS id,
       g.edge_m, g.geom,
       ST_X(ST_Centroid(g.geom)) AS cx,
       ST_Y(ST_Centroid(g.geom)) AS cy,
       count(d.geom)::integer AS n,
       ln(1 + count(d.geom)::float8) AS y
FROM grid AS g LEFT JOIN devils AS d ON ST_Contains(g.geom, d.geom)
GROUP BY g.edge_m, g.i, g.j, g.geom;
ALTER TABLE s4_cells ADD PRIMARY KEY (id);
CREATE VIEW s4_neighbours AS
SELECT a.edge_m, a.id AS i, b.id AS j
FROM s4_cells AS a
JOIN s4_cells AS b
  ON b.edge_m = a.edge_m AND b.id <> a.id
 AND ST_DWithin(ST_Point(a.cx, a.cy, 3577), ST_Point(b.cx, b.cy, 3577), 2 * a.edge_m);
SELECT c.edge_m, count(*) AS cells, sum(c.n) AS records_in_cells,
       (SELECT count(*) FROM occurrences AS o JOIN species AS s ON s.id = o.species_id
        WHERE s.scientific_name = 'Sarcophilus harrisii') - sum(c.n) AS records_in_no_cell,
       max(c.n) AS max_n
FROM s4_cells AS c
GROUP BY c.edge_m
ORDER BY c.edge_m;
-- result: s4-cells · 2 rows · sha256 29efe20842eb
edge_m | cells | records_in_cells | records_in_no_cell | max_n
5000 | 1094 | 1972 | 28 | 256
10000 | 296 | 1990 | 10 | 475

In s4.sql, with the blocks it runs after.

The island gives 1,094 cells of 5 km and 296 of 10 km. 1,972 of the 2,000 records fall in a 5 km cell and 1,990 in a 10 km one; the fullest cell holds 256 records at 5 km and 475 at 10 km. ST_HexagonGrid returns the hexagons covering the island's bounding box, with their column and row (i, j), so WHERE ST_Intersects keeps the ones touching land (Grids, g6). ST_Contains puts a record in the one cell that contains it.

Before building weights from a distance band, look at the distances. In a hexagon grid with edge e, centroids are √3·e apart (the touching cells), then 3·e (the next ring). A threshold of 2·e sits between the two, with no pair near it, so every cell gets exactly its touching cells and a rounding error can't move a pair across. The block also counts duplicate centroids and cells with no neighbour (an "island"), which would break the statistic.

s4-geometry · PostgreSQL with PostGIS 

WITH pairs AS (
  SELECT a.edge_m, ST_Distance(ST_Point(a.cx, a.cy, 3577), ST_Point(b.cx, b.cy, 3577)) AS d_m
  FROM s4_cells AS a
  JOIN s4_cells AS b
    ON b.edge_m = a.edge_m AND b.id > a.id
   AND ST_DWithin(ST_Point(a.cx, a.cy, 3577), ST_Point(b.cx, b.cy, 3577), 4 * a.edge_m)
), classes AS (
  SELECT edge_m, round(d_m::numeric, 2) AS d_m,
         dense_rank() OVER (PARTITION BY edge_m ORDER BY round(d_m::numeric, 2)) AS class
  FROM pairs
)
SELECT c.edge_m,
       2 * c.edge_m AS threshold_m,
       (SELECT min(d_m) FROM classes WHERE edge_m = c.edge_m AND class = 1) AS class_1_m,
       (SELECT min(d_m) FROM classes WHERE edge_m = c.edge_m AND class = 2) AS class_2_m,
       count(*) - count(DISTINCT (c.cx, c.cy)) AS duplicate_centroids,
       count(*) FILTER (WHERE NOT EXISTS (SELECT 1 FROM s4_neighbours AS nb WHERE nb.i = c.id)) AS islands,
       (SELECT count(*) FROM pairs AS p WHERE p.edge_m = c.edge_m AND abs(p.d_m - 2 * p.edge_m) <= 1) AS pairs_near_threshold
FROM s4_cells AS c
GROUP BY c.edge_m
ORDER BY c.edge_m;
-- result: s4-geometry · 2 rows · sha256 099b07470829
edge_m | threshold_m | class_1_m | class_2_m | duplicate_centroids | islands | pairs_near_threshold
5000 | 10000 | 8660.25 | 15000.00 | 0 | 0 | 0
10000 | 20000 | 17320.51 | 30000.00 | 0 | 0 | 0

In s4.sql, with the blocks it runs after.

The grid passes at both sizes: the nearest centroids are 8,660.25 m and 17,320.51 m apart, the next class 15,000 and 30,000 m, with no duplicate centroid, no island and no pair within 1 m of the threshold. dense_rank() numbers the distinct distances in order, so class 1 is the smallest and class 2 the next; rounding to 0.01 m first keeps 8,660.254 and 8,660.2540001 in one class.

Lab

  1. Global Moran's I in plain SQL, with row-standardized weights (w_ij = 1/k_i, where k_i is cell i's number of neighbours) and its expected value, variance and z-score under the normality assumption. These are the formulas PySAL's esda.Moran uses (EI, VI_norm, z_norm):

    s4-moran-sql · PostgreSQL with PostGIS 

    WITH z AS (
  SELECT id, edge_m, y - avg(y) OVER (PARTITION BY edge_m) AS z FROM s4_cells
), w AS (
  SELECT edge_m, i, j, 1 / count(*) OVER (PARTITION BY i)::float8 AS w FROM s4_neighbours
), pairs AS (
  SELECT w.edge_m, w.i, w.j, w.w AS w_ij, t.w AS w_ji
  FROM w JOIN w AS t ON t.i = w.j AND t.j = w.i
), sums AS (
  SELECT p.edge_m,
         sum(p.w_ij * zi.z * zj.z) AS wzz,
         sum(p.w_ij) AS s0,
         sum((p.w_ij + p.w_ji) ^ 2) / 2 AS s1
  FROM pairs AS p JOIN z AS zi ON zi.id = p.i JOIN z AS zj ON zj.id = p.j
  GROUP BY p.edge_m
), s2 AS (
  SELECT edge_m, sum(ws ^ 2) AS s2
  FROM (SELECT edge_m, i, sum(w_ij) + sum(w_ji) AS ws FROM pairs GROUP BY edge_m, i) AS cell
  GROUP BY edge_m
), cells AS (
  SELECT edge_m, count(*)::float8 AS n, sum(z * z) AS zz FROM z GROUP BY edge_m
), moran AS (
  SELECT c.edge_m, c.n,
         c.n / s.s0 * s.wzz / c.zz AS i,
         -1 / (c.n - 1) AS ei,
         (c.n ^ 2 * s.s1 - c.n * s2.s2 + 3 * s.s0 ^ 2) / ((c.n - 1) * (c.n + 1) * s.s0 ^ 2) - (1 / (c.n - 1)) ^ 2 AS var_norm
  FROM cells AS c JOIN sums AS s USING (edge_m) JOIN s2 USING (edge_m)
)
SELECT edge_m, n::integer AS n, i, ei, var_norm, (i - ei) / sqrt(var_norm) AS z
FROM moran
ORDER BY edge_m;
    -- result: s4-moran-sql · 2 rows · sha256 7cd3a5f317d7
edge_m | n | i | ei | var_norm | z
5000 | 1094 | 0.27329612511604184 | -0.0009149130832570906 | 0.000326399057870225 | 15.177857980449323
10000 | 296 | 0.2035877773186644 | -0.003389830508474576 | 0.0012597916525627741 | 5.831415676740585

    In s4.sql, with the blocks it runs after.

    Read it one CTE at a time. z takes each grid's mean from its values. w gives every pair its row-standardized weight: count(*) OVER (PARTITION BY i) is k_i. pairs puts w_ij and w_ji on one row, since they differ when two neighbours have different numbers of neighbours. sums makes the numerator and S0, and S1 (half the sum of (w_ij + w_ji)²); s2 makes S2 (the sum over cells of each cell's row total plus column total, squared). moran combines them into I, E[I] and the variance under normality. At 5 km, I is 0.2733, far above E[I] (−0.0009): z is 15.18. At 10 km, I is 0.2036 and z 5.83.

  2. The same statistic in PySAL, from the same cells. DistanceBand builds the binary weights; transformation="r" row-standardizes them. z_rand is the z-score under the randomization assumption, which ArcGIS reports:

    s4-moran-pysal · Python

    Runs after: py-connect

    import esda
import numpy as np
from libpysal.weights import DistanceBand

results = {}
with connect() as conn:
    edges = [e for (e,) in conn.execute("SELECT DISTINCT edge_m FROM s4_cells ORDER BY edge_m")]
    for edge_m in edges:
        rows = conn.execute(
            "SELECT id, y, cx, cy FROM s4_cells WHERE edge_m = %s ORDER BY id", (edge_m,)
        ).fetchall()
        ids = [r[0] for r in rows]
        y = np.array([r[1] for r in rows])
        xy = np.array([(r[2], r[3]) for r in rows])
        w = DistanceBand(xy, threshold=2 * edge_m, binary=True, ids=ids)
        assert w.islands == [], f"cells with no neighbour: {w.islands}"
        results[edge_m] = esda.Moran(y, w, transformation="r", permutations=0)

print("edge_m | n | I | E[I] | Var_norm | z_norm | z_rand")
for edge_m, mi in results.items():
    print(f"{edge_m} | {mi.n} | {mi.I:.6f} | {mi.EI:.6f} | {mi.VI_norm:.6e} | {mi.z_norm:.6f} | {mi.z_rand:.6f}")
    -- result: s4-moran-pysal · python · sha256 30d9b36db807
edge_m | n | I | E[I] | Var_norm | z_norm | z_rand
5000 | 1094 | 0.273296 | -0.000915 | 3.263991e-04 | 15.177858 | 15.246590
10000 | 296 | 0.203588 | -0.003390 | 1.259792e-03 | 5.831416 | 5.847466

    In s4.py and s4.ipynb, with the blocks it runs after.

    Every printed digit equals the SQL's. z_rand is a little larger than z_norm (15.25 against 15.18 at 5 km): the randomization variance uses the values' own kurtosis (how heavy their tails are) where the normality one assumes a normal distribution's.

  3. Equal printed digits could still hide different neighbours. This block checks the pairs themselves: per grid, it builds PySAL's distance band from the same centroids, compares its pairs with s4_neighbours' one by one, and stops with an error naming the grid and the first five differing pairs if one set has a pair the other lacks. Then it prints the statistics at full precision (repr gives every digit a float holds), which the course's own test compares with the SQL's, within 5 × 10⁻⁷, every time it runs:

    s4-agree · Python

    Runs after: py-connect

    import sys

import esda
import numpy as np
from libpysal.weights import DistanceBand

with connect() as conn:
    edges = [e for (e,) in conn.execute("SELECT DISTINCT edge_m FROM s4_cells ORDER BY edge_m")]
    grids = {}
    for edge_m in edges:
        rows = conn.execute(
            "SELECT id, y, cx, cy FROM s4_cells WHERE edge_m = %s ORDER BY id", (edge_m,)
        ).fetchall()
        sql_pairs = set(conn.execute("SELECT i, j FROM s4_neighbours WHERE edge_m = %s", (edge_m,)).fetchall())
        grids[edge_m] = (rows, sql_pairs)

table = []
for edge_m, (rows, sql_pairs) in grids.items():
    w = DistanceBand(np.array([(r[2], r[3]) for r in rows]), threshold=2 * edge_m, binary=True, ids=[r[0] for r in rows])
    pysal_pairs = {(i, j) for i, js in w.neighbors.items() for j in js}
    differ = sorted(sql_pairs ^ pysal_pairs)
    if differ:
        shown = ", ".join(f"{p} only in {'SQL' if p in sql_pairs else 'PySAL'}" for p in differ[:5])
        sys.exit(f"{edge_m} m grid: {len(differ)} neighbour pairs differ; first: {shown}")
    mi = esda.Moran(np.array([r[1] for r in rows]), w, transformation="r", permutations=0)
    table.append((edge_m, mi.n, len(pysal_pairs), mi.I, mi.EI, mi.VI_norm, mi.z_norm))

print("neighbour pairs equal s4_neighbours on every grid")
print("edge_m | n | pairs | I | E[I] | Var_norm | z_norm")
for edge_m, n, pairs, *stats in table:
    print(" | ".join([str(edge_m), str(n), str(pairs), *(repr(float(s)) for s in stats)]))
    -- result: s4-agree · python · sha256 44e33d16750b
neighbour pairs equal s4_neighbours on every grid
edge_m | n | pairs | I | E[I] | Var_norm | z_norm
5000 | 1094 | 6246 | 0.2732961251160256 | -0.0009149130832570906 | 0.00032639905787023167 | 15.177857980448268
10000 | 296 | 1620 | 0.2035877773186628 | -0.003389830508474576 | 0.0012597916525627275 | 5.831415676740648

    In s4.py and s4.ipynb, with the blocks it runs after.

    Both grids pass: the 6,246 pairs at 5 km and the 1,620 at 10 km are the same pairs in both programs (each counted once from each side, so a pair of cells is two rows of s4_neighbours). The statistics differ from the SQL's only in their last digits, by about 10⁻¹² at most.

  4. The permutation test on the 5 km cells. esda's global Moran takes no seed, so NumPy's is set first; 999 shuffles:

    s4-permutations · Python

    Runs after: py-connect

    import esda
import numpy as np
from libpysal.weights import DistanceBand

with connect() as conn:
    rows = conn.execute("SELECT id, y, cx, cy FROM s4_cells WHERE edge_m = 5000 ORDER BY id").fetchall()
ids = [r[0] for r in rows]
y = np.array([r[1] for r in rows])
w = DistanceBand(np.array([(r[2], r[3]) for r in rows]), threshold=10000, binary=True, ids=ids)

np.random.seed(2026)
mi = esda.Moran(y, w, transformation="r", permutations=999)
print(f"I: {mi.I:.6f}  shuffles: {mi.permutations}")
print(f"shuffled I: mean {mi.EI_sim:.6f}  sd {mi.seI_sim:.6f}  max {mi.sim.max():.6f}")
print(f"pseudo p: {mi.p_sim:.4f}  z_sim: {mi.z_sim:.2f}")
    -- result: s4-permutations · python · sha256 e4f919f6c5a8
I: 0.273296  shuffles: 999
shuffled I: mean -0.000776  sd 0.018192  max 0.065027
pseudo p: 0.0010  z_sim: 15.07

    In s4.py and s4.ipynb, with the blocks it runs after.

    In 999 shuffles no I came near the real one: the largest was 0.0650, against 0.2733. So the pseudo p-value is the smallest possible, 0.001, and the shuffles' spread (mean −0.0008, standard deviation 0.0182) gives a z of 15.07, close to the formulas'.

Result

The grid checks pass at both sizes (the nearest centroids 8,660.25 m and 17,320.51 m apart, the next class 15,000 and 30,000 m, no duplicate centroid, no island, no pair near the threshold). Moran's I is 0.2733 on the 1,094 cells of 5 km (E[I] −0.0009, z 15.18) and 0.2036 on the 296 cells of 10 km (z 5.83): SQL and PySAL agree to every printed digit, on the same neighbour pairs. ArcGIS Pro's z would be z_rand, 15.25. In 999 shuffles no I came near the real one (the largest 0.0650), so the pseudo p-value is the smallest possible, 0.001, and the shuffles' spread gives z 15.07, close to the formulas'. Larger cells give a smaller I and a much smaller z: the modifiable areal unit problem (Grids, g6) reaches the statistic. And since Bias in records, s1, showed the records are mostly one atlas and one museum's samples, the honest sentence is "devil records cluster strongly at both scales", not "devils cluster".

Exercise

y = ln(1 + n) still lets a cell with 256 records weigh more than one with 1. Ask a plainer question: do cells with any record at all sit together? Compute Moran's I, E[I] and the z-score under normality in SQL for y = 1 when a cell has a record and 0 when it has none, on both grids, with the same row-standardized weights, and show how many cells have a record. Is the pattern weaker or stronger than for ln(1 + n), and what does that say?

Show answer

s4-ex · PostgreSQL with PostGIS 

WITH v AS (
  SELECT id, edge_m, CASE WHEN n > 0 THEN 1.0 ELSE 0.0 END::float8 AS y FROM s4_cells
), z AS (
  SELECT id, edge_m, y - avg(y) OVER (PARTITION BY edge_m) AS z FROM v
), w AS (
  SELECT edge_m, i, j, 1 / count(*) OVER (PARTITION BY i)::float8 AS w FROM s4_neighbours
), pairs AS (
  SELECT w.edge_m, w.i, w.j, w.w AS w_ij, t.w AS w_ji
  FROM w JOIN w AS t ON t.i = w.j AND t.j = w.i
), sums AS (
  SELECT p.edge_m,
         sum(p.w_ij * zi.z * zj.z) AS wzz,
         sum(p.w_ij) AS s0,
         sum((p.w_ij + p.w_ji) ^ 2) / 2 AS s1
  FROM pairs AS p JOIN z AS zi ON zi.id = p.i JOIN z AS zj ON zj.id = p.j
  GROUP BY p.edge_m
), s2 AS (
  SELECT edge_m, sum(ws ^ 2) AS s2
  FROM (SELECT edge_m, i, sum(w_ij) + sum(w_ji) AS ws FROM pairs GROUP BY edge_m, i) AS cell
  GROUP BY edge_m
), cells AS (
  SELECT z.edge_m, count(*)::float8 AS n, sum(z.z * z.z) AS zz, sum(v.y)::integer AS with_records
  FROM z JOIN v USING (id)
  GROUP BY z.edge_m
), moran AS (
  SELECT c.edge_m, c.n, c.with_records,
         c.n / s.s0 * s.wzz / c.zz AS i,
         -1 / (c.n - 1) AS ei,
         (c.n ^ 2 * s.s1 - c.n * s2.s2 + 3 * s.s0 ^ 2) / ((c.n - 1) * (c.n + 1) * s.s0 ^ 2) - (1 / (c.n - 1)) ^ 2 AS var_norm
  FROM cells AS c JOIN sums AS s USING (edge_m) JOIN s2 USING (edge_m)
)
SELECT edge_m, n::integer AS n, with_records,
       round(i::numeric, 6) AS i,
       round(ei::numeric, 6) AS ei,
       round(((i - ei) / sqrt(var_norm))::numeric, 2) AS z
FROM moran
ORDER BY edge_m;
-- result: s4-ex · 2 rows · sha256 1840f4611e5b
edge_m | n | with_records | i | ei | z
5000 | 1094 | 307 | 0.208633 | -0.000915 | 11.60
10000 | 296 | 178 | 0.229060 | -0.003390 | 6.55

In s4.sql, with the blocks it runs after.

Check

The answer: s4-moran-sql with one change, the value in the first CTE.

The answer has 307 of the 1,094 cells of 5 km with a record, and 178 of the 296 of 10 km. At 5 km, presence clusters less than the counts did (I 0.2086 against 0.2733, z 11.60 against 15.18): part of what ln(1 + n) measured was cells with many records next to each other, and with presence alone a cell of 256 records counts as much as a cell of 1. It still clusters strongly: cells with a record sit beside cells with a record. At 10 km it goes the other way (I 0.2291 against 0.2036, z 6.55 against 5.83): most cells there have a record, and what clusters is the contrast between areas where someone recorded a devil and areas where nobody did. Both versions answer a question about records: where people looked and published, not where devils are. The choice of value, like the choice of cell size, changes the number, so a Moran's I is reported with the value it measures, its cells and its weights.

Dialects

Lab schema: the PostGIS tables

critwire Learn lab database, PostGIS schema learn_curriculum

Every column, in words · The diagram alone