krigekit#

Python wrapper for a Fortran kriging and SGSIM engine.

The shared library (libkriging.so / kriging.dll) must be compiled from the Fortran sources in the fortran/ directory and placed in this package directory before use. See the README for build instructions.

Public API#

Spatial usage:

k = Kriging(ndim=2, nvar=1)

Space-time usage:

k = SpaceTimeKriging(nvar=1)
k.set_st_model(...)

Classes (returned by Kriging factory): Kriging — spatial kriging / co-kriging / SGSIM SpaceTimeKriging — 3-D + time kriging / SGSIM
Convenience functions: ordinary_kriging — one-shot point kriging cokriging — one-shot co-kriging sequential_gaussian_simulation — one-shot SGSIM spacetime_kriging — one-shot ST kriging spacetime_cokriging — one-shot ST co-kriging

Submodules#

Classes#

`Kriging`	Python interface to the Fortran t_kriging spatial kriging/simulation engine.
`SpaceTimeKriging`	Python interface to the Fortran t_kriging_st space-time kriging engine.
`IndicatorKriging`	Multiple Indicator Kriging and Sequential Indicator Simulation,

Functions#

`ordinary_kriging`(→ tuple[numpy.ndarray, numpy.ndarray])	One-shot ordinary kriging with a single isotropic (or anisotropic) variogram.
`cokriging`(→ tuple[numpy.ndarray, numpy.ndarray])	One-shot ordinary co-kriging with multiple variables.
`sequential_gaussian_simulation`(→ numpy.ndarray)	Sequential Gaussian Simulation.
`spacetime_kriging`(→ tuple[np.ndarray, np.ndarray])	One-shot ordinary space-time kriging (single variable).
`spacetime_cokriging`(→ tuple[np.ndarray, np.ndarray])	One-shot ordinary space-time co-kriging.

Package Contents#

class krigekit.Kriging(ndim: int = 2, nvar: int = 1, ndrift: int = 0, unbias: int = 1, nsim: int = 0, anisotropic_search: bool = False, weight_correction: bool = False, use_old_weight: bool = False, store_weight: bool = False, cross_validation: bool = False, write_mat: bool = False, neglect_error: bool = True, varying_vgm: bool = False, std_ck: bool = False, verbose: bool = False, pf_cache: bool = False, weight_file: str = '', bounds: tuple | None = None, seed: int | None = None)#

Python interface to the Fortran t_kriging spatial kriging/simulation engine.

Array convention#

All coordinate arrays use (nobs, ndim) shape — rows are points, columns are spatial dimensions. This matches NumPy, pandas, and scikit-learn conventions. The wrapper transparently transposes to Fortran’s (ndim, nobs) before calling the library.

Typical workflow#

>>> k = Kriging(ndim=2, nvar=1)
>>> k.set_obs(ivar=1, coord=obs_coord, value=obs_value, nmax=20)
>>> k.set_grid(coord=grid_coord)
>>> k.set_vgm(ivar=1, jvar=1, vtype="sph", nugget=0, sill=1.0, a_major=1000, a_minor1=500, a_minor2=50)
>>> k.set_search(ivar=1)
>>> k.solve()
>>> estimate, variance = k.get_results()
>>> del k    # release memory

For sequential Gaussian simulation add nsim=N to the constructor and call set_sim() after set_grid().

set_obs(ivar: int, coord: numpy.ndarray, value: numpy.ndarray, variance: numpy.ndarray | None = None, nmax: int | None = None, maxdist: float | None = None, sk_mean: float = 0.0)#

Set observations for variable ivar.

Drift values are set separately via set_obs_drift() after this call, when ndrift > 0.

Parameters:

ivar (int) – Variable index, 1-based.
coord (ndarray, shape (nobs, ndim)) – Duplicate coordinate tuples within the same variable are rejected. Observation coordinates. Rows are points, columns are spatial dimensions — standard Python/NumPy convention. The wrapper transposes to Fortran’s (ndim, nobs) internally.
value (ndarray, shape (nobs,)) – Observed values.
variance (ndarray, shape (nobs,), optional) – Per-observation measurement error variance added to the diagonal of the covariance matrix. Defaults to zeros (no measurement error).
nmax (int, optional) – Maximum number of neighbours. Default: use all observations.
maxdist (float, optional) – Maximum search distance. Default: unlimited.
sk_mean (float) – Global mean for simple kriging (unbias=0). Default 0.0.

set_obs_drift(ivar: int, drift: numpy.ndarray)#

Set external drift values at observation locations for variable ivar.

Call after set_obs() for the same ivar. Only needed when ndrift > 0 was passed to the constructor.

Parameters:

ivar (int) – Variable index, 1-based.
drift (ndarray, shape (nobs, ndrift)) – Drift values. Rows are observations, columns are drift functions. Transposed to (ndrift, nobs) internally before calling Fortran.

update_obs_value(ivar: int, value: numpy.ndarray)#

Replace observation values for variable ivar in-place.

Coordinates and the kd-tree are unchanged. The primary use case is weight reuse: after solving once (with store_weight=True or use_old_weight=True), call this method with new observed values at the same locations and call solve() again to get updated estimates without recomputing search neighbourhoods or the LHS factorization.

Parameters:

ivar (int) – Variable index, 1-based.
value (ndarray, shape (nobs,)) – New observed values. Length must match the nobs passed to the previous set_obs() call for this variable.

set_vgm(ivar: int, jvar: int, vtype: str, nugget: float = 0.0, sill: float = 1.0, a_major: float = 1.0, a_minor1: float | None = None, a_minor2: float | None = None, azimuth: float = 0.0, dip: float = 0.0, plunge: float = 0.0, append: bool = True, product: bool = False)#

Add one nested variogram structure for the (ivar, jvar) pair. Call multiple times with append=True to build a nested (multi-structure) model. Pass append=False to clear any previously set structures for the pair before adding this one (useful when reusing a Kriging object with a different variogram).

Parameters:

ivar (int) – Variable indices (1-based). Use ivar=jvar for auto-variograms, ivar≠jvar for cross-variograms. The LMC constraint b12² ≤ b11 × b22 must be satisfied for each nested structure.
jvar (int) – Variable indices (1-based). Use ivar=jvar for auto-variograms, ivar≠jvar for cross-variograms. The LMC constraint b12² ≤ b11 × b22 must be satisfied for each nested structure.
vtype (str) – Variogram type: one of sph, exp, gau, pow, lin, hol, bsq, cir, nug, cyc, dco. cyc is the GP periodic (exp-sine-squared) kernel where a_major is the period; dco is a damped cosine.
nugget (float) – Nugget contribution of this structure (default 0).
sill (float) – Partial sill of this structure (default 1).
a_major (float) – Range along the major axis (default 1).
a_minor1 (float, optional) – Range along the first minor axis. Defaults to a_major (isotropic in the horizontal plane).
a_minor2 (float, optional) – Range along the second minor axis. Defaults to a_minor1.
azimuth (float) – Rotation angles in degrees (default 0).
dip (float) – Rotation angles in degrees (default 0).
plunge (float) – Rotation angles in degrees (default 0).
append (bool) – True (default) appends this structure to any existing ones. False clears the current model before adding.
product (bool) –
False (default) adds this structure to the previous ones (standard additive nesting). True multiplies this structure with the immediately preceding one, forming a product group. The Schur product of two positive-definite covariances is also positive-definite, so products are always valid.

Example — damped cosine with independent decay and period:
```
k.set_vgm(1, 1, vtype="exp", sill=1.0, a_major=5.0)   # decay over 5 yr
k.set_vgm(1, 1, vtype="hol", sill=1.0, a_major=0.5,   # 1-yr period
          product=True)
```

Example

>>> k.set_vgm(1, 1, vtype="sph", nugget=0.0, sill=1.0, a_major=500.0)
>>> k.set_vgm(1, 1, vtype="nug", nugget=0.1, sill=0.0, a_major=1.0)
>>> k.set_vgm(1, 1, vtype="sph", nugget=0.0, sill=0.9, a_major=500.0)

set_vgm_block(ib: int, ivar: int, jvar: int, vtype: str, nugget: float = 0.0, sill: float = 1.0, a_major: float = 1.0, a_minor1: float | None = None, a_minor2: float | None = None, azimuth: float = 0.0, dip: float = 0.0, plunge: float = 0.0)#

Add one nested variogram structure for a specific block ib.

Requires varying_vgm=True in the constructor and set_grid() to have been called first (because the number of blocks must be known before the per-block variogram array can be allocated in Fortran).

Call multiple times for the same ib to build a nested model.

Parameters:

ib (int) – Block index (1-based).
ivar (int) – Variable indices (1-based).
jvar (int) – Variable indices (1-based).
vtype (str) – Variogram type: sph, exp, gau, pow, lin, hol, bsq, cir, nug, cyc, or dco.
nugget (float) – Nugget contribution (default 0).
sill (float) – Partial sill (default 1).
a_major (float) – Range along the major axis (default 1).
a_minor1 (float, optional) – First minor-axis range (defaults to a_major).
a_minor2 (float, optional) – Second minor-axis range (defaults to a_minor1).
azimuth (float) – Rotation angles in degrees (default 0).
dip (float) – Rotation angles in degrees (default 0).
plunge (float) – Rotation angles in degrees (default 0).

set_grid(coord: numpy.ndarray | None = None, rangescale: numpy.ndarray | None = None, localnugget: numpy.ndarray | None = None)#

Set the estimation grid for point kriging (one node per block).

For block kriging use set_grid_block(). For cross-validation use set_grid_cv(). Drift is set separately via set_grid_drift() when ndrift > 0.

Parameters:

coord (ndarray, shape (ngrid, ndim)) – Grid coordinates. Rows are grid nodes, columns are spatial dimensions.
rangescale (ndarray, shape (ngrid,), optional) – Per-block variogram range scaling factor. Values > 1 increase the effective range, useful to account for data sparsity. Default: 1.0 for all blocks.
localnugget (ndarray, shape (ngrid,), optional) – Additional nugget added per block to model local uncertainty. Default: 0.0 for all blocks.

set_grid_block(coord: numpy.ndarray, block_type: int, nblockpnt: numpy.ndarray, pointweight: numpy.ndarray | None = None, blocksize: numpy.ndarray | None = None, rangescale: numpy.ndarray | None = None, localnugget: numpy.ndarray | None = None)#

Set the estimation grid for block kriging.

Drift is set separately via set_grid_drift() when ndrift > 0.

Parameters:

coord (ndarray, shape (ngrid, ndim)) – Sub-node coordinates across all blocks (total ngrid = sum(nblockpnt)).
block_type (int) – -4 = Gaussian quadrature nodes (auto-generated); >0 = user-supplied sub-nodes (coord contains sub-node positions).
nblockpnt (ndarray of int, shape (nblock,)) – Number of sub-nodes per block.
pointweight (ndarray, shape (sum(nblockpnt),), optional) – Weight of each sub-node. Uniform weights (1/nblockpnt) used if omitted.
blocksize (ndarray, shape (nblock,ndim), optional) – Block size in each dimension when block_type == -4.
rangescale (ndarray, shape (nblock,), optional) – Per-block variogram range scaling. Default: 1.0.
localnugget (ndarray, shape (nblock,), optional) – Per-block additional nugget. Default: 0.0.

set_grid_cv()#

Set up the grid for cross-validation mode.

No coordinate argument is needed — Fortran derives the grid from the observation coordinates automatically. Call instead of set_grid() when cross_validation=True was passed to the constructor.

set_grid_drift(drift: numpy.ndarray, ivar: int | None = None)#

Set external drift values at grid/block locations.

Call after set_grid(), set_grid_block(), or set_grid_cv(). Only needed when ndrift > 0.

Parameters:

drift (ndarray, shape (nblocks, ndrift)) – Drift values. Rows are blocks, columns are drift functions. Note: use nblocks (number of blocks), not ngrid (number of sub-nodes), even for block kriging. Transposed to (ndrift, nblocks) internally before calling Fortran.
ivar (int, optional) – Target-variable index (1-based) whose RHS receives this drift. None (default) broadcasts the same drift to all target variables — the usual case when external drift is independent of which variable is being estimated.

Note

ivar here refers to the target variable (which variable’s estimate uses this drift in its RHS), not the source variable. This is the opposite end from set_obs_drift(), whose ivar identifies the source variable (whose observations form the F-matrix column).

set_sim(randpath: numpy.ndarray | None = None, sample: numpy.ndarray | None = None)#

Set up Sequential Gaussian Simulation parameters.

Call after set_grid() and before set_search(). Only needed when nsim > 0.

Parameters:

randpath (ndarray of int, shape (nblocks,), optional) – 1-based random visiting order for the block loop. When None, Fortran generates a random permutation.
sample (ndarray, shape (nblocks, nvar, nsim), optional) – Pre-drawn samples. When None, Fortran generates them via the object’s set_sim override: N(0,1) for Kriging, U(0,1) for IndicatorKriging.

set_search(ivar: int = 1, anis1: float = 1.0, anis2: float = 1.0, azimuth: float = 0.0, dip: float = 0.0, plunge: float = 0.0, sector_search: bool = False)#

Build the KD-tree and configure the search ellipse for variable ivar. Call once per variable after set_obs() (and set_sim() for SGSIM).

Parameters:

ivar (int) – Variable index (1-based).
anis1 (float) – Horizontal anisotropy ratio (minor / major range). 1.0 = isotropic.
anis2 (float) – Vertical anisotropy ratio (vertical / major range). 1.0 = isotropic.
azimuth (float) – Azimuth of the major axis (degrees, clockwise from North).
dip (float) – Dip angle (degrees, positive downward).
plunge (float) – Plunge angle (degrees).
sector_search (bool) – Enable sector (quadrant in 2D, octant in 3D) search limiting candidates per sector. If True, the search space is divided into quadrants/octants centered on the prediction location. Candidates are sorted by distance, and at most nmax (from set_obs()) are selected from each sector. This ensures a balanced spatial distribution of neighbours and prevents clustering artifacts. The maximum total neighbours selected is 4 * nmax in 2D or 8 * nmax in 3D. If search anisotropy is enabled, coordinates are rotated and scaled according to the anisotropy parameters before sector assignment.

set_grad(coord1: numpy.ndarray, coord2: numpy.ndarray, grad_value: numpy.ndarray, ivar: int = 1, variance: numpy.ndarray | None = None, drift_ext: numpy.ndarray | None = None)#

Set gradient observation pairs (Delhomme 1979: “Kriging in hydrology”).

Each pair (coord1[i], coord2[i]) approximates the directional gradient at a boundary as a finite difference. The constraint Z(xs1) - Z(xs2) = grad_value[i] is enforced as a hard equality. For a no-flow (zero normal gradient) boundary use grad_value = 0.

Call set_grad() after set_search() and before solve().

Parameters:

coord1 (ndarray, shape (ngrad, ndim)) – Positive-side virtual node coordinates.
coord2 (ndarray, shape (ngrad, ndim)) – Negative-side virtual node coordinates.
grad_value (ndarray, shape (ngrad,)) – Known gradient values. Use 0 for no-flow boundaries.
ivar (int, default 1) – Variable index (1-based) the gradient pairs constrain. For cokriging, specifies which variable’s gradient is observed.
variance (ndarray, shape (ngrad,), optional) – Gradient observation variance (default 0 = exact constraint). A non-zero value relaxes the constraint, analogous to obs nugget.
drift_ext (ndarray, shape (ngrad, ndrift), optional) – External drift differences f_ext(xs1) - f_ext(xs2) for each pair. Required when ndrift > 0; omit for ordinary kriging.

solve(nthread: int = 0, ncache: int | None = None)#

Run the kriging or SGSIM loop over all blocks. Calls prepare(), then the parallel block loop internally.

Parameters:

nthread (int, optional) – Number of OpenMP threads to use for this call. 0 (default) leaves the OMP runtime setting unchanged. 1 forces single-threaded execution (useful for reproducible results or when calling solve() from inside another parallel region).
ncache (int, optional) – Number of per-thread multi-slot hcache entries to use for this solve call. None keeps the compiled/object default, 0 disables the hcache, and 1 gives a one-slot hcache for cache-overhead comparisons. The single-entry ctx%cache and optional persistent factor cache are unaffected.

property solver_stats: dict#

Solver statistics from the most recent solve() call.

Returns a dict with five integer counts, reset to zero at the start of every solve():

fail: Blocks where both Cholesky and SSYTRF failed; solution set to NaN (only possible when neglect_error is True).
chol_fact: Fresh Cholesky factorizations performed — O(n³) each, one per unique neighbourhood on a cache miss.
chol_reuse: Blocks solved via a cached Cholesky factorization — O(n²) each. A large value relative to chol_fact means the neighbourhood cache is working effectively.
ssytrf_fact: SSYTRF (Bunch-Kaufman LDL^T) factorizations performed — O(n³) each, once per unique neighbourhood. Non-zero means Cholesky failed for at least one neighbourhood (e.g. a non-SPD system).
ssytrf_reuse: Blocks solved by a cached SSYTRF via SSYTRS — O(n²) each.

get_results(copy: bool = False, squeeze: bool = True)#

Retrieve the kriging estimates and variances after solve().

Fortran fills estimate(nsim, nblocks) directly into a Fortran-contiguous Python-owned buffer.

Parameters:

copy (bool, default False) – If True, return C-contiguous copies for downstream NumPy/Pandas use. If False, return views / Fortran-order arrays when possible.
squeeze (bool, default True) – If True, return a 1-D estimate when nsim == 1.

Returns:

estimate (ndarray) – (ngrid, nvar, nsim). Shape (ngrid,) when (nsim == 1 or ==1) and squeeze; otherwise shape (ngrid, nvar, nsim).
variance (ndarray, shape (nblocks,)) – Kriging variance at each block.

Example

>>> est, var = k.get_results()
>>> kriging_estimate = est[0]          # shape (nblocks,)
>>> sim_realisation1 = est[0]          # same for nsim=1

get_result_array() → numpy.ndarray#

Return all results as a NumPy structured (record) array — one row per block.

The array contains block centroid coordinates alongside every estimate, simulation realization, and variance produced by the last solve().

Fields#

Coordinates (always present):: x, y [, z] — block centroid coordinates.
Estimates / simulations:: Kriging (nsim == 0), nvar == 1: estimate Kriging (nsim == 0), nvar > 1: est_v1, est_v2, … SGSIM (nsim > 0), nvar == 1: sim_1, sim_2, …, sim_{nsim} SGSIM (nsim > 0), nvar > 1: v1_s1, v1_s2, …, v{nvar}_s{nsim}
Variances (always present, diagonal of conditional covariance):: nvar == 1: variance nvar > 1: var_v1, var_v2, …

returns:: np.ndarray with named fields (structured array / record array) – Shape (nblocks,). Access a column with arr['estimate'], convert to a plain 2-D array with np.column_stack([arr[f] for f in arr.dtype.names]).

Example

>>> k.solve()
>>> ra = k.get_result_array()
>>> ra.dtype.names
('x', 'y', 'estimate', 'variance')
>>> ra['estimate']          # 1-D array, shape (nblocks,)
>>> import pandas as pd
>>> df = pd.DataFrame(ra)   # direct conversion to DataFrame

get_result_df() → pandas.DataFrame#

Return all results as a pandas.DataFrame.

Wraps get_result_array(); column names and field descriptions are identical to those documented there.

Returns:: pandas.DataFrame – One row per block; columns match the fields of the structured array returned by get_result_array().

Example

>>> k.solve()
>>> df = k.get_result_df()
>>> df.columns.tolist()
['x', 'y', 'estimate', 'variance']

get_factor() → dict#

Return the persistent LHS factorization cached after the last solve().

The Cholesky factorization of the kriging covariance matrix K and the related Schur-complement matrices are computed once per solve() call (or reused across blocks with the same neighbour set). Starting from the second solve() call on unchanged observations and variogram, the cached factors allow the Fortran engine to skip kriging_setup entirely.

This method exposes those matrices and the assembled linear system so that users can inspect or verify the factorization inputs.

Returns:

dict with keys
``valid`` (bool) – True if a persistent factor exists (i.e. at least one solve() has been completed and observations/variogram have not changed since).
``npp`` (int) – Number of neighbours in the LHS matrix (size of K).
``p`` (int) – Number of drift + unbiasedness columns (size of the Schur complement).
``L`` (ndarray, shape (npp, npp)) – Lower-triangular Cholesky factor of K (stored column-major by Fortran, returned as a C-contiguous array).
``kinv_drift`` (ndarray, shape (npp, max(1, p))) – K⁻¹ F (K inverse applied to the drift matrix F).
``schur`` (ndarray, shape (max(1, p), max(1, p))) – Cholesky factor of the Schur complement F’ K⁻¹ F.
``matA`` (ndarray, shape (npp + p, npp + p)) – Assembled linear-system LHS before factorization.
``rhsB`` (ndarray, shape (nvar, npp + p)) – Assembled linear-system RHS before solving.

Example

>>> k = Kriging(ndim=2, nvar=1, ndrift=1, unbias=0)
>>> # ... set_obs, set_vgm, set_grid, set_obs_drift, set_grid_drift ...
>>> k.solve()
>>> f = k.get_factor()
>>> if f['valid']:
...     L = f['L']            # Cholesky factor of covariance matrix
...     kinv = f['kinv_drift']  # K^{-1} F
...     schur = f['schur']    # Cholesky of Schur complement

Notes

The factor is invalidated (valid = False) whenever set_obs() or set_vgm() is called. Call solve() again to repopulate.

get_estimate_all(copy: bool = False)#

Return multivariable estimates / simulations for all variables.

Populated when nvar > 1. For co-kriging without simulation, the leading dimension is 1.

Parameters:: copy (bool, default False) – If True, return a C-contiguous copy. If False, return the Fortran-contiguous output buffer filled by the Fortran core.
Returns:: np.ndarray, shape (nblock, nvar, max(nsim, 1)) – Values of all variables. out[ib, kvar, isim] is the value at block ib+1 for variable kvar+1 in realization isim+1.

get_variance_all(copy: bool = False)#

Return the conditional covariance matrix for all variables.

Returns:: np.ndarray, shape (nblock, nvar, nvar) – Conditional covariance matrix at each block. The diagonal contains each variable’s kriging variance, and out[:, 0, 0] matches the variance returned by get_results().

free_weight_store()#: Release the in-memory weight store, freeing its memory.

set_weights(weights: dict) → None#

Load kriging weights into the in-memory store so solve() reuses them.

Activated when use_old_weight=True and no weight_file is given. solve() then applies the supplied neighbour indices and kriging weights directly — skipping the kriging-system solve — and restores the stored variance. This is the in-memory equivalent of the use_old_weight=True + factor-file workflow.

Typical workflow#

>>> # First run: solve and capture weights + variance
>>> k1 = Kriging(ndim=2, nvar=1, store_weight=True)
>>> k1.set_obs(...); k1.set_grid(...); k1.set_vgm(...); k1.set_search(...)
>>> k1.solve()
>>> w = k1.get_weights()           # {'nnear', 'inear', 'weight', 'variance'}
>>>
>>> # Second run: same grid/vgm, new obs values, reuse weights
>>> k2 = Kriging(ndim=2, nvar=1, use_old_weight=True)   # no weight_file
>>> k2.set_obs(...new_values...)
>>> k2.set_grid(...); k2.set_vgm(...); k2.set_search(...)
>>> k2.set_weights(w)              # populate the in-memory store
>>> k2.solve()                     # fast: skips kriging system solve
>>> est, var = k2.get_results()

param weights:

Dict as returned by get_weights(), with keys:

nnearndarray (nblock, ngroups), int32: Number of active neighbours per block and group.
inearndarray (nblock, ngroups, nmax), int32: 1-based neighbour indices.
weightndarray (nblock, nvar, ngroups, nmax), float64: Kriging weights. For nvar==1 the array may be 3-D (nblock, ngroups, nmax) (the shape returned by get_weights()).
orderndarray (nblock,), optional: Random visiting order for SGSIM. Defaults to None.
variancendarray (nblock,) or (nblock, nvar, nvar), optional: Per-block conditional variance. Included automatically when the dict was produced by get_weights(). Defaults to zeros if absent.

type weights:

dict

Notes

Call after set_obs(), set_grid(), set_vgm(), and set_search(). The Fortran-side use_old_weight flag is set automatically; you may also pass use_old_weight=True to the constructor to declare intent explicitly.

get_weights() → dict#

Return the stored kriging weights and neighbour indices.

alloc_weight_store() must have been called before solve().

Returns:

dict with keys
``nnear`` (ndarray, shape (nblock, ngroups), dtype int32) – Number of active neighbours for each block and group. ngroups = ngroups_base when set_grad has not been called, or ngroups_base + nvar when gradient data is present. ngroups_base = nvar (kriging) or 2*nvar (SGSIM). Group layout:

indices 0..nvar-1 — obs groups (variable 1..nvar) indices nvar..2*nvar-1 — sim groups (SGSIM only) indices ngroups_base..ngroups-1 — grad groups (present only when set_grad called)
``inear`` (ndarray, shape (nblock, ngroups, nmax), dtype int32) – 1-based neighbour indices. Entries beyond nnear[ib, ig] are zero.
``weight`` (ndarray, shape (nblock, nvar, ngroups, nmax), dtype float64) – Kriging weights. Entries beyond nnear[ib, ig] are zero. Shape is (nblock, ngroups, nmax) when nvar == 1.
``variance`` (ndarray, shape (nblock,) for nvar==1, else (nblock, nvar, nvar)) – Per-block conditional kriging variance stored alongside the weights. Present only when the compiled library supports the variance store (i.e. built with the current source). Pass this dict directly to set_weights() to get a full round-trip.

get_info()#

krigekit.ordinary_kriging(obs_coord: numpy.ndarray, obs_value: numpy.ndarray, grid_coord: numpy.ndarray, vgm_spec: dict | list[dict], nmax: int | None = None, maxdist: float | None = None, search_anis1: float = 1.0, search_anis2: float = 1.0, search_azimuth: float = 0.0, rangescale: float | None = None, localnugget: float | None = None, nthread=0, ncache: int | None = None) → tuple[numpy.ndarray, numpy.ndarray]#

One-shot ordinary kriging with a single isotropic (or anisotropic) variogram.

Parameters:

obs_coord (ndarray, shape (nobs, ndim)) – Observation coordinates. Rows are points, columns are spatial dimensions.
obs_value (ndarray, shape (nobs,)) – Observation values.
grid_coord (ndarray, shape (ngrid, ndim)) – Grid coordinates to estimate.
vgm_spec (dict or list of dict) – One variogram structure dict, or a list of dicts for nested models. Each dict is passed as keyword arguments to Kriging.set_vgm() (keys: vtype, nugget, sill, a_major, and optionally a_minor1, a_minor2, azimuth, dip, plunge).
nmax (int) – Maximum number of neighbours.
maxdist (float, optional) – Maximum search distance.
search_anis1 (float) – Anisotropy ratios for search ellipse (1.0 = isotropic).
search_anis2 (float) – Anisotropy ratios for search ellipse (1.0 = isotropic).
search_azimuth (float) – Azimuth of search ellipse major axis (degrees from North).
nthread (int) – max OMP threads for this call (0 or absent = OMP default)
ncache (int, optional) – Per-thread hcache slots for this solve. None uses the default.

Returns:

estimate (ndarray, shape (ngrid,))
variance (ndarray, shape (ngrid,))

Example

>>> est, var = ordinary_kriging(
...     obs_coord, obs_value, grid_coord,
...     vgm_spec=dict(vtype="sph", nugget=100, sill=900, a_major=1000, a_minor1=500),
...     nmax=20)

krigekit.cokriging(obs_coords: list[numpy.ndarray], obs_values: list[numpy.ndarray], grid_coord: numpy.ndarray, vgm_spec: dict, nmax: int | None = None, rangescale: float | None = None, localnugget: float | None = None, nthread: int = 0, ncache: int | None = None, std_ck: bool = False) → tuple[numpy.ndarray, numpy.ndarray]#

One-shot ordinary co-kriging with multiple variables.

Parameters:

obs_coords (list of ndarray, each shape (nobs_i, ndim)) – Observation coordinates per variable. Rows are points.
obs_values (list of ndarray, each shape (nobs_i,)) – Observation values per variable.
grid_coord (ndarray, shape (ngrid, ndim)) – Grid coordinates.
vgm_spec (dict) – Mapping (ivar, jvar) to a variogram dict or list of dicts. Each dict is passed as keyword arguments to Kriging.set_vgm(). Both (i,j) and (j,i) can be provided; if only (i,j) is given, (j,i) will mirror it automatically (handled inside Fortran set_vgm).
nmax (int) – Maximum neighbours per variable.
nthread (int) – max OMP threads for this call (0 or absent = OMP default)
ncache (int, optional) – Per-thread hcache slots for this solve. None uses the default.
std_ck (bool) – Use standard Ordinary Kriging.

Returns:

estimate (ndarray, shape (ngrid,))
variance (ndarray, shape (ngrid,))

Example

>>> est, var = cokriging(
...     obs_coords=[coord1, coord2],
...     obs_values=[val1, val2],
...     grid_coord=grid,
...     vgm_spec={
...         (1,1): dict(vtype="sph", nugget=100, sill=900, a_major=1000, a_minor1=500),
...         (2,2): dict(vtype="sph", nugget=50,  sill=450, a_major=1000, a_minor1=500),
...         (1,2): dict(vtype="sph", nugget=0,   sill=600, a_major=1000, a_minor1=500),
...     })

krigekit.sequential_gaussian_simulation(obs_coord: numpy.ndarray, obs_value: numpy.ndarray, grid_coord: numpy.ndarray, vgm_spec: str, nsim: int, nmax: int | None = None, randpath: numpy.ndarray | None = None, sample: numpy.ndarray | None = None, seed: int | None = None, rangescale: float | None = None, localnugget: float | None = None, nthread: int = 0, ncache: int | None = None) → numpy.ndarray#

Sequential Gaussian Simulation.

Parameters:

obs_coord (ndarray, shape (nobs, ndim)) – Observation coordinates. Rows are points, columns are spatial dimensions.
obs_value (ndarray, shape (nobs,)) – Observation values.
grid_coord (ndarray, shape (ngrid, ndim)) – Grid coordinates.
vgm_spec (dict or list of dict) – One or more nested variogram structure dicts, each passed as keyword arguments to Kriging.set_vgm().
nsim (int) – Number of realisations.
nmax (int) – Maximum neighbours (includes previously simulated nodes).
seed (int, optional) – Random seed for reproducibility.
nthread (int) – max OMP threads for this call (0 or absent = OMP default)
ncache (int, optional) – Per-thread hcache slots for this solve. None uses the default.

Returns:

simulations (ndarray, shape (nsim, ngrid)) – Each row is one realisation in the original (non-randomised) block order.

class krigekit.SpaceTimeKriging(nvar: int = 1, ndrift: int = 0, unbias: int = 1, nsim: int = 0, anisotropic_search: bool = False, weight_correction: bool = False, use_old_weight: bool = False, store_weight: bool = False, cross_validation: bool = False, write_mat: bool = False, neglect_error: bool = True, verbose: bool = False, weight_file: str = '', bounds: tuple | None = None, seed: int | None = None)#

Python interface to the Fortran t_kriging_st space-time kriging engine.

Supports 3D spatial + 1D temporal data, sum-metric and product-sum covariance models, ordinary/simple kriging, co-kriging, ST gradient constraints, and SGSIM (primary variable only, conditioned on secondary observations).

Coordinate convention#

Observation coord arrays use (nobs, ndim+1) shape — the first ndim columns are spatial (x, y [, z]) and the last column is time. ndim may be 2 (x, y, t) or 3 (x, y, z, t) and is inferred from the first set_obs() call. set_grid() accepts either the same combined (ngrid, ndim+1) format or split (ngrid, ndim) + time arrays.

Typical workflow (single variable, sum-metric)#

>>> k = SpaceTimeKriging(nvar=1)
>>> k.set_st_model(model='sum_metric', transform='bounded', at=5.0)
>>> k.set_obs(ivar=1, coord=obs_coord_st, value=obs_value,
...           nmax=30, maxdist=5000)
>>> k.set_vgm(ivar=1, jvar=1, vtype="sph", nugget=0, sill=0.8, a_major=1000, a_minor1=500, a_minor2=200)
>>> k.set_vgm_temporal(ivar=1, jvar=1, spec="exp 0 0.6 10.0")
>>> k.set_vgm_joint_sills(ivar=1, jvar=1, sills=[0.4])
>>> k.set_grid(coord=grid_coord, time=grid_time)
>>> k.set_search(ivar=1)
>>> k.solve()
>>> estimate, variance = k.get_results()
>>> del k

set_st_model(model: str = 'sum_metric', transform: str = 'linear', at: float = 1.0, time_nugget: float = 0.0, time_sill: float = 1.0, k_ps: float = 0.0)#

Set global space-time model parameters. Must be called before set_vgm.

Parameters:

model (str) – 'sum_metric' or 'product_sum'.
transform (str) – Variogram type used for f(dt): 'nug' | 'sph' | 'exp' | 'gau' | 'pow' | 'bsq' | 'cir' | 'lin'. Aliases: 'linear' → 'lin', 'bounded' → 'exp', 'power' → 'pow'.
at (float) – Joint temporal scale (same time units as observations).
time_nugget (float) – Nugget jump in f(dt) for the sum-metric model (applied for dt ≠ 0).
time_sill (float) – Upper scale in f(dt): f(dt) = time_nugget + time_sill * (1 - corefunc(|dt| / at)) for dt ≠ 0; f(0) = 0 always.
k_ps (float) – Product-sum coefficient k (model='product_sum' only).

set_obs(ivar: int, coord: numpy.ndarray, value: numpy.ndarray, time: numpy.ndarray | None = None, variance: numpy.ndarray | None = None, nmax: int | None = None, maxdist: float | None = None, sk_mean: float = 0.0)#

Load observations for variable ivar. Duplicate checks include all coordinate columns including time, so repeated spatial locations are allowed at different times.

Accepts two coordinate formats — pick whichever matches your workflow:

Combined format (default):

coord : (nobs, ndim+1) — first ndim columns are spatial (x[,y[,z]]),
        last column is time; ndim must be 1, 2, or 3.
time  : omitted (None)

Split format (explicit time array):

coord : (nobs, ndim)  spatial coordinates only
time  : (nobs,)       observation times

ndim is inferred from the first set_obs() call and must be consistent across all subsequent calls on the same object.

Parameters:

value ((nobs,) observed values)
variance ((nobs,) measurement error variance (default: zeros))
nmax (max neighbours)
maxdist (max search radius in km-equivalent space (same units as h_ST))
sk_mean (global mean for simple kriging (unbias=0); default 0)

update_obs_value(ivar: int, value: numpy.ndarray)#

Replace observation values for variable ivar in-place.

Coordinates and the KD-tree are unchanged. After solving once with stored weights, call this method with new observed values and solve again to reuse the existing neighbourhoods and weights.

set_obs_drift(ivar: int, drift: numpy.ndarray)#: Set external drift values at observations for variable ivar. drift shape: (nobs, ndrift) — transposed internally.

set_grad(coord1: numpy.ndarray, coord2: numpy.ndarray, grad_value: numpy.ndarray, ivar: int = 1, variance: numpy.ndarray | None = None, drift_ext: numpy.ndarray | None = None)#

Set time-aware ST gradient observation pairs.

coord1 and coord2 must both have shape (ngrad, ndim+1) with columns x, y [, z], t (same ndim as set_obs()). The constraint is Z(coord1[i]) - Z(coord2[i]) = grad_value[i]. Because time is part of each endpoint coordinate, targets at other times are penalized by the temporal covariance model.

set_vgm(ivar: int, jvar: int, vtype: str, nugget: float = 0.0, sill: float = 1.0, a_major: float = 1.0, a_minor1: float | None = None, a_minor2: float | None = None, azimuth: float = 0.0, dip: float = 0.0, plunge: float = 0.0)#

Add one spatial nested structure to vgm(ivar, jvar). Call multiple times for nested models.

Same parameters as Kriging.set_vgm() — see that docstring.

set_vgm_temporal(ivar: int, jvar: int, vtype: str, nugget: float = 0.0, sill: float = 1.0, at_k: float = 1.0)#

Add one temporal nested structure to vgm(ivar, jvar). Call multiple times for nested models.

vtype : variogram type (e.g. ‘sph’, ‘exp’, ‘gau’) nugget : nugget contribution of this structure sill : partial sill of this structure at_k : temporal practical range (same time units as observations)

set_vgm_joint_sills(ivar: int, jvar: int, *sills: float)#

Set joint sills for the sum-metric model.

Pass one float per spatial nested structure of vgm(ivar, jvar). Must be called after all set_vgm() calls for (ivar, jvar).

Example:: k.set_vgm_joint_sills(1, 1, 0.05, 0.07)

set_grid(coord: numpy.ndarray, time: numpy.ndarray | None = None, rangescale: numpy.ndarray | None = None, localnugget: numpy.ndarray | None = None)#

Set point estimation targets.

Accepts two coordinate formats — pick whichever matches your workflow:

Combined format (consistent with set_obs()):

coord : (ngrid, ndim+1) — first ndim columns are spatial (x[,y[,z]]),
        last column is time; ndim must be 1, 2, or 3.
time  : omitted (None)

Split format (explicit time array):

coord : (ngrid, ndim)  spatial coordinates only
time  : (ngrid,)       prediction times

Parameters:

rangescale ((ngrid,) local range scale factors (default: ones))
localnugget ((ngrid,) local nugget additions (default: zeros))

set_grid_cv()#: Cross-validation mode: predict at observation locations.

set_grid_drift(drift: numpy.ndarray)#: Drift values at estimation grid. drift shape: (ngrid, ndrift).

set_sim(randpath: numpy.ndarray | None = None, sample: numpy.ndarray | None = None)#: Prepare SGSIM random path and pre-drawn N(0,1) samples. Call after set_grid() and set_obs() but before set_search(). When randpath/sample are None, Fortran generates them internally.

set_search(ivar: int, time_at: float = 1.0, anis1: float = 1.0, anis2: float = 1.0, azimuth: float = 0.0, dip: float = 0.0, plunge: float = 0.0, sector_search: bool = False)#

Build the space-time KD-tree and configure the search ellipse for variable ivar.

Call after set_obs() (and after set_sim() for ivar=1 in SGSIM).

Parameters:

ivar (int) – Variable index (1-based).
time_at (float) – Temporal scale factor (default 1.0) to convert the time axis into km-equivalent search units: the search-tree time coordinate is t * time_at. This ensures that L2 distance in the 4D search space matches the sum-metric space-time distance: h_ST = sqrt(h_S^2 + (time_at * dt)^2). Normally, you should pass the same value as at in set_st_model().
anis1 (float) – Spatial minor/major anisotropy ratio (default 1.0).
anis2 (float) – Spatial vertical/major anisotropy ratio (default 1.0).
azimuth (float) – Azimuth of the spatial major axis in degrees (default 0.0, clockwise from North).
dip (float) – Dip angle of the spatial major axis in degrees (default 0.0, positive downward).
plunge (float) – Plunge angle of the spatial major axis in degrees (default 0.0).
sector_search (bool) – Enable sector (octant) search limiting candidates per sector. If True, candidate neighbours are partitioned into 8 spatial octants centered on the prediction location. At most nmax (from set_obs()) candidates are selected per octant. This ensures a balanced spatial distribution of neighbours and prevents clustering artifacts. The maximum total neighbours selected is 8 * nmax. If search anisotropy is enabled, spatial coordinates are rotated/scaled according to the anisotropy parameters before sector assignment.

solve(nthread: int = 0, ncache: int | None = None)#

Run the ST kriging or SGSIM loop.

Parameters:

nthread (int, optional) – Maximum number of OpenMP threads. 0 (default) lets the OpenMP runtime choose (respects OMP_NUM_THREADS).
ncache (int, optional) – Number of per-thread multi-slot hcache entries for this solve. None keeps the compiled/object default, 0 disables hcache, and 1 builds a one-slot hcache for overhead comparisons.

get_factor() → dict#: Return cached factor matrices and the assembled linear system.

property solver_stats: dict#

Solver statistics from the most recent solve() call.

Returns a dict with three integer counts that are reset to zero at the start of every solve():

chol_ok: Blocks solved by Cholesky factorization (either a fresh factorize or a cache hit that reused a previously computed Cholesky factor).
ssytrf_fact: Number of SSYTRF (Bunch-Kaufman LDL^T) factorizations performed. Each one is O(n³) but occurs only once per unique neighbourhood. A non-zero value means Cholesky failed for at least one neighbourhood; a value equal to 1 with global search means the factorization was done once and cached for all blocks.
ssytrf_reuse: Blocks solved by a cached SSYTRF factorization using SSYTRS, which is O(n²). When this is large relative to ssytrf_fact the SSYTRF caching is working effectively.

Example — global neighbourhood, Cholesky fails, 10 000 grid blocks:

k.solve()
s = k.solver_stats
# Expected: chol_ok=0, ssytrf_fact=1, ssytrf_reuse=9999

get_results(copy: bool = False, squeeze: bool = True) → tuple[np.ndarray, np.ndarray]#

Retrieve kriging estimate and variance.

Parameters:

copy (bool, default False) – If True, return C-contiguous copies for downstream NumPy/Pandas use. If False, return views / Fortran-order arrays when possible.
squeeze (bool, default True) – If True, return a 1-D estimate when nsim == 1.

Returns:

estimate (ndarray, shape (ngrid,) when nsim == 1 and squeeze;) – otherwise shape (ngrid, nsim) [block first]
variance (ndarray, shape (ngrid,))

krigekit.spacetime_kriging(obs_coord: numpy.ndarray, obs_value: numpy.ndarray, grid_coord: numpy.ndarray, grid_time: numpy.ndarray, spatial_spec: dict | list[dict], temporal_spec: dict | list[dict], joint_sills: list[float], model: str = 'sum_metric', transform: str = 'linear', at: float = 1.0, time_nugget: float = 0.0, time_sill: float = 1.0, nmax: int = 20, maxdist: float | None = None, search_anis1: float = 1.0, search_anis2: float = 1.0, search_azimuth: float = 0.0, k_ps: float = 0.0, nthread: int = 0, ncache: int | None = None) → tuple[np.ndarray, np.ndarray]#

One-shot ordinary space-time kriging (single variable).

Parameters:

obs_coord ((nobs, ndim+1) observation coordinates — first ndim cols spatial, last col time)
obs_value ((nobs,) observed values)
grid_coord ((ngrid, ndim) prediction spatial coordinates)
grid_time ((ngrid,) prediction times)
spatial_spec (dict or list[dict] spatial variogram structure(s))
temporal_spec (dict or list[dict] temporal variogram structure(s))
joint_sills (list[float] joint sills (sum-metric only))
model ('sum_metric' or 'product_sum')
transform ('nug', 'sph', 'exp', 'gau', 'pow', 'bsq', 'cir', or 'lin')
at (joint temporal scale (also used as time_at for the KD-tree))
time_nugget (temporal variogram nugget/sill for set_vgm_temporal)
time_sill (temporal variogram nugget/sill for set_vgm_temporal)
nmax (max neighbours)
maxdist (max search radius in km-equivalent space (h_ST units))
nthread (max OMP threads for this call (0 = OMP default))
ncache (per-thread hcache slots for this solve; None uses default)

Returns:

estimate ((ngrid,))
variance ((ngrid,))

krigekit.spacetime_cokriging(obs_coords: list[np.ndarray], obs_values: list[np.ndarray], grid_coord: numpy.ndarray, grid_time: numpy.ndarray, spatial_specs: dict, temporal_specs: dict, joint_sills: dict, model: str = 'sum_metric', transform: str = 'linear', at: float = 1.0, time_nugget: float = 0.0, time_sill: float = 1.0, nmax: int = 20, maxdist: float | None = None, nthread: int = 0, ncache: int | None = None) → tuple[np.ndarray, np.ndarray]#

One-shot ordinary space-time co-kriging.

Parameters:

obs_coords (list of (nobs_i, ndim+1) arrays, one per variable — first ndim cols spatial, last col time)
obs_values (list of (nobs_i,) arrays)
grid_coord ((ngrid, ndim))
grid_time ((ngrid,))
spatial_specs (dict (ivar,jvar) -> dict or list[dict])
temporal_specs (dict (ivar,jvar) -> dict or list[dict])
joint_sills (dict (ivar,jvar) -> list[float])
nthread (max OMP threads for this call (0 = OMP default))
ncache (per-thread hcache slots for this solve; None uses default)

Returns:

estimate ((ngrid,))
variance ((ngrid,))

class krigekit.IndicatorKriging(ncat: int, nvar: int | None = None, ndim: int = 2, ndrift: int = 0, unbias: int = 1, nsim: int = 0, anisotropic_search: bool = False, weight_correction: bool = False, use_old_weight: bool = False, store_weight: bool = False, cross_validation: bool = False, write_mat: bool = False, neglect_error: bool = True, varying_vgm: bool = False, std_ck: bool = False, verbose: bool = False, pf_cache: bool = False, weight_file: str = '', bounds: tuple | None = None, seed: int | None = None)#

Bases: krigekit.kriging.Kriging

Multiple Indicator Kriging and Sequential Indicator Simulation, with optional secondary co-variate support.

Extends Kriging — all setup, solve, and results methods are inherited unchanged. The differences are:

The Fortran object is a t_kriging_indicator (created via krige_ind_create), which overrides prepare, sim_draw, and post_solve to implement indicator-specific behaviour.
ncat names the K indicator categories. nvar defaults to ncat (pure MIS) but can be set larger to add secondary continuous co-variates for co-kriging MIS (see below).
set_categorical_obs() is a convenience method that converts raw category labels into K binary indicator datasets.

Parameters:

ncat (int) – Number of categories K. Indicator variables occupy ivar = 1..ncat.
nvar (int, optional) – Total number of co-kriging variables. Defaults to ncat (pure MIS). Set nvar = ncat + M to include M secondary continuous variables (ivar = ncat+1 .. nvar). Secondary variables contribute to the kriging weights but are excluded from the CDF draw and probability normalisation.
ndim (int) – Number of spatial dimensions (2 or 3).
nsim (int) – 0 = estimation (returns probabilities); >0 = SIS (returns one-hot draws).

:param All other keyword arguments are passed through to Kriging.:

Notes

For SIS (nsim > 0), call set_sim() after set_grid().

Co-kriging MIS example (K=3 categories + 1 secondary variable):

ik = IndicatorKriging(ncat=3, nvar=4, ndim=2)
ik.set_categorical_obs(coord, cats, nmax=20)   # loads ivar=1,2,3
ik.set_obs(ivar=4, coord=sec_coord, value=sec_val)  # secondary
for iv in range(1, 5):
    for jv in range(1, 5):
        ik.set_vgm(ivar=iv, jvar=jv, ...)
ik.set_grid(coord=grid_coord)
for k in range(1, 5):
    ik.set_search(ivar=k)
ik.solve()
probs, var = ik.get_results()   # shape (ngrid, 3) — secondary excluded

set_sim(randpath: numpy.ndarray | None = None, sample: numpy.ndarray | None = None)#

Set up Sequential Indicator Simulation parameters.

Overrides set_sim(). When sample is None, delegates to the parent with no sample so that the Fortran set_sim_indicator override generates U(0, 1) draws directly; no sample array is created in Python. When sample is supplied, validates that every value lies in [0, 1] then delegates to the parent.

Parameters:

randpath (ndarray of int, shape (nblocks,), optional) – Random visiting order (1-based). Generated with a random permutation if omitted.
sample (ndarray, shape (nblocks, nvar, nsim), optional) – Pre-drawn U(0, 1) samples. Every value must lie in [0, 1]. When None, Fortran generates U(0, 1) via set_sim_indicator.

set_categorical_obs(coord: numpy.ndarray, categories: numpy.ndarray, category_labels=None, nmax: int = 32, maxdist: float | None = None, variance: numpy.ndarray | None = None)#

Build K binary indicator datasets from raw category labels.

For each category k (ivar = 1..K), creates a binary indicator array I_k = (categories == category_labels[k-1]) and calls set_obs(ivar=k, ...)().

Parameters:

coord (ndarray, shape (nobs, ndim)) – Sample coordinates.
categories (array-like, shape (nobs,)) – Integer or string category label at each sample.
category_labels (list, optional) – Ordered list of K unique labels — one per indicator variable. If omitted, sorted unique values from categories are used. The order determines which label maps to ivar=1, 2, …, K.
nmax (int) – Maximum neighbours per indicator variable.
maxdist (float) – Maximum search distance (0 = unlimited).
variance (ndarray, shape (nobs,) or None) – Measurement error variance at each sample (0 if omitted).

set_indicator_vgm(vtype: str = 'sph', nugget: float = 0.0, sill: float = 1.0, a_major: float = 1.0, a_minor1: float | None = None, a_minor2: float | None = None, azimuth: float = 0.0, dip: float = 0.0, plunge: float = 0.0, cross: str = 'same', proportions: numpy.ndarray | None = None)#

Set indicator variograms for all K² pairs in one call.

This is a convenience wrapper around set_vgm() that handles the three standard approaches for the cross-variogram sills.

Parameters:

vtype – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
nugget – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
sill – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
a_major – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
a_minor1 – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
a_minor2 – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
azimuth – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
dip – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
plunge – Shared variogram model — same shape and range are used for all pairs. sill is the base partial sill (see cross for how it is applied to off-diagonal pairs).
cross ({"same", "independent", "proportional"}) –
Controls how cross-variogram sills are computed:
- "same" (default) — sill is identical for all K² pairs. Simple and robust; post_solve normalisation absorbs the error.
- "independent" — auto-variogram sill as given (or from proportions); cross-variogram sill = 0 (no co-kriging). Equivalent to running K separate ordinary kriging systems.
- "proportional" — auto sills ∝ p_k (1 − p_k) from proportions; cross sills = √(s_k · s_l). Satisfies the LMC positive-definiteness condition for each nested structure. Requires proportions.
proportions (array-like of shape (K,), optional) –
Observed category proportions [p_1, …, p_K]. When given:
- auto sill for category k = p_k · (1 − p_k) (indicator variance)
- cross sill for k ≠ l = √(s_k · s_l) (if cross="proportional")
Ignored when cross="same".

Notes

Theoretical sills differ per category (p_k varies), and cross-sills are theoretically negative (−p_k · p_l for jointly exhaustive indicators). cross="same" uses one positive value for all pairs, relying on post_solve_indicator normalisation to produce a valid probability simplex. cross="proportional" is more principled; cross="independent" is the most conservative (no coupling between indicator variables).

Examples

Same sill for all 16 pairs (current practice):

ik.set_indicator_vgm(vtype="sph", nugget=0.02, sill=0.19,
                     a_major=500, a_minor1=80, azimuth=90)

Category-specific sills from observed proportions, proportional cross-terms:

props = np.array([0.18, 0.23, 0.21, 0.38])
ik.set_indicator_vgm(vtype="sph", nugget=0.02, sill=0.19,
                     a_major=500, a_minor1=80, azimuth=90,
                     cross="proportional", proportions=props)

Independent kriging (zero cross-covariance):

ik.set_indicator_vgm(vtype="sph", nugget=0.02, sill=0.19,
                     a_major=500, a_minor1=80, azimuth=90,
                     cross="independent", proportions=props)