Spatial Structure¶

Does a model's point-wise bias have exploitable spatial structure — does the error at one station tell you anything about the error at the next, and over what distance?

The empirical variogram answers this directly. Its two headline numbers are the ones that decide whether a bias field can be interpolated at all:

quantity	reads as
range	the distance beyond which errors are effectively unrelated
nugget	the irreducible point-scale component — variance that no amount of interpolation removes
relative nugget (`nugget / sill`)	how much of the field is noise. High means heavy smoothing and gauges that will not be reproduced exactly

Use the log-ratio, not the raw difference¶

Fit on log(model / obs) rather than model - obs. A multiplicative bias field is usually the stationary one; the raw difference typically keeps climbing with lag instead of reaching a clean sill, which means the stationarity assumption underpinning any interpolation does not hold.

fit_bias_variogram reports ch_n and mat_n — the Cressie–Hawkins and Matheron estimates near the origin — precisely so this is checkable rather than assumed. A ratio near 1 means the field is well behaved at short range; far from 1 means heavy tails that the classical estimator will exaggerate.

The nugget is fixed, not fitted¶

fit_exponential_variogram holds the nugget at a data-derived near-origin estimate and fits only (psill, range).

This is deliberate and worth understanding before changing it. On a typical station network there are too few sub-10 km pairs to anchor the intercept, and a free three-parameter fit collapses the nugget toward zero — which then reads as a smoothly interpolable bias field when the data say no such thing.

Quote a range with its resampling band, or not at all

With few stations the range is poorly identified — worst when the field is nugget-dominated, where many (range, sill) pairs fit a flattish cloud about equally well. bootstrap_variogram_params gives the p5/p50/p95 band. A range reported without it is a number with no error bar, and the band is often wide enough to change the conclusion.

That function takes a seed and builds its own generator rather than accepting one. This is load-bearing: callers typically invoke it once per field, and each call is meant to start from the same stream so results are independently reproducible. Threading a shared generator through instead silently changes every result after the first.

Example¶

import numpy as np
from pyproj import Transformer
from modverif.spatial import fit_bias_variogram, bootstrap_variogram_params

# Coordinates must be in a PROJECTED, metric CRS -- the lags are distances in km.
T = Transformer.from_crs(4326, 2193, always_xy=True)   # e.g. NZTM
gx, gy = T.transform(lons, lats)
x_km, y_km = np.asarray(gx) / 1000.0, np.asarray(gy) / 1000.0

fit = fit_bias_variogram(x_km, y_km, np.log(model / obs))
band = bootstrap_variogram_params(x_km, y_km, np.log(model / obs))

if fit['fit_ok']:
    print(f"range {fit['range_km']:.0f} km  (p5-p95 {band['range_km'][0]:.0f}-{band['range_km'][2]:.0f})")
    print(f"nugget effect {100 * fit['rel_nugget']:.0f}%")

API¶

Spatial-structure diagnostics for scattered model-minus-observation error.

The question this answers is whether a model's point-wise bias has exploitable spatial structure: does the error at one station tell you anything about the error at the next one, and over what distance? The empirical variogram is the standard tool, and its fitted range and nugget say directly how far a bias field can be interpolated and how much of it is irreducible point-scale noise.

Intended on the log-ratio log(model / obs) rather than the raw difference: a multiplicative bias field is usually the stationary one, and comparing the Cressie--Hawkins and Matheron estimators near the origin is a cheap test of which behaves better (see fit_bias_variogram's ch_n / mat_n).

The fitted length scale is called range_km throughout -- parameter and result key alike. Nothing here named rng is anything but a numpy.random.Generator.

`best_kmeans(xy, k, rng, n_restart=20)` ¶

k-means on point coordinates, keeping the lowest-inertia labelling over several restarts.

k-means converges to a local optimum that depends on initialisation, so a single run is a coin toss dressed as an answer. Restarts make the labelling reproducible in practice rather than only in principle.

NOTE: The restart seeds are drawn in one vectorised call. Replacing that with per-restart scalar draws consumes the caller's generator differently and changes every downstream result.

Parameters:

Name	Type	Description	Default
`xy`	`ndarray`	`(n, 2)` coordinates, in a projected metric CRS.	required
`k`	`int`	Number of clusters.	required
`rng`	`Generator`	Caller's generator, used only to seed the restarts.	required
`n_restart`	`int`	Restarts to attempt.	`20`

Returns:

Type	Description
`ndarray or None`	Cluster label per point, or None if every restart collapsed to fewer than `k` non-empty clusters -- which is the honest answer for a k the point set cannot support, and callers must handle it rather than assume an array.

`bootstrap_variogram_params(x_km, y_km, z, estimator='ch', n_boot=400, drop_frac=0.1, seed=0)` ¶

Resampling sensitivity of the fitted range and nugget.

Refits on n_boot random subsamples, each dropping drop_frac of the points, and reports p5/p50/p95 bands. With few points the range is poorly identified -- worst when the field is nugget-dominated, where many (range, sill) pairs fit a flattish cloud about equally well -- so its band is usually wide. The nugget, anchored to near-origin data, is much tighter. A range quoted without this band is a number with no error bar.

Subsampling without replacement, deliberately: resampling with replacement would create zero-distance duplicate pairs and corrupt the near-origin estimate.

WARNING: This function takes a seed and constructs its own generator; it does not accept one. That is intentional and load-bearing. Callers commonly invoke it more than once per run (say, for a log-ratio field and a raw-difference field), and each call is meant to start from the same stream so the two are independently reproducible. Threading one shared generator through instead -- the usual tidy-up -- silently changes every result after the first call.

Parameters:

Name	Type	Description	Default
`x_km`	`ndarray`	Point coordinates in a projected, metric CRS, in km.	required
`y_km`	`ndarray`	Point coordinates in a projected, metric CRS, in km.	required
`z`	`ndarray`	Field value per point.	required
`estimator`	`(ch, matheron)`	Passed through to `fit_bias_variogram`.	`'ch'`
`n_boot`	`int`	Number of subsample refits to attempt.	`400`
`drop_frac`	`float`	Fraction of points dropped per subsample; at least 6 points are always kept.	`0.1`
`seed`	`int`	Seed for this call's generator.	`0`

Returns:

Type	Description
`dict or None`	`range_km` and `nugget`, each a `(p5, p50, p95)` tuple, plus `drop_frac` and `n_ok` (how many refits converged). None if fewer than half the refits converged, since a band from a minority of fits would misrepresent the uncertainty rather than describe it.

`ch_gamma(abs_diffs)` ¶

Cressie--Hawkins robust semivariance from a set of |z_i - z_j| pair values.

Built on the mean square-root difference rather than the mean squared difference, so a single outlying pair cannot dominate the estimate. Prefer this to matheron_gamma on heavy-tailed fields; their ratio near the origin is itself a diagnostic.

Parameters:

Name	Type	Description	Default
`abs_diffs`	`ndarray`	Absolute pairwise differences of the field.	required

Returns:

Type	Description
`float`	Semivariance, or NaN if no pairs were supplied.

`empirical_variogram(distances, abs_diffs, estimator='ch', min_pairs_per_bin=MIN_PAIRS_PER_BIN, max_lag_pct=MAX_LAG_PCT)` ¶

Bin scattered pairs by separation distance and estimate semivariance in each bin.

The bin count is derived from the number of pairs rather than fixed, so the same call works for a few dozen stations and for a few thousand without retuning.

Parameters:

Name	Type	Description	Default
`distances`	`ndarray`	Condensed pairwise distances (as from `scipy.spatial.distance.pdist`), in km.	required
`abs_diffs`	`ndarray`	Absolute pairwise field differences, in the same condensed order as `distances`.	required
`estimator`	`(ch, matheron)`	`'ch'` selects `ch_gamma`, anything else `matheron_gamma`.	`'ch'`
`min_pairs_per_bin`	`int`	Bins with fewer pairs are dropped entirely rather than reported noisily.	`MIN_PAIRS_PER_BIN`
`max_lag_pct`	`float`	Percentile of `distances` at which to cap the lag axis.	`MAX_LAG_PCT`

Returns:

Name	Type	Description
`centers`	`ndarray`	Bin centre distances, km.
`gammas`	`ndarray`	Semivariance per retained bin.
`counts`	`ndarray`	Pair count per retained bin -- the natural fit weight.
`max_lag`	`float`	The lag cap actually used, km.

`exponential_variogram(h, nugget, psill, range_km)` ¶

Isotropic exponential semivariogram model.

Parameters:

Name	Type	Description	Default
`h`	`ndarray`	Lag distances, km.	required
`nugget`	`float`	Semivariance at zero separation -- the irreducible point-scale component.	required
`psill`	`float`	Partial sill; the model asymptotes to `nugget + psill`.	required
`range_km`	`float`	e-folding length scale, km. Reported under the same name in the fit result dicts.	required

Returns:

Type	Description
`ndarray`	Modelled semivariance at each lag.

`fit_bias_variogram(x_km, y_km, z, estimator='ch', near_origin_km=NEAR_ORIGIN_KM)` ¶

Full variogram pipeline for a scattered field: empirical estimate plus a nugget-fixed fit.

Single entry point so a diagnostic and any downstream interpolation cannot disagree about which variogram they are using.

Parameters:

Name	Type	Description	Default
`x_km`	`ndarray`	Point coordinates in a projected, metric CRS, in km.	required
`y_km`	`ndarray`	Point coordinates in a projected, metric CRS, in km.	required
`z`	`ndarray`	Field value per point -- typically `log(model / obs)`.	required
`estimator`	`(ch, matheron)`	Semivariance estimator, and which near-origin estimate becomes the fixed nugget.	`'ch'`
`near_origin_km`	`float`	Separation below which pairs count toward the nugget estimate.	`NEAR_ORIGIN_KM`

Returns:

Type Description

dict

The fit parameters (nugget, psill, sill, range_km, rel_nugget), the empirical points (centers, gammas, counts, max_lag), and diagnostics (var_z, ch_n, mat_n, n_near, fit_ok). Fit parameters are NaN when fit_ok is False.

ch_n / mat_n is a stationarity read: a ratio near 1 means the field is well behaved at short range, and a ratio far from 1 means heavy tails that a Matheron estimate will exaggerate.

`fit_exponential_variogram(centers, gammas, counts, var_z, max_lag, d_min, nugget)` ¶

Fit an isotropic exponential variogram with the nugget fixed, not free.

Fixing it is deliberate. With a typical station network there are too few sub-10 km pairs to anchor the intercept, and a free three-parameter fit collapses the nugget toward zero -- which then reads as a smoothly interpolable bias field when the data say no such thing. Only (psill, range) are fitted, weighted by pair count.

Parameters:

Name	Type	Description	Default
`centers`	`ndarray`	Empirical variogram points from `empirical_variogram`.	required
`gammas`	`ndarray`	Empirical variogram points from `empirical_variogram`.	required
`counts`	`ndarray`	Empirical variogram points from `empirical_variogram`.	required
`var_z`	`float`	Variance of the field; bounds the sill.	required
`max_lag`	`float`	Upper bound for the fitted range, km.	required
`d_min`	`float`	Smallest pairwise distance; lower bound for the fitted range.	required
`nugget`	`float`	The data-derived near-origin estimate to hold fixed. Clipped into `[0, var_z]`.	required

Returns:

Type	Description
`dict or None`	`nugget`, `psill`, `sill`, `range_km`, `rel_nugget`; or None if the fit did not converge or there were too few points to attempt one.

`matheron_gamma(abs_diffs)` ¶

Classical Matheron semivariance from a set of |z_i - z_j| pair values.

Parameters:

Name	Type	Description	Default
`abs_diffs`	`ndarray`	Absolute pairwise differences of the field.	required

Returns:

Type	Description
`float`	Semivariance, or NaN if no pairs were supplied.

`morans_i(zc, weights, s0, n, denom)` ¶

Moran's I for a pre-centred field under a given weight matrix.

Takes s0, n and denom as arguments rather than deriving them because callers evaluate many weight matrices, or many permutations, against the same field -- recomputing the invariants each time is the dominant cost.

Parameters:

Name	Type	Description	Default
`zc`	`ndarray`	Field with its mean already removed.	required
`weights`	`ndarray`	`(n, n)` spatial weight matrix, zero on the diagonal.	required
`s0`	`float`	Sum of `weights`.	required
`n`	`int`	Number of locations.	required
`denom`	`float`	`zc @ zc`, the field's total squared deviation.	required

Returns:

Type	Description
`float`	Moran's I. Its null expectation is `-1 / (n - 1)`, not zero -- a small positive value can still be below chance.

`morans_i_at_points(values, x, y, rng, n_perm=999, min_points=8)` ¶

Moran's I and its permutation p-value for scattered points, with inverse-distance weights.

The convenience layer over morans_i and morans_i_permutation, which take pre-computed weight matrices and invariants because callers evaluating many bands or many permutations should not recompute them. That is the right shape for a hot loop and the wrong shape for the common case: "is this scattered field spatially clustered?"

This exists because the low-level pair was, in practice, easier to re-implement than to call -- which is a design smell worth fixing rather than documenting.

Parameters:

Name	Type	Description	Default
`values`	`ndarray`	Field value per point.	required
`x`	`ndarray`	Point coordinates in a projected metric CRS. Units set the weight scale; only relative distances matter.	required
`y`	`ndarray`	Point coordinates in a projected metric CRS. Units set the weight scale; only relative distances matter.	required
`rng`	`Generator`	Caller's generator, consumed once per permutation.	required
`n_perm`	`int`	Permutations for the null.	`999`
`min_points`	`int`	Below this, `(NaN, NaN)` is returned rather than a number -- Moran's I on a handful of points is dominated by the weight matrix's geometry, not the field.	`8`

Returns:

Name	Type	Description
`morans_i`	`float`	Observed statistic. Its null expectation is `-1 / (n - 1)`, not zero.
`p_value`	`float`	One-sided permutation p-value for positive autocorrelation.

`morans_i_permutation(zc, weight_mats, n, denom, rng, n_perm=999)` ¶

Permutation null for Moran's I across one or more weight matrices.

WARNING: One shuffle is evaluated against every weight matrix, by design. Testing each matrix with its own independent permutations would draw n_perm x len(weight_mats) times instead of n_perm, and -- more importantly -- would destroy the correlation between the bands' nulls. A correlogram's bands are not independent tests of independent quantities; they are one field viewed at several scales, and a multiple-comparison correction applied across independently-generated nulls is answering a different question from the one asked.

This is why the function takes a list of weight matrices rather than being called once per matrix. Refactoring it into a per-matrix helper changes both the draw count and the result.

Parameters:

Name	Type	Description	Default
`zc`	`ndarray`	Field with its mean already removed.	required
`weight_mats`	`sequence of (np.ndarray, float)`	`(weights, s0)` pairs -- e.g. one per distance band of a correlogram, or a single entry for a global statistic.	required
`n`	`int`	Number of locations.	required
`denom`	`float`	`zc @ zc`.	required
`rng`	`Generator`	Caller's generator. Consumed once per permutation, in order.	required
`n_perm`	`int`	Number of permutations.	`999`

Returns:

Name	Type	Description
`observed`	`ndarray`	Moran's I per weight matrix on the real data.
`null`	`ndarray`	`(n_perm, len(weight_mats))` null statistics.

Spatial Structure¶

Use the log-ratio, not the raw difference¶

The nugget is fixed, not fitted¶

Example¶

API¶

best_kmeans(xy, k, rng, n_restart=20) ¶

bootstrap_variogram_params(x_km, y_km, z, estimator='ch', n_boot=400, drop_frac=0.1, seed=0) ¶

ch_gamma(abs_diffs) ¶

empirical_variogram(distances, abs_diffs, estimator='ch', min_pairs_per_bin=MIN_PAIRS_PER_BIN, max_lag_pct=MAX_LAG_PCT) ¶

exponential_variogram(h, nugget, psill, range_km) ¶

fit_bias_variogram(x_km, y_km, z, estimator='ch', near_origin_km=NEAR_ORIGIN_KM) ¶

fit_exponential_variogram(centers, gammas, counts, var_z, max_lag, d_min, nugget) ¶

matheron_gamma(abs_diffs) ¶

morans_i(zc, weights, s0, n, denom) ¶

morans_i_at_points(values, x, y, rng, n_perm=999, min_points=8) ¶

morans_i_permutation(zc, weight_mats, n, denom, rng, n_perm=999) ¶

`best_kmeans(xy, k, rng, n_restart=20)` ¶

`bootstrap_variogram_params(x_km, y_km, z, estimator='ch', n_boot=400, drop_frac=0.1, seed=0)` ¶

`ch_gamma(abs_diffs)` ¶

`empirical_variogram(distances, abs_diffs, estimator='ch', min_pairs_per_bin=MIN_PAIRS_PER_BIN, max_lag_pct=MAX_LAG_PCT)` ¶

`exponential_variogram(h, nugget, psill, range_km)` ¶

`fit_bias_variogram(x_km, y_km, z, estimator='ch', near_origin_km=NEAR_ORIGIN_KM)` ¶

`fit_exponential_variogram(centers, gammas, counts, var_z, max_lag, d_min, nugget)` ¶

`matheron_gamma(abs_diffs)` ¶

`morans_i(zc, weights, s0, n, denom)` ¶

`morans_i_at_points(values, x, y, rng, n_perm=999, min_points=8)` ¶

`morans_i_permutation(zc, weight_mats, n, denom, rng, n_perm=999)` ¶