Generating Categorical Attributes with Spatial Autocorrelation

Building type is drawn from the right marginal distribution — sixty percent residential, thirty commercial, ten industrial — and the map looks like static. Real cities have residential districts, not residential confetti.

Part of Attribute Correlation Modeling: that page covers correlations between attributes at a location. This one is correlation of a single attribute across locations, which needs a different mechanism and is usually the more visible of the two when it is missing.

Root Cause: The Marginal Is Right and the Joint Is Absent

Drawing each feature’s category independently from a marginal distribution reproduces that distribution exactly. It also guarantees that neighbouring features are independent, which is the one property the real data does not have.

The gap shows up in a statistic almost nobody computes on synthetic output and everybody notices by eye: the join count, the number of adjacent pairs sharing a category. Under independent draws with category proportion p, the expected share of same-category joins is simply Σpᵢ², and in real spatial data it is far higher — often twice that, sometimes more. The generator is not slightly off; it is producing a structurally different object that happens to have the right histogram.

The consequence is not cosmetic. Anything downstream that aggregates over neighbourhoods — catchment analysis, district-level summaries, any model with a spatial lag term — will behave differently on confetti than on clustered data, and the difference is systematic rather than noisy.

Clustering ratio for three category generators against the observed reference Three generators run along the horizontal axis, all producing the same marginal distribution of sixty percent residential, thirty commercial and ten industrial. The vertical axis is the clustering ratio: the share of adjacent pairs sharing a category, divided by what that share would be if every draw were independent. A value of one means no spatial structure whatsoever. Independent draws sit at essentially exactly one, which is the definition rather than a coincidence — they reproduce the histogram and nothing else. A latent Gaussian field cut at quantiles with a short correlation range lifts the ratio partway. The same construction with a longer range reaches the observed reference line drawn across the chart. The note underneath makes the point that the marginal distribution is identical in all three cases, so any check that looks only at category proportions passes all of them equally, and the join count is the statistic that separates them. Identical marginals; only the join count tells them apart 0.0× 1.0× 2.0× clustering ratio (1 = independent) observed data — 1.42× 0.98× independent draws 1.62× latent field, range 1.2 1.90× latent field, range 3.0 26×26 cells, rook adjacency, marginal fixed at 60/30/10 percent in every case, fixed seed. Any check that looks only at category proportions passes all three equally. The join count is the statistic that separates a map from confetti, and it costs one pass over the adjacency list.
Same marginal distribution, three generators: independent draws sit at the theoretical floor while the observed data clusters far above it.

Prerequisite Check: Measure the Target Before Generating Anything

python
def join_counts(labels: dict, adjacency: dict) -> dict:
    """Share of adjacent pairs sharing a category, and per-category breakdown."""
    same = 0
    total = 0
    per_cat: dict = {}
    for a, neighbours in adjacency.items():
        for b in neighbours:
            if a >= b:
                continue                    # count each undirected pair once
            total += 1
            if labels[a] == labels[b]:
                same += 1
                per_cat[labels[a]] = per_cat.get(labels[a], 0) + 1
    return {"same_share": same / total, "per_category": per_cat, "pairs": total}


def independence_baseline(proportions: dict) -> float:
    """What the same-category share would be under independent draws."""
    return sum(p * p for p in proportions.values())

Run both on the reference data. The ratio between them is the calibration target, and it is the single number the generator has to hit — a same-share of 0.58 against a baseline of 0.46 is a clustering ratio of 1.26, and that is what the generated output should reproduce.

Per-category breakdown matters too. Industrial land clusters far more strongly than retail in most cities, and a generator tuned to a single global ratio will over-cluster the diffuse categories while under-clustering the concentrated ones.

Fix: Generate a Correlated Latent Field, Then Threshold It

The reliable construction is indirect. Rather than drawing categories, draw a smooth continuous field with the right spatial correlation and cut it into categories at quantiles chosen to hit the marginal.

1 — Build the correlated latent field

python
def latent_field(coords, range_m: float, rng) -> list[float]:
    """A Gaussian field smoothed to a target correlation range."""
    raw = [rng.normal() for _ in coords]
    out = []
    for i, (x, y) in enumerate(coords):
        num = den = 0.0
        for j, (u, v) in enumerate(coords):
            d2 = (x - u) ** 2 + (y - v) ** 2
            if d2 > (3 * range_m) ** 2:
                continue                      # the kernel is negligible past 3 ranges
            w = math.exp(-d2 / (2 * range_m * range_m))
            num += w * raw[j]
            den += w
        out.append(num / math.sqrt(den))       # keep unit variance
    return out

The range_m parameter is the whole control surface: it is the distance over which categories stay similar, it maps directly onto the observed clustering ratio, and it is the only thing that needs calibrating.

2 — Threshold at quantiles that reproduce the marginal

python
def assign_categories(field: list[float], proportions: dict) -> list[str]:
    """Cut the latent field at quantiles so the marginal comes out exactly right."""
    order = sorted(range(len(field)), key=lambda i: field[i])
    labels = [None] * len(field)
    cursor = 0
    for cat, share in proportions.items():
        take = round(share * len(field))
        for i in order[cursor:cursor + take]:
            labels[i] = cat
        cursor += take
    for i in order[cursor:]:                   # rounding remainder
        labels[i] = list(proportions)[-1]
    return labels

Cutting by rank rather than by a fixed threshold makes the marginal exact by construction, which means the calibration loop only has one free parameter left. That separation — marginal exact, clustering tuned — is what keeps the fit tractable.

3 — Calibrate the range against the observed ratio

python
def calibrate_range(coords, adjacency, proportions, target_ratio, rng,
                    lo=50.0, hi=5000.0, tol=0.01) -> float:
    """Bisect on the correlation range until the join-count ratio matches."""
    baseline = independence_baseline(proportions)
    for _ in range(24):
        mid = math.sqrt(lo * hi)               # geometric bisection: range spans decades
        labels = assign_categories(latent_field(coords, mid, rng), proportions)
        ratio = join_counts(dict(enumerate(labels)), adjacency)["same_share"] / baseline
        if abs(ratio - target_ratio) < tol:
            return mid
        lo, hi = (mid, hi) if ratio < target_ratio else (lo, mid)
    return math.sqrt(lo * hi)

Geometric rather than arithmetic bisection, because a plausible correlation range spans two or three orders of magnitude and arithmetic bisection wastes most of its iterations at the top of the interval.

Calibrating the latent correlation range against a target clustering ratio The horizontal axis is the correlation range of the latent Gaussian field, in cells, on a scale spanning an order of magnitude. The vertical axis is the resulting clustering ratio measured on the thresholded categories. The curve starts just above one at the shortest range, where the field is essentially white noise and the output is indistinguishable from independent draws, rises steeply through the middle of the range, and begins to flatten at the long end as the map turns into a few large blocks. A horizontal line marks the target ratio taken from the observed data. Four points mark successive geometric-bisection probes, each labelled with its iteration number, closing on the range where the curve crosses the target. The note underneath explains the choice of geometric rather than arithmetic bisection and warns about the flat end of the curve, where a small change in target implies a large change in range and the fit becomes ill-conditioned. One free parameter: the marginal is exact, the range does the fitting 1 2 3 4 5 0.50 1.00 1.50 2.00 2.50 latent correlation range (cells) clustering ratio target 1.42× probe 1 probe 2 probe 3 probe 4 22×22 cells, one shared noise draw across all ranges so the curve is a function of the range alone, fixed seed. Bisection is geometric because a plausible range spans decades and arithmetic bisection wastes most of its iterations at the top of the interval. Watch the flat end: where the curve levels off, a small change in the target implies a large change in the range, and the fit stops being well determined.
Clustering ratio as a function of the latent correlation range, with the bisection path to a target ratio marked.

Verification Step: Check Both Halves Separately

python
def test_marginal_is_exact(labels, proportions, tol=0.005):
    got = tally_shares(labels)
    for cat, want in proportions.items():
        assert abs(got[cat] - want) < tol, f"{cat}: {got[cat]:.3f} vs {want:.3f}"


def test_clustering_matches_target(labels, adjacency, proportions, target, tol=0.05):
    ratio = (join_counts(labels, adjacency)["same_share"]
             / independence_baseline(proportions))
    assert abs(ratio - target) < tol, f"clustering ratio {ratio:.3f} vs target {target:.3f}"


def test_per_category_clustering(labels, adjacency, observed_per_cat, tol=0.1):
    """A global match can hide two categories that are wrong in opposite directions."""
    got = join_counts(labels, adjacency)["per_category"]
    for cat, want in observed_per_cat.items():
        assert abs(got[cat] - want) / want < tol, f"{cat} clusters wrongly"


def test_not_degenerate(labels, adjacency):
    """Too much clustering is also a failure: three giant blobs is not a city."""
    assert count_connected_components(labels, adjacency) > 20

The last test is the one people leave out. Calibration pressure runs in one direction — the symptom was too little clustering — and it is entirely possible to overshoot into a map with three enormous single-category regions, which reproduces the join-count target and looks nothing like the reference.

Edge Cases & Gotchas

Anisotropic clustering. Land use follows corridors — along rivers, along transport lines — so an isotropic kernel under-clusters along the corridor and over-clusters across it. Replacing the distance term with an anisotropic metric is a two-line change and is worth it wherever the underlying geography has strong directionality.

One map, three adjacency definitions, three clustering ratios Three adjacency definitions run along the horizontal axis, all applied to exactly the same generated map with exactly the same categories in exactly the same places. The vertical axis is the clustering ratio, the share of adjacent pairs sharing a category divided by the independence baseline. Rook contiguity, counting only the four cells sharing an edge, gives the highest value. Queen contiguity, which adds the four diagonal neighbours, gives a lower one, because diagonal neighbours are further apart and less likely to match. A twelve-neighbour disc reaching two cells out gives a lower value again, for the same reason extended further. Nothing about the map changed between the three bars; only the question being asked changed. The note underneath draws the operational conclusion: the adjacency definition is part of the target, not an implementation detail, and a generated release calibrated under one definition and validated under another will appear to miss a target it actually hit. The same map, three definitions of 'adjacent', three answers 0.00× 1.00× 2.00× clustering ratio 1.87× rook (4) 1.81× queen (8) 1.74× k = 12 disc One 24×24 map, identical categories in identical places, measured three ways with a fixed seed. The adjacency definition is part of the target, not an implementation detail. A release calibrated under one definition and validated under another will appear to miss a target it actually hit — and the discrepancy is large enough to send somebody looking for a modelling bug that is not there.
The same generated map measured three ways: the adjacency definition changes the number by more than most calibration tolerances allow.

Categories with hard spatial constraints. Industrial zoning that cannot occur inside a protected area is a constraint, not a tendency, and it belongs as a mask applied to the ranking rather than as pressure on the latent field. Trying to express a hard constraint through correlation produces a generator that violates it occasionally, which is the worst outcome.

The adjacency definition changes the number. Join counts computed over queen contiguity, rook contiguity or k-nearest neighbours give different values for the same map. The generated and reference figures have to use the same definition, and that definition belongs in the contract alongside the ratio itself.

Very large feature counts. The naïve field construction is quadratic. Past a few tens of thousands of features it needs a spatial index or a grid-based approximation — the same chunking argument as everywhere else, with the extra requirement that the kernel support must not exceed the chunk overlap or the correlation will break at chunk boundaries.

Why Not a Potts Model

The other standard answer is a Markov random field — a Potts model, sampled with Gibbs or Swendsen-Wang, where the energy penalises adjacent cells with different categories. It is the more principled construction and it is usually the wrong practical choice here.

The marginal is not controllable. In the latent-field construction the marginal is exact by construction, because the categories are assigned by rank. In a Potts model the marginal emerges from the interaction strength and the external field together, so hitting a target proportion means fitting a second parameter, and the two parameters interact — raising the interaction strength to increase clustering also shifts the proportions, and the fit becomes two-dimensional for no gain.

Convergence is not observable. A Gibbs sampler on a Potts model has to reach equilibrium before the output means anything, and near the critical interaction strength — which is exactly where realistic clustering levels sit — mixing is slow and there is no reliable stopping rule. A run that has not converged produces plausible-looking output, which is the worst property a generator can have.

Reproducibility is harder. The latent-field construction is a fixed sequence of arithmetic from a seed. A Gibbs sampler’s output depends on the iteration count, the sweep order and the initialisation, all of which have to be pinned and recorded to get a byte-identical release.

The Potts model earns its place when the interaction itself is the object of study, or when the category structure has genuinely local rules that a smooth latent field cannot express — a category that can only occur adjacent to another, for instance. For reproducing an observed clustering level at a target marginal, the latent field is less elegant, more controllable, and finishes in one pass.