critwire · Learn · Learn's data

Hot and cold spots s5

Prerequisites: Global Moran's I (s4)

Where do values cluster, and how many of those places survive a careful test? Getis-Ord Gi* in plain SQL and in PySAL, local Moran's I with permutations, its quadrants in SQL, and the false discovery rate, on the Tasmanian devil's records counted into hexagons. Everything runs in your own PostgreSQL with PostGIS and the Python environment, on the lab tables (the setup page).

Concept

Global Moran's I, s4, gave one number for the whole island; this module asks where values cluster, and how many of those places survive a careful test.

Worked example

Gi* in plain SQL. Like s2 and s3, the blocks read s4's s4_cells and s4_neighbours: every download of this module runs s4-cells first, right after connecting, so s5.sql and s5.py work on their own; to run the blocks one at a time, run s4-cells once first.

The block keeps its result as a temporary view, s5_gistar, one row per 5 km cell: id, z, the two-sided p, and fdr, true where the cell survives the false discovery rate. The exercise below reads it, and so does the capstone (c1), instead of repeating the query. A temporary view lasts until the session ends:

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.

s5-gistar · PostgreSQL with PostGIS 

CREATE OR REPLACE TEMP VIEW s5_gistar AS
WITH c AS (
  SELECT id, y FROM s4_cells WHERE edge_m = 5000
), stats AS (
  SELECT count(*)::float8 AS n, avg(y) AS ybar, sqrt(avg(y * y) - avg(y) ^ 2) AS s FROM c
), w AS (
  SELECT i, j FROM s4_neighbours WHERE edge_m = 5000
  UNION ALL
  SELECT id, id FROM c
), sums AS (
  SELECT w.i, sum(c.y) AS wy, count(*)::float8 AS wsum
  FROM w JOIN c ON c.id = w.j
  GROUP BY w.i
), z AS (
  SELECT u.i, (u.wy - st.ybar * u.wsum) / (st.s * sqrt((st.n * u.wsum - u.wsum ^ 2) / (st.n - 1))) AS z
  FROM sums AS u, stats AS st
), p AS (
  SELECT i, z, erfc(abs(z) / sqrt(2)) AS p,
         row_number() OVER (ORDER BY erfc(abs(z) / sqrt(2)), i) AS rank,
         count(*) OVER () AS m
  FROM z
), cut AS (
  SELECT coalesce(max(rank) FILTER (WHERE p <= rank * 0.05 / m), 0) AS k FROM p
)
SELECT p.i AS id, p.z, p.p, p.rank <= cut.k AS fdr
FROM p, cut;
SELECT count(*) AS cells,
       count(*) FILTER (WHERE z > 0 AND p <= 0.05) AS hot_05,
       count(*) FILTER (WHERE z < 0 AND p <= 0.05) AS cold_05,
       count(*) FILTER (WHERE z > 0 AND fdr) AS hot_fdr,
       count(*) FILTER (WHERE z < 0 AND fdr) AS cold_fdr,
       round(max(z)::numeric, 6) AS max_z,
       round(min(z)::numeric, 6) AS min_z
FROM s5_gistar;
-- result: s5-gistar · 1 row · sha256 84ede64c9841
cells | hot_05 | cold_05 | hot_fdr | cold_fdr | max_z | min_z
1094 | 99 | 0 | 39 | 0 | 10.632449 | -1.289294

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

Read it one CTE at a time. stats makes n, ȳ and S, the population standard deviation (sqrt(avg(y²) − ȳ²), dividing by n, not n − 1, as the formula requires). w adds each cell to its own neighbours (the star) with UNION ALL, so sums gets Σ_j w_ij y_j and, with binary weights, Σ_j w_ij = Σ_j w_ij² = the number of cells in the neighbourhood. z is the Gi* formula. p turns z into a two-sided p-value with erfc (new in PostgreSQL 16) and ranks the cells by it, i breaking ties. cut is the Benjamini-Hochberg step: the largest rank whose p-value is at most rank · 0.05 / m (0 when no cell passes). The last query counts the cells.

99 cells are hot at p ≤ 0.05 and 39 survive the false discovery rate; the highest z is 10.63. There is no cold spot at all, and the lowest z is only −1.29: 787 of the 1,094 cells hold no record (s4's exercise), so an empty neighbourhood is the typical one and can't stand out below the mean.

Lab

  1. The same Gi* in PySAL: G_Local with star=True and binary weights (transform="B"), from the same cells and the same band. esda's own p_norm is one-sided, so the block computes the two-sided p-value from Zs; esda's fdr returns the p-value cut:

    s5-gistar-pysal · Python

    Runs after: py-connect

    import esda
import numpy as np
from libpysal.weights import DistanceBand
from scipy.stats import norm

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)

g = esda.G_Local(y, w, transform="B", star=True, permutations=0)
p = 2 * norm.sf(np.abs(g.Zs))
cut = esda.fdr(p, 0.05)
print("cells | hot_05 | cold_05 | hot_fdr | cold_fdr | max_z | min_z")
print(f"{len(y)} | {np.sum((g.Zs > 0) & (p <= 0.05))} | {np.sum((g.Zs < 0) & (p <= 0.05))} | "
      f"{np.sum((g.Zs > 0) & (p <= cut))} | {np.sum((g.Zs < 0) & (p <= cut))} | {g.Zs.max():.6f} | {g.Zs.min():.6f}")
    -- result: s5-gistar-pysal · python · sha256 42a8f2bf8f68
cells | hot_05 | cold_05 | hot_fdr | cold_fdr | max_z | min_z
1094 | 99 | 0 | 39 | 0 | 10.632449 | -1.289294

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

    The line equals the SQL's, digit for digit: the course's test compares the two every time the labs run.

  2. Local Moran's I with 999 permutations (seeded), by quadrant: all cells, significant at p ≤ 0.05, and after the false discovery rate. transformation="r" row-standardizes the weights, as in s4; seed=2026 makes the permutations repeatable:

    s5-lisa · 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)

lisa = esda.Moran_Local(y, w, transformation="r", permutations=999, seed=2026)
cut = esda.fdr(lisa.p_sim, 0.05)
print(f"false discovery rate cut: p <= {cut:.4f}")
print("quadrant | cells | p <= 0.05 | after FDR")
for q, name in [(1, "high-high"), (2, "low-high"), (3, "low-low"), (4, "high-low")]:
    in_q = lisa.q == q
    print(f"{name} | {in_q.sum()} | {np.sum(in_q & (lisa.p_sim <= 0.05))} | {np.sum(in_q & (lisa.p_sim <= cut))}")
    -- result: s5-lisa · python · sha256 ff0ae1985d8b
false discovery rate cut: p <= 0.0021
quadrant | cells | p <= 0.05 | after FDR
high-high | 185 | 68 | 16
low-high | 230 | 34 | 2
low-low | 557 | 0 | 0
high-low | 122 | 29 | 29

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

    The false discovery rate's cut is p ≤ 0.0021, so only cells with a pseudo p-value of 0.001 or 0.002 survive: the floor 999 permutations allow, and the step above it. 68 high-high cells are significant at p ≤ 0.05 and 16 after the correction. 29 high-low outliers survive too: single cells of many records among empty ones, the shape of an atlas grid point or a sampling site. No low-low cell is significant.

  3. The quadrants need no permutations: they come from the sign of each cell's deviation from the mean and of its spatial lag, both plain SQL. This block computes them from s4_cells and s4_neighbours, with the records each quadrant holds; the course's test checks that its cells per quadrant equal esda's cells column:

    s5-lisa-sql · PostgreSQL with PostGIS 

    WITH c AS (
  SELECT id, n, y - avg(y) OVER () AS z FROM s4_cells WHERE edge_m = 5000
), lag AS (
  SELECT nb.i AS id, avg(c.z) AS lag
  FROM s4_neighbours AS nb JOIN c ON c.id = nb.j
  WHERE nb.edge_m = 5000
  GROUP BY nb.i
), q AS (
  SELECT c.id, c.n,
         CASE WHEN c.z > 0 AND l.lag > 0 THEN 1
              WHEN c.z <= 0 AND l.lag > 0 THEN 2
              WHEN c.z <= 0 AND l.lag <= 0 THEN 3
              ELSE 4 END AS q
  FROM c JOIN lag AS l USING (id)
)
SELECT (ARRAY['high-high', 'low-high', 'low-low', 'high-low'])[q] AS quadrant,
       count(*) AS cells,
       sum(n) AS records
FROM q
GROUP BY q
ORDER BY q;
    -- result: s5-lisa-sql · 4 rows · sha256 61ccb6a48dc8
quadrant | cells | records
high-high | 185 | 1613
low-high | 230 | 0
low-low | 557 | 0
high-low | 122 | 359

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

    avg(c.z) over a cell's neighbours is the row-standardized lag: each of its k_i neighbours weighs 1/k_i. The quadrants are numbered as esda's q (1 high-high, 2 low-high, 3 low-low, 4 high-low), with a value exactly at the mean counted as low, as esda does. The counts equal esda's: 185, 230, 557, 122. The records column shows what "low" means here: both low quadrants hold no record at all. The mean of y is below ln 2, the value of a cell with one record, so a single record lifts a cell above the mean, "high" means "has a record" and "low" means "empty".

Result

Gi* finds 99 hot cells at p ≤ 0.05 and 39 after the false discovery rate, with the highest z 10.63, and no cold spot: SQL and PySAL give the same line. Local Moran's I puts 185 cells in high-high, 230 in low-high, 557 in low-low and 122 in high-low, in SQL as in esda; 68 high-high cells are significant at p ≤ 0.05 and 16 after the correction, with 29 high-low outliers, and no low-low cell. The correction removes more than half of Gi*'s hot spots and three in four of LISA's high-high cells. What is left are the places the atlas and the museum's teams recorded devils most: hot spots of recording effort, which a report or a map labels as such, never as where devils are most common.

Exercise

How much of the record do the 39 hot spots that survive the false discovery rate hold, and who recorded it? For each dataset with a devil record in a 5 km cell, count its records in those 39 cells and its records in any cell, each with its share of the column's total. Compare the shares: does one dataset weigh more in the hot spots than in the whole sample?

Show answer

s5-ex · PostgreSQL with PostGIS

Runs after: s5-gistar

WITH devils AS (
  SELECT o.dataset_key, 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'
), placed AS (
  SELECT d.dataset_key, g.fdr AND g.z > 0 AS hot
  FROM devils AS d
  JOIN s4_cells AS c ON c.edge_m = 5000 AND ST_Contains(c.geom, d.geom)
  JOIN s5_gistar AS g ON g.id = c.id
)
SELECT ds.title,
       count(*) FILTER (WHERE p.hot) AS hot_records,
       round(100.0 * count(*) FILTER (WHERE p.hot) / sum(count(*) FILTER (WHERE p.hot)) OVER (), 1) AS hot_pct,
       count(*) AS records_in_cells,
       round(100.0 * count(*) / sum(count(*)) OVER (), 1) AS pct
FROM placed AS p JOIN datasets AS ds USING (dataset_key)
GROUP BY ds.title
ORDER BY hot_records DESC, records_in_cells DESC, ds.title;
-- result: s5-ex · 8 rows · sha256 a2b623ef575f
title | hot_records | hot_pct | records_in_cells | pct
Tasmanian Natural Values Atlas | 870 | 78.5 | 1402 | 71.1
Tasmanian Museum and Art Gallery provider for OZCAM | 237 | 21.4 | 555 | 28.1
Northern Territory Museum and Art Gallery provider for OZCAM | 1 | 0.1 | 1 | 0.1
Museums Victoria provider for OZCAM | 0 | 0.0 | 5 | 0.3
MVZ Mammal Collection (Arctos) | 0 | 0.0 | 3 | 0.2
… 3 more rows

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

Check

The answer: each record placed in its cell as in s4-cells, then joined to s5_gistar; the block runs after s5-gistar, so the view exists.

The 39 hot spots, 3.6 % of the 1,094 cells, hold 1,108 of the 1,972 records in cells (870 + 237 + 1), 56 %. The atlas supplies 78.5 % of them against 71.1 % of all records in cells, and the museum's dataset 21.4 % against 28.1 %: the hot spots are, even more than the sample, the atlas's. Five of the eight datasets have no record in a hot spot. So the map of hot spots is largely a map of where the atlas's observers worked; before reading it as anything about devils, a report says which datasets made it, as s1 did for the whole sample. FILTER (WHERE p.hot) counts one kind of row inside a group, and sum(…) OVER () over the grouped counts gives the column's total for the shares.

Dialects

Lab schema: the PostGIS tables

critwire Learn lab database, PostGIS schema learn_curriculum

Every column, in words · The diagram alone