Computing Moran’s I on Masked Spatial Data

Compute global Moran’s I on the true grid and on the masked grid using the same spatial weights matrix $W$, then confirm the relative deviation is under 10% — evidence that masking preserved the spatial autocorrelation that makes the map analytically useful.

Core Calculation # permalink

Global Moran’s I summarises spatial autocorrelation across $n$ cells with values $x_i$, mean $\bar{x}$, and a spatial weights matrix $W = [w_{ij}]$ encoding which cells are neighbours:

The right-hand fraction is a weighted covariance of each cell against its neighbours, normalised by the total variance. $I$ ranges roughly from $-1$ (perfect dispersion, a checkerboard) through $\approx 0$ (no autocorrelation, spatial randomness) to $+1$ (strong clustering, neighbours alike). Under the null hypothesis of no autocorrelation, the expected value is:

which is why independent noise — the kind injected by coordinate jittering and noise-injection methods or a differential-privacy mechanism — drives measured $I$ toward zero: the noise itself carries no spatial structure.

Utility check # permalink

Let $I_{\text{true}}$ and $I_{\text{mask}}$ be Moran’s I on the original and masked grids under identical $W$. The preservation criterion is the relative deviation:

Worked small example # permalink

Take a $2\times2$ grid with rook contiguity (edge neighbours only). Cell layout and true values:

Neighbour pairs (rook): (1,2), (1,3), (2,4), (3,4) and their reverses, so $\sum\sum w_{ij} = 8$. Mean $\bar{x} = 2.5$; deviations $z = (1.5, 0.5, -0.5, -1.5)$. Denominator $\sum z_i^2 = 2.25 + 0.25 + 0.25 + 2.25 = 5.0$.

Weighted cross-product (each undirected pair counted twice): pair (1,2): $1.5 \times 0.5 = 0.75$; (1,3): $1.5 \times -0.5 = -0.75$; (2,4): $0.5 \times -1.5 = -0.75$; (3,4): $-0.5 \times -1.5 = 0.75$. Sum over directed pairs $= 2 \times (0.75 - 0.75 - 0.75 + 0.75) = 0$.

This tiny grid is deliberately non-autocorrelated (high next to low). Real heatmaps with genuine hotspots yield strongly positive $I$; the code below demonstrates preservation on such a grid.

Python Implementation # permalink

The function is numpy-only so it runs without geospatial dependencies; a note shows the esda / libpysal equivalent. Grids are assumed to be regular rasters in a projected CRS (EPSG:3857 web-Mercator) so that row/column adjacency is a valid neighbour relation.

import numpy as np
from numpy.typing import NDArray


def rook_weights(rows: int, cols: int) -> NDArray[np.float64]:
    """
    Build a binary rook-contiguity weights matrix for an rows x cols grid.

    Cells are flattened in row-major (C) order. Rook contiguity links cells
    sharing an edge (up/down/left/right). Use the SAME matrix on the true and
    masked grids so the Moran's I comparison is apples-to-apples.
    """
    n = rows * cols
    w = np.zeros((n, n), dtype=float)
    for r in range(rows):
        for c in range(cols):
            i = r * cols + c
            if c + 1 < cols:               # right neighbour
                j = i + 1
                w[i, j] = w[j, i] = 1.0
            if r + 1 < rows:               # down neighbour
                j = i + cols
                w[i, j] = w[j, i] = 1.0
    return w


def morans_i(values: NDArray[np.float64], w: NDArray[np.float64]) -> float:
    """
    Global Moran's I for a flattened grid under weights matrix w.

    Args:
        values: Flattened cell values (row-major), same order used to build w.
        w: Spatial weights matrix, shape (n, n).

    Returns:
        Global Moran's I in roughly [-1, 1]. Near 0 means spatial randomness.
    """
    x = np.asarray(values, dtype=float).ravel()
    n = x.size
    z = x - x.mean()                       # centre by the grid mean
    w_sum = w.sum()
    if w_sum == 0 or np.all(z == 0):
        return 0.0
    # Weighted cross-product of neighbour deviations: z^T W z
    numerator = z @ w @ z
    denominator = np.sum(z ** 2)
    return (n / w_sum) * (numerator / denominator)


# --- Example: 20x20 clustered grid, masked with independent noise ---
rng = np.random.default_rng(0)
rows = cols = 20
# Build an autocorrelated field: a smooth gradient + mild local noise.
yy, xx = np.mgrid[0:rows, 0:cols]
true_grid = (50 + 4 * xx + 4 * yy + rng.normal(0, 3, (rows, cols))).ravel()

w = rook_weights(rows, cols)
i_true = morans_i(true_grid, w)

# Masked grid: independent per-cell noise (e.g. a DP or jittering step).
masked_grid = true_grid + rng.normal(0, 8.0, size=true_grid.size)
i_mask = morans_i(masked_grid, w)

print(f"I_true = {i_true:.4f}, I_mask = {i_mask:.4f}")

# esda/libpysal equivalent (optional):
#   from libpysal.weights import lat2W
#   from esda.moran import Moran
#   w_pysal = lat2W(rows, cols, rook=True)
#   Moran(true_grid, w_pysal).I

Key decisions: the weights matrix is built once and reused for both grids so the comparison is valid; values are centred by the grid mean inside morans_i; and the z @ w @ z quadratic form is the exact vectorised numerator of the Moran’s I definition, avoiding double loops.

Verification # permalink

Confirm the masked map preserved spatial autocorrelation within the 10% band:

def verify_morans_preservation(
    i_true: float, i_mask: float, max_reldev: float = 0.10
) -> dict:
    """Check relative deviation of masked Moran's I against the true value."""
    if i_true == 0:
        raise ValueError("I_true is ~0; grid has no autocorrelation to preserve.")
    reldev = abs(i_mask - i_true) / abs(i_true)
    return {
        "i_true": round(i_true, 4),
        "i_mask": round(i_mask, 4),
        "relative_deviation": round(reldev, 4),
        "passed": bool(reldev <= max_reldev),
    }

print(verify_morans_preservation(i_true, i_mask))
# e.g. {'i_true': 0.83, 'i_mask': 0.78, 'relative_deviation': 0.06, 'passed': True}

• Same weights matrix used on both grids.
• $I_{\text{true}}$ $I_{\text{true}}$ meaningfully positive (there is structure to preserve).
• Relative deviation $\le 0.10$ Relative deviation $\le 0.10$ → autocorrelation preserved.
• If deviation $> 0.20$ If deviation $> 0.20$, relax the mask or coarsen the grid.

Edge Cases & Adjustments # permalink

• Edge effects. Corner and border cells have fewer neighbours, biasing $I$ on small grids. Prefer row-standardised weights (each row of $W$ divided by its neighbour count) so under-connected edge cells are not systematically down-weighted, and interpret small grids cautiously.
• Weight choice changes the value. Rook vs. queen contiguity, distance-band, or k-nearest-neighbour weights all yield different $I$. This is fine for a before/after comparison as long as the identical $W$ is applied to both grids; never compare an $I$ computed under different weight schemes.
• Sparse / mostly-empty grids. When most cells are zero, the variance denominator is dominated by a few hotspots and $I$ becomes unstable. Aggregate to a coarser resolution before computing, or restrict the analysis to the populated extent.
• CRS gotcha. Contiguity assumes a regular grid in a projected CRS. Computing neighbour relations directly on WGS84 (EPSG:4326) degrees distorts adjacency near the poles; reproject to a metric CRS such as EPSG:3857 or a local UTM zone before building $W$.

FAQ # permalink

What does Moran’s I actually tell you about a masked map?

It measures spatial autocorrelation — how similar neighbouring cells are. On a masked or noised map it reveals whether the clustering that gives the map its value survived the privacy transformation. A value near the original means neighbours remain alike; a value collapsing toward zero means the mask scattered the pattern into spatial noise.

Why does independent noise drive Moran’s I toward zero?

Purely independent noise has an expected Moran’s I of $-1/(n-1)$, essentially zero for a large grid. Adding independent per-cell noise blends this near-zero signal into the true autocorrelated signal, diluting the measured $I$ downward in proportion to the noise-to-signal variance ratio.

Rook or queen contiguity for the weights matrix?

Rook links edge-sharing cells (4 neighbours); queen also links diagonals (8 neighbours). Queen is usually preferred for raster heatmaps because it captures diagonal clustering, but what matters most is applying the exact same $W$ to the true and masked grids.

What relative deviation in Moran’s I is acceptable?

A relative deviation under 10% is the common preservation target; within it, spatial clustering is considered intact for uses such as hotspot detection. Beyond about 20%, treat the masked map as having materially altered spatial structure and relax the privacy parameters or coarsen the grid.

Related

• Utility-Preservation Metrics for Masked Maps — parent guide to the full metric family
• Measuring Utility Loss in Differentially Private Heatmaps — the wider metric suite that includes this check
• Accuracy vs. Utility Tradeoffs in Geospatial DP — framing the privacy-utility balance Moran’s I helps quantify
• Coordinate Jittering & Noise-Injection Methods — the masking step whose spatial impact this metric measures
• Differential Privacy for Location Data — the mechanisms that add the independent noise being evaluated

← Back to Utility-Preservation Metrics for Masked Maps