Geostatistical Modeling Spatial Statistics with Python

Choosing an Analysis Scale for Spatial Statistics

TL;DR: Bound the choice from three directions — the process range from a variogram (analyse at or below it), the decision unit if one is fixed by policy, and the data floor set by minimum counts and disclosure rules. Take the coarsest floor and the finest ceiling; analyse inside that window, at the scale where the sensitivity curve plateaus. Then publish the unit, the constraint that bound it, and the band.

Why This Matters

A sensitivity sweep tells you how much your answer moves; it does not tell you which answer to report. That decision is where most analyses quietly become indefensible — a unit is inherited from whatever shapefile was to hand, and the justification is written afterwards. Making the choice explicitly, from constraints that exist before the data are aggregated, is what separates a stated modelling assumption from an accident. The mechanics of the sweep itself are in measuring MAUP sensitivity across aggregation levels, and the theory in spatial scale and the modifiable areal unit problem.

Environment and Version Pinning

bash
pip install "geopandas>=1.0" "numpy>=1.24" "pandas>=2.0" "scikit-gstat>=1.0" "libpysal>=4.9"
python
import geopandas as gpd
import numpy as np
import pandas as pd
import skgstat as skg

Step-by-Step Implementation

1. Estimate the process range from point support

The variogram range is the most defensible upper bound available, because it is a property of the phenomenon rather than of the units. Fit it on the unaggregated observations — the whole point is to learn the scale before you impose one.

python
coords = np.column_stack([pts.geometry.x.to_numpy(), pts.geometry.y.to_numpy()])
values = pts["x"].to_numpy()

V = skg.Variogram(
    coords, values,
    model="spherical",
    n_lags=18,
    maxlag=0.4,             # fraction of the maximum separation distance
    normalize=False,
)
nugget, sill, rng_m = V.parameters[2], V.parameters[1], V.parameters[0]
print(f"range   = {rng_m:,.0f} m")
print(f"sill    = {sill:.3f}")
print(f"nugget  = {nugget:.3f}  ({100 * nugget / (nugget + sill):.0f}% of total)")
text
range   = 3,180 m
sill    = 0.742
nugget  = 0.061  (8% of total)

Read this as: observations more than about 3.2 km apart carry essentially no information about each other. Cells wider than that average across an entire correlation structure, which is why autocorrelation statistics fall apart above the range. The full treatment of what these three parameters mean is in estimating nugget, sill and range parameters.

The analysis window between the data floor and the process range A horizontal scale of cell sizes from 100 metres to 10 kilometres. A red zone on the left marks cells too fine to meet the minimum count. A red zone on the right marks cells coarser than the 3.2 kilometre variogram range. The green window between them, roughly 600 metres to 3.2 kilometres, is where the analysis unit should sit. A marker shows a 1 kilometre choice inside the window and a 5 kilometre administrative unit outside it. The window is what the constraints leave you 100 m600 m1 km 3.2 km5 km10 km Below the data floor too few observations per zone disclosure control bites Analysis window counts sufficient, structure preserved choose where the sweep curve plateaus 600 m – 3.2 km Beyond the process range each cell averages a whole correlation cycle neighbouring cells go independent 1 km grid chosen inside the window 5 km wards model fine, report coarse When the decision unit falls outside the window, fit inside it and aggregate the predictions afterwards

2. Compute the data floor

The lower bound is arithmetic, not judgement. Sweep cell sizes downward until the fraction of zones meeting the minimum count falls below what you can accept.

python
def coverage(points, cell, min_n=30, side=10_000.0):
    col = np.floor(points.geometry.x.to_numpy() / cell).astype(np.int64)
    row = np.floor(points.geometry.y.to_numpy() / cell).astype(np.int64)
    counts = pd.Series(1, index=pd.MultiIndex.from_arrays([col, row])).groupby(level=[0, 1]).size()
    ok = (counts >= min_n).sum()
    return {"cell_m": cell, "zones": len(counts), "zones_ok": int(ok),
            "pct_ok": 100.0 * ok / len(counts),
            "median_n": float(counts.median())}

floor = pd.DataFrame([coverage(pts, c, min_n=30) for c in (200, 400, 600, 800, 1000, 1500)])
print(floor.to_string(index=False, float_format=lambda v: f"{v:.1f}"))
text
 cell_m  zones  zones_ok  pct_ok  median_n
    200   2401       0.0     0.0       1.0
    400    625      12.0     1.9       6.0
    600    289      96.0    33.2      13.0
    800    169     147.0    87.0      23.0
   1000    100     100.0   100.0      40.0
   1500     49      49.0   100.0      81.0

With a 30-observation minimum, the floor sits between 800 m and 1 km. Anything finer is a study of the densest corner of the study area rather than of the study area.

Coverage against cell size, and where the data floor falls A bar chart plots the percentage of zones holding at least thirty observations for cell sizes of 200, 400, 600, 800, 1000 and 1500 metres. The bars read 0, 1.9, 33.2, 87.0, 100 and 100 percent, with median counts per cell of 1, 6, 13, 23, 40 and 81. The three finest bars are drawn in red as below the floor, the 800 metre bar in amber, and the two coarsest in green. A bracket marks the data floor between 800 metres and 1 kilometre. The data floor is a cliff, not a slope 10075 50250 zones reaching n ≥ 30 (%) 0.0%1.9%33.2% 87.0%100%100% Clustered sampling makes it a cliff at 200 m the median cell holds 1 point at 400 m only 1.9% of cells reach n ≥ 30 coverage then jumps 33% → 87% in one step 200 m400 m600 m 800 m1 km1.5 km cell size 1613 234081 median n per cell data floor: 800 m – 1 km Above 1 km coverage saturates — the upper bound comes from the 3.2 km variogram range, not from this curve.

3. State the decision unit, if there is one

Some units are not modifiable in any useful sense. If a funding formula, a statutory boundary or a service catchment fixes the reporting geography, that unit is a given — but it is a reporting constraint, not necessarily a modelling one.

python
DECISION_UNIT_M = 5000.0     # ward-scale reporting

window_lo, window_hi = 1000.0, rng_m
inside = window_lo <= DECISION_UNIT_M <= window_hi
print(f"analysis window : {window_lo:,.0f} m – {window_hi:,.0f} m")
print(f"decision unit   : {DECISION_UNIT_M:,.0f} m — "
      f"{'inside' if inside else 'OUTSIDE the window'}")
text
analysis window : 1,000 m – 3,180 m
decision unit   : 5,000 m — OUTSIDE the window

This is the common and awkward case, and it has one correct answer: fit inside the window, then aggregate the predictions to the decision unit. Fitting directly on 5 km wards would bake the aggregation bias into every coefficient, where no amount of later reporting can separate it out.

4. Pick the point on the curve

Inside the window, prefer the scale where the sensitivity curve flattens — the point past which extra aggregation buys stability at the cost of resolution.

python
ladder = scale_ladder(pts, pearson_r, [1000, 1250, 1500, 2000, 2500, 3000])
ladder["delta"] = ladder["value"].diff().abs()
print(ladder.to_string(index=False, float_format=lambda v: f"{v:.3f}"))
text
 cell_m  zones  dropped  value  delta
   1000    100        0  0.771    NaN
   1250     64        0  0.812  0.041
   1500     44        0  0.838  0.026
   2000     25        0  0.887  0.049
   2500     16        0  0.913  0.026
   3000     11        0  0.930  0.017

No plateau appears here — the curve is still climbing at the top of the window, which is itself the finding: this correlation is dominated by the scale effect throughout the admissible range. The right call is the finest unit in the window (1 km), because it is the least aggregated estimate that meets the data floor, with the climb reported as a caveat.

The three constraints and what each one settles Three columns compare the process range, the data floor and the decision unit. Each column lists the direction it bounds, the evidence that establishes it, and the failure that follows from ignoring it: structure destroyed by over-aggregation, noise and disclosure breaches from under-aggregation, and a result nobody can act on. Three constraints, three different kinds of evidence bounds evidence ignore it and Process range upper — cells no coarser variogram fitted on unaggregated point support range = 3.2 km here autocorrelation collapses every cell averages a full cycle Data floor lower — cells no finer coverage sweep against a minimum count per zone 1 km for n ≥ 30 here noise swamps signal and disclosure rules are breached Decision unit fixed — reporting geography policy, statute or operational catchment 5 km wards here nobody can act on it or bias is fitted into coefficients Coarsest floor ≤ unit ≤ finest ceiling. An empty window means the question needs different data, not a different unit.

5. Record the decision alongside the band

python
decision = {
    "unit_m": 1000.0,
    "bound_by": "data floor (min 30 obs/zone at 1 km)",
    "process_range_m": round(float(rng_m)),
    "window_m": [1000.0, round(float(rng_m))],
    "value_at_unit": 0.771,
    "zoning_95ci": [0.735, 0.802],
    "across_scale_range": [0.489, 0.982],
    "note": "curve still climbing at window top; finest admissible unit chosen",
}
pd.Series(decision).to_json("analysis_scale_decision.json", indent=2)

Persisting this object next to the results is what makes the choice reviewable a year later, when nobody remembers why the grid was 1 km.

Interpreting the Output

A window that is wide (floor well below the range) is comfortable: pick the plateau and move on. A window that is narrow means the data density and the process scale are barely compatible, and the result will be sensitive to small changes in either — say so explicitly. A window that is empty, where the data floor exceeds the process range, is the important case: there are not enough observations to resolve the structure at the scale the structure exists. No aggregation rescues that, and reporting a single correlation from it is misleading. Collect denser data, or switch to a point-support model that borrows strength across the whole field rather than within zones.

Critical Best Practices

Estimate the range before you aggregate, never after

A variogram fitted on zone means describes the zones, not the process, and its range is contaminated by the very aggregation you are trying to calibrate. Fit on the raw observations. If you only have pre-aggregated data, you cannot establish the ceiling from within the dataset and must borrow a range from literature or from a comparable fine-grained study — and should say so.

Fit fine, report coarse

When the decision unit is coarser than the window, resist the pull to model on it. Fit at the window scale, predict, then aggregate the predictions. The aggregation is then a documented, reversible step rather than a silent property of the coefficients.

Set the minimum count from the statistic, not by habit

Thirty is a convention, not a rule. A mean needs fewer observations per zone than a variance; a correlation needs more than either; a regression with four covariates needs considerably more. Derive the floor from the statistic being computed and note the derivation.

Re-derive the window when the data change

The floor moves whenever sample density changes — a new survey round, a different study area, a filter applied upstream. The ceiling moves whenever the phenomenon changes. A window computed once and reused across projects is an inherited assumption of exactly the kind this whole exercise exists to eliminate.

Troubleshooting

Symptom Likely cause Fix
Variogram range larger than the study area Unremoved trend inflating semivariance at long lags Detrend first — see stationarity and trend analysis
Data floor exceeds the process range Sample density too low for the structure present Densify sampling, or model on point support instead of aggregating
Coverage jumps from 2% to 100% between two cell sizes Highly clustered sampling, so most cells are empty until they are large Report coverage as a curve, and consider irregular zones sized by density
Plateau appears then reverses at coarse cells Too few zones for a stable statistic Ignore rungs with fewer than ~20 zones; the estimate there is noise
Different variables imply different windows Each variable has its own range Take the intersection of windows, or analyse each at its own scale and state that
Range estimate unstable across variogram models Weak or noisy structure near the origin Compare spherical, exponential and Matern fits; use the most conservative range

Next Steps

Once the unit is fixed, build the neighbour structure on it — spatial weight matrices covers the choices that follow, and measuring MAUP sensitivity across aggregation levels is the sweep to rerun whenever the data change.

Frequently Asked Questions

Should the analysis unit be smaller or larger than the variogram range?

Smaller. A unit at or below the range preserves the spatial structure you are trying to measure, because neighbouring units remain correlated and the autocorrelation statistics have something to detect. Once cells exceed the range, each cell averages over a full correlation cycle, neighbouring cells become effectively independent, and measures such as Moran’s I collapse toward their expected value regardless of the underlying process.

What if the decision unit is coarser than the process range?

Analyse at the finer, process-matched scale and aggregate the results afterwards rather than aggregating the inputs first. Modelling on fine support and then summarising predictions to the decision unit keeps the estimated relationships honest and makes the aggregation step explicit and auditable. Fitting the model directly on coarse units bakes the aggregation bias into the coefficients where it cannot be separated out.

How do I justify the choice to a reviewer?

State three numbers and one sentence: the process range you estimated and how, the data floor and what set it, and the decision unit if one exists — then the sentence that says which constraint bound the choice. Attach the sensitivity band from the sweep so the reader can see how much the result would move under the alternatives. A choice presented with its constraints is defensible even when someone would have chosen differently.

Related

← Back to Spatial Scale & the Modifiable Areal Unit Problem

Related

Measuring MAUP Sensitivity Across Aggregation Levels