krigekit#

Python wrapper for a Fortran kriging and SGSIM engine.

PyPI binary wheels include the compiled Fortran library. A source build is needed only for development, custom compiler settings, or unsupported platforms. See the installation guide for details.

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

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,
`SpaceTimeVariogramModel`	Joint model composed from spatial and temporal marginal models.
`VgmComponent`	One nested variogram component.
`VgmStructure`	An ordered set of nested variogram components for one variable pair.
`VgmStructureST`	Theoretical space-time variogram: spatial/temporal structures + coupling.
`ObservationSet`	Observation data and search configuration for one variable.
`VariogramModel`	Python-side variogram model with `Kriging.set_vgm`-style specs.
`VariogramSystem`	Multivariable variogram system for cokriging workflows.
`IndicatorVariogramSystem`	Variogram system for `K` mutually exclusive categorical indicators.

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.
`calc_anisotropic_lag`(lag, *[, anis1, anis2, azimuth, ...])	Calculate engine-compatible equivalent major-axis lag distance.
`calc_lag_vectors`(coord0, coord1[, pairwise])	Return signed lag vectors `coord1 - coord0`.

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)
>>> 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, nmax=20)
>>> 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().

param ndim:

Number of spatial dimensions (2 or 3).

type ndim:

int

param nvar:

Number of variables (1 for ordinary/simple kriging, >1 for cokriging).

type nvar:

int

param ndrift:

Number of external drift functions (0 = no drift).

type ndrift:

int

param unbias:

1 = ordinary kriging (sum-of-weights = 1 constraint); 0 = simple kriging (no constraint, uses sk_mean).

type unbias:

int

param nsim:

Number of simulations. 0 = kriging only; >0 = SGSIM.

type nsim:

int

param anisotropic_search:

Use anisotropic search ellipse for neighbour lookup.

type anisotropic_search:

bool

param weight_correction:

Force kriging weights to be non-negative and sum to 1.

type weight_correction:

bool

param use_old_weight:

Read pre-computed weights from weight_file instead of solving.

type use_old_weight:

bool

param store_weight:

Write computed weights to weight_file while also estimating blocks.

type store_weight:

bool

param cross_validation:

Leave-one-out cross-validation mode.

type cross_validation:

bool

param write_mat:

Write matrix for debugging.

type write_mat:

bool

param neglect_error:

Ignore solver errors and set failed block to NaN instead of aborting.

type neglect_error:

bool

param varying_vgm:

Use a different variogram per estimation block (spatially varying anisotropy). When True, call set_vgm_block() for each block after set_grid(). Defaults to False (single global model).

type varying_vgm:

bool

param std_ck:

Co-kriging unbiasedness formulation (only relevant when nvar > 1 and unbias=1).

False (default) — Isaaks & Srivastava: single combined constraint (Σw₁ + Σw₂ = 1) plus a local-mean correction applied post-solve. Matches the GSLIB/legacy behaviour.
True — standard cokriging: separate per-variable constraints (Σwᵢ = 1 for own variable, Σwⱼ = 0 for others), equivalent to the gstat/ISATIS formulation. Use this to match gstat output.

type std_ck:

bool

param verbose:

Print progress messages.

type verbose:

bool

param pf_cache:

Enable the persistent between-solve factorization cache. When True, the Cholesky factor of K is stored after the first solve() and reused on subsequent calls when the neighbour set and variogram are unchanged (speeds up repeated solves on the same observation grid). Defaults to False; enable only when you plan to call solve() multiple times and need the speedup.

type pf_cache:

bool

param weight_file:

Path to the weight file (required when use_old_weight or store_weight).

type weight_file:

str

param bounds:

(lower, upper) clipping bounds for the estimate. None means no clipping (uses Fortran defaults: [-huge, +huge]).

type bounds:

tuple(float, float) or None

param seed:

Random seed.

type seed:

int, optional

set_obs(ivar: int, coord: numpy.ndarray, value: numpy.ndarray, variance: numpy.ndarray | 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).
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_nscore(ivar: int = 1, zmin: float | None = None, zmax: float | None = None, ltail='linear', utail='linear', ltpar: float = 1.0, utpar: float = 1.0, weights: numpy.ndarray | None = None)#

Enable the normal-score transform for variable ivar.

Builds a normal-score (Gaussian anamorphosis) transform from the current observation values and replaces them with their normal scores, so that the subsequent solve() runs in Gaussian space; estimates or simulated realisations are automatically back-transformed to data units. The transform lives in the Fortran engine, so every C-API client gets identical, reproducible results.

Call after set_obs() (and before solve()). The variogram supplied with set_vgm() should be that of the normal scores (unit sill). For kriging (nsim == 0), estimates and variances are back-transformed to data units with Gaussian quadrature.

Parameters:

ivar (int) – Variable index, 1-based.
zmin (float, optional) – Lower/upper bounds used by the back-transform tail extrapolation. Default to the minimum/maximum of the observed data (no extrapolation beyond the data range).
zmax (float, optional) – Lower/upper bounds used by the back-transform tail extrapolation. Default to the minimum/maximum of the observed data (no extrapolation beyond the data range).
ltail ({"linear", "power", "hyperbolic"} or int) – Tail-extrapolation model below z[0] and above z[-1] (1=linear, 2=power, 4=hyperbolic; hyperbolic is upper-tail only and requires positive data).
utail ({"linear", "power", "hyperbolic"} or int) – Tail-extrapolation model below z[0] and above z[-1] (1=linear, 2=power, 4=hyperbolic; hyperbolic is upper-tail only and requires positive data).
ltpar (float) – Power / hyperbolic tail parameter (ignored for linear tails).
utpar (float) – Power / hyperbolic tail parameter (ignored for linear tails).
weights (ndarray, shape (nobs,), optional) – Declustering weights for the empirical CDF. Default: equal weights.

set_uscore(ivar: int = 1, zmin: float | None = None, zmax: float | None = None, ltail='linear', utail='linear', ltpar: float = 1.0, utpar: float = 1.0, weights: numpy.ndarray | None = None)#

Enable the uniform quantile transform for variable ivar.

Builds a weighted empirical CDF from the current observation values and replaces those values with their cumulative probabilities in [0, 1]. The subsequent solve() runs in uniform-score space; estimates or simulated values are automatically back-transformed through the empirical quantile table to data units.

Call after set_obs() (and before solve()). The variogram supplied with set_vgm() should be fit on the uniform scores. For kriging (nsim == 0), estimates and variances are back-transformed to data units with Gaussian quadrature.

Parameters are the same as set_nscore(): zmin/zmax bound the back-transform tails, ltail/utail choose tail models, and weights supplies optional declustering weights.

set_quantile(*args, **kwargs)#: Alias for set_uscore().

transform_value_to_score(value, ivar: int = 1)#

Transform data-unit value through the active score transform.

Call set_nscore() or set_uscore() first. The returned scores are normal scores for set_nscore and uniform CDF scores for set_uscore.

transform_score_to_value(score, ivar: int = 1)#

Back-transform score-space score through the active transform.

Call set_nscore() or set_uscore() first. Scores are interpreted according to the active transform for ivar.

back_transform_score(score, ivar: int = 1)#: Alias for transform_score_to_value().

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, nmax: int | None = None, maxdist: float | None = None, 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 of the major axis below horizontal (degrees, positive downward).
plunge (float) – Plunge angle (degrees).
nmax (int, optional) – Maximum number of neighbours. Default: use all observations.
maxdist (float, optional) – Maximum search distance. Default: unlimited.
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 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) – Total size of the shared factor cache pool (slots across all threads). None keeps the compiled/object default. 0 disables factorization reuse entirely — both the shared hcache and the single-entry ctx%cache adjacency cache are switched off, so every block is factorized afresh (chol_reuse/ssytrf_reuse stay zero); use this for a clean no-cache baseline. Positive values below 4 are promoted to 4 (the minimum meaningful pool is 4 slots = 1 bucket). The optional persistent factor cache (pf_cache, set on the object) is controlled separately.

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: observation groups (variable 1..nvar)
- indices nvar..2*nvar-1: simulation groups (SGSIM only)
- indices ngroups_base..ngroups-1: gradient 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) – Total shared hcache pool 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) – Total shared hcache pool 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) – Total shared hcache pool 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)
>>> 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, nmax=30, maxdist=5000)
>>> 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, 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))
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, product: bool = False)#

Add one spatial marginal structure to vgm(ivar, jvar).

Parameters match Kriging.set_vgm(), including product. When product=True, this structure is multiplied with the immediately preceding spatial structure in covariance space.

Call repeatedly for nested models. The first structure cannot be a product member.

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

Add one temporal marginal structure.

Parameters:

ivar (int) – One-based variable-pair indices.
jvar (int) – One-based variable-pair indices.
vtype (str) – Variogram type such as "sph", "exp", "gau" or "hol".
nugget (float, optional) – Nugget contribution of this structure.
sill (float, optional) – Partial sill of this structure.
at_k (float, optional) – Temporal practical range or period parameter, in the same units as the observation time coordinate.
product (bool, optional) – Multiply this structure with the immediately preceding temporal structure in covariance space instead of adding it.

Notes

Call this method repeatedly for nested models. The first structure cannot be a product member.

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, nmax: int | None = None, maxdist: float | None = None, 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 below horizontal, degrees positive downward (default 0.0).
plunge (float) – Plunge angle of the spatial major axis in degrees (default 0.0).
nmax (int, optional) – Maximum number of neighbours. Default: use all observations.
maxdist (float, optional) – Maximum search radius in km-equivalent space (same units as h_ST).
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 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) – Total size of the shared factor cache pool (slots across all threads). None keeps the compiled/object default. 0 disables factorization reuse entirely — both the shared hcache and the single-slot adjacency cache are switched off, so every block is factorized afresh (chol_reuse/ssytrf_reuse stay zero); use this for a clean no-cache baseline. Positive values below 4 are promoted to the four-slot minimum pool. The persistent factor cache (pf_cache, set on the object) is controlled separately.

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 (total shared hcache pool 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 (total shared hcache pool 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).
Indicator encoding and the K x K coregionalization are built with IndicatorVariogramSystem and transferred via its apply(); this engine wrapper no longer owns that construction.

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).
**kwargs – 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):

system = IndicatorVariogramSystem(categories=[1, 2, 3])
system.set_categorical_obs(coord, cats)
system.set_indicator_vgm(vtype="sph", a_major=1000,
                         sill_strategy="theoretical",
                         cross_strategy="closure")
ik = IndicatorKriging(ncat=3, nvar=4, ndim=2)
system.apply(ik)                            # indicator block (ivar 1..3)
ik.set_obs(ivar=4, coord=sec_coord, value=sec_val)   # secondary
ik.set_vgm(ivar=4, jvar=4, ...)             # secondary auto/cross models
ik.set_grid(coord=grid_coord)
for k in range(1, 5):
    ik.set_search(ivar=k, nmax=20)
ik.solve()
probs, var = ik.get_results()   # shape (ngrid, 3) — secondary excluded

Parameters:

ndim (int) – Number of spatial dimensions (2 or 3).
nvar (int) – Number of variables (1 for ordinary/simple kriging, >1 for cokriging).
ndrift (int) – Number of external drift functions (0 = no drift).
unbias (int) – 1 = ordinary kriging (sum-of-weights = 1 constraint); 0 = simple kriging (no constraint, uses sk_mean).
nsim (int) – Number of simulations. 0 = kriging only; >0 = SGSIM.
anisotropic_search (bool) – Use anisotropic search ellipse for neighbour lookup.
weight_correction (bool) – Force kriging weights to be non-negative and sum to 1.
use_old_weight (bool) – Read pre-computed weights from weight_file instead of solving.
store_weight (bool) – Write computed weights to weight_file while also estimating blocks.
cross_validation (bool) – Leave-one-out cross-validation mode.
write_mat (bool) – Write matrix for debugging.
neglect_error (bool) – Ignore solver errors and set failed block to NaN instead of aborting.
varying_vgm (bool) – Use a different variogram per estimation block (spatially varying anisotropy). When True, call set_vgm_block() for each block after set_grid(). Defaults to False (single global model).
std_ck (bool) –
Co-kriging unbiasedness formulation (only relevant when nvar > 1 and unbias=1).
- False (default) — Isaaks & Srivastava: single combined constraint (Σw₁ + Σw₂ = 1) plus a local-mean correction applied post-solve. Matches the GSLIB/legacy behaviour.
- True — standard cokriging: separate per-variable constraints (Σwᵢ = 1 for own variable, Σwⱼ = 0 for others), equivalent to the gstat/ISATIS formulation. Use this to match gstat output.
verbose (bool) – Print progress messages.
pf_cache (bool) – Enable the persistent between-solve factorization cache. When True, the Cholesky factor of K is stored after the first solve() and reused on subsequent calls when the neighbour set and variogram are unchanged (speeds up repeated solves on the same observation grid). Defaults to False; enable only when you plan to call solve() multiple times and need the speedup.
weight_file (str) – Path to the weight file (required when use_old_weight or store_weight).
bounds (tuple(float, float) or None) – (lower, upper) clipping bounds for the estimate. None means no clipping (uses Fortran defaults: [-huge, +huge]).
seed (int, optional) – Random seed.

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

Return indicator results, excluding secondary covariate channels.

The kriging engine stores all nvar estimates internally because secondary variables participate in the cokriging system. Public indicator results contain only the first ncat variables, matching the probability or one-hot category array documented by this class.

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.

class krigekit.SpaceTimeVariogramModel(spatial=None, temporal=None)#

Bases: krigekit.variogram_base._VariogramModelBase

Joint model composed from spatial and temporal marginal models.

This mirrors the Fortran vgm_struct_st layout: spatial corresponds to cs, temporal corresponds to ct, and this object owns only full space-time observations, coupling parameters, and transfer metadata.

Create a space-time model from optional marginal models.

set_vgm(*args, **kwargs)#: Add a structure to the spatial marginal and return self.

set_vgm_temporal(vtype, nugget=0.0, sill=1.0, at_k=1.0, product=False, **kwargs)#: Add a structure to the temporal marginal and return self.

experimental(store=True, **kwargs)#: Calculate the full space-time cloud with stored spatial anisotropy.

calc_spacetime_variogram(spatial_lag, temporal_lag, params=None)#

Evaluate the fitted product-sum space-time semivariogram.

Parameters:

spatial_lag (array-like) – Broadcastable spatial and temporal lag arrays.
temporal_lag (array-like) – Broadcastable spatial and temporal lag arrays.
params (array-like, optional) – (a, b, p, spatial_range, temporal_range). If omitted, use spacetime_params_ from fit_spacetime_product_sum() or set_spacetime_params().

calc_spacetime_variogram_between(coord0, coord1, time0, time1, *, pairwise=False, params=None)#

Evaluate the space-time variogram between coordinates.

Spatial lag vectors are rotated and scaled using the anisotropy stored by set_spacetime_anisotropy(). The resulting distance is the equivalent lag along the major axis, in the original coordinate units.

set_spacetime_anisotropy(*, anis1=1.0, anis2=1.0, azimuth=0.0, dip=0.0, plunge=0.0)#

Set spatial anisotropy for product-sum fitting and evaluation.

anis1 and anis2 are minor/major range ratios, matching the Fortran engine. Angles use the engine convention: azimuth clockwise from north and dip positive downward.

set_spacetime_params(params, *, spatial_vtype=None, temporal_vtype=None)#

Manually set product-sum space-time parameters.

params is (a, b, p, spatial_range, temporal_range). Valid covariance conversion requires p <= 0, a + p > 0 and b + p > 0.

fit_spacetime_product_sum(avgvgm=None, *, spatial_vtype='sph', temporal_vtype='gau', starts=None, bounds=None, spatial_col=None, temporal_col=('time_lag', 'mean'), variogram_col=('variogram', 'mean'), count_col=('variogram', 'count'), weight_cap_quantile=0.9, min_marginal_sill=0.0001, options=None, avg_kwargs=None)#

Fit a constrained product-sum model to averaged space-time bins.

The fitted form is

a*g_s(h_s) + b*g_t(h_t) + p*g_s(h_s)*g_t(h_t),

where both marginal variograms have unit sill. Multiple starting points are fitted with SLSQP; the successful result with the smallest weighted objective is stored.

Returns:: VariogramModel – self with spacetime_params_, spacetime_fit_result_ and spacetime_fit_results_ populated.

calc_spacetime_sum_metric_variogram(spatial_lag, temporal_lag, params=None)#

Evaluate a fitted sum-metric space-time semivariogram.

params contains spatial_scale, temporal_scale, one joint sill per spatial structure, and the joint temporal scale at.

fit_spacetime_sum_metric(spatial_model=None, temporal_model=None, avgvgm=None, *, transform='lin', time_nugget=0.0, time_sill=1.0, p0=None, bounds=None, spatial_col=('distance', 'mean'), temporal_col=('time_lag', 'mean'), variogram_col=('variogram', 'mean'), count_col=('variogram', 'count'), weight_cap_quantile=0.9, max_nfev=20000, avg_kwargs=None)#

Fit sum-metric coupling while retaining fitted marginal shapes.

The fit estimates a spatial marginal scale, temporal marginal scale, one non-negative joint sill per spatial structure, and the joint temporal scale at. Allowing the two marginal scales to adjust avoids forcing separately fitted boundary marginals to explain the entire interior space-time lag surface.

fit(avgvgm=None, *, model='product_sum', **kwargs)#

Fit a space-time coupling model, returning a FitResult.

model selects the coupling form:

"product_sum" fits the constrained a*g_s(h_s) + b*g_t(h_t) + p*g_s(h_s)*g_t(h_t) model; its parameters are (a, b, p, spatial_range, temporal_range).
"sum_metric" fits a spatial marginal scale, a temporal marginal scale, one joint sill per spatial structure, and the joint temporal scale at.

Remaining keyword arguments are forwarded to the underlying fitter, and the fitted parameters are stored on the model for to_spacetime_kriging_specs() / to_sum_metric_kriging_specs(). FitResult.summary() reports the labelled parameter table; variance and p-values are not estimated for these constrained/weighted joint fits, so those columns are NaN.

to_sum_metric_kriging_specs()#: Return fitted sum-metric marginal, coupling, and search parameters.

to_spacetime_kriging_specs(*, z_scale=None, spatial_nugget=0.0, temporal_nugget=0.0)#: Convert fitted product-sum parameters to engine-ready dictionaries.

class krigekit.VgmComponent#

One nested variogram component.

The public representation deliberately matches the flat arguments accepted by krigekit.Kriging.set_vgm(). The Fortran engine groups the anisotropy fields into its internal vgm_aniso type after transfer.

validate()#: Validate finite parameters and strictly positive ranges.

copy(**changes)#: Return a validated copy, optionally replacing selected fields.

property display_name#

the name if set, otherwise the model type.

Type:: Human-readable label

set_anisotropy(*, a_major=None, a_minor1=None, a_minor2=None, ratio_minor1=None, ratio_minor2=None, anis1=None, anis2=None, azimuth=None, dip=None, plunge=None)#

Update ranges and angles in place.

anis1 and anis2 are aliases for the minor/major ratios used by the kriging search API.

anisotropy_dict()#: Return ranges, ratios, and angles used by geometry helpers.

calc_anisotropic_distance(lag)#: Return equivalent major-axis distance for lag vector(s).

calc_covariance(distance)#: Evaluate this component’s covariance at scalar lag distance(s).

calc_covariance_lag(lag)#: Evaluate covariance for coordinate lag vector(s), with anisotropy.

property cov0#: Covariance at zero lag, including this component’s nugget.

calc_variogram(distance)#: Evaluate gamma(h) = C(0) - C(h) at scalar lag distance(s).

calc_variogram_lag(lag)#: Evaluate the semivariogram for coordinate lag vector(s).

to_flat_dict()#: Return the flat engine representation (excludes the name metadata).

classmethod from_flat_dict(spec)#: Construct a component from a strict flat engine-style mapping.

class krigekit.VgmStructure(components=None, name=None)#

An ordered set of nested variogram components for one variable pair.

Create a structure from components or flat component specs.

property ncomponent#: Number of nested components (always len(self.components)).

copy()#: Return an independent copy with copied components.

clear()#: Remove all components and return self.

validate()#: Validate every component and return self.

set_vgm(vtype, nugget=0.0, sill=1.0, a_major=1.0, a_minor1=None, a_minor2=None, azimuth=0.0, dip=0.0, plunge=0.0, append=True, product=False, name=None)#

Add one nested component.

Parameters mirror krigekit.Kriging.set_vgm() (without ivar and jvar). Pass append=False to clear existing components first, or product=True to multiply this component with the preceding one in covariance space.

set_structure_params(index=0, **params)#

Update fields on one component, with validation.

index is a zero-based component position. Accepts any VgmComponent field; the component is rebuilt and revalidated.

set_anisotropy(index=None, **params)#

Update anisotropy on one, several, or all components.

index accepts None (all components), one zero-based integer, or a sequence of integers. Keyword arguments are forwarded to VgmComponent.set_anisotropy() (a_minor1, ratio_minor1 / anis1, azimuth, …).

covariance(distance)#: Evaluate the nested/product covariance at scalar lag distance(s).

property cov0#: Covariance at zero lag, including nugget and product groups.

variogram(distance)#: Evaluate gamma(h) = C(0) - C(h) at scalar lag distance(s).

calc_covariance(coord0, coord1, pairwise=False)#

Evaluate covariance between coordinates, applying anisotropy.

coord0 and coord1 have shape (dim,) or (n, dim). By default matching rows are compared and a single point is broadcast; pairwise=True returns the full (n0, n1) matrix.

calc_variogram(coord0, coord1, pairwise=False)#: Evaluate semivariogram values between coordinates with anisotropy.

fit(data, *, kind='auto', p0=None, x_col=('distance', 'mean'), y_col=('variogram', 'mean'), sigma_col=None, weight_col=None, weights=None, bounds=None, fit_nugget=True, inplace=True, **kwargs)#

Fit sills, major ranges, and nugget to an averaged variogram.

Returns a krigekit.variogram_fitting.FitResult. Fitted values are written into copies of the current components, so model type, anisotropy, and name are preserved. kind selects the fit; only the "isotropic" path (default "auto" for a distance/variogram table) is implemented here – use VariogramModel.fit_anisotropy for directional fits.

to_kriging_specs(replace=False)#

Return flat specs accepted by Kriging.set_vgm.

With replace=True the first spec carries append=False so it clears any existing model for the pair; later specs append.

to_temporal_specs()#

Return specs accepted by SpaceTimeKriging.set_vgm_temporal.

The one-dimensional a_major value is renamed to at_k; spatial anisotropy fields are omitted because temporal marginals are 1-D.

apply_to(kriging, ivar, jvar, replace=True)#: Apply this structure to a krigekit.Kriging object.

class krigekit.VgmStructureST(spatial=None, temporal=None, *, name=None)#

Theoretical space-time variogram: spatial/temporal structures + coupling.

The active coupling form is named by model ("product_sum" or "sum_metric"). Use set_product_sum() or set_sum_metric() to populate it, then calc_variogram() to evaluate and to_kriging_specs() to emit engine-ready dictionaries.

Create an empty structure or wrap existing marginal structures.

property cs#: Alias for the spatial structure, matching vgm_struct_st%cs.

property ct#: Alias for the temporal structure, matching vgm_struct_st%ct.

property ncomponent_spatial#: Number of nested spatial components.

property ncomponent_temporal#: Number of nested temporal components.

copy()#: Return an independent copy of this structure and its parameters.

set_product_sum(params, *, spatial_vtype, temporal_vtype, anisotropy=None)#

Populate a product-sum coupling and rebuild cs/ct.

params is (a, b, p, spatial_range, temporal_range) with unit-sill marginals; a valid covariance conversion requires p <= 0, a + p > 0 and b + p > 0. cs/ct are rebuilt as single-component structures whose sills are a + p and b + p.

set_sum_metric(spatial, temporal, params, *, transform, time_nugget=0.0, time_sill=1.0)#

Populate a sum-metric coupling from fitted marginal structures.

spatial and temporal are VgmStructure marginal shapes, and params is (spatial_scale, temporal_scale, *joint_sills, at) with one joint sill per spatial component.

validate()#: Validate the active coupling parameters and return self.

calc_variogram(spatial_lag, temporal_lag)#: Evaluate the active space-time semivariogram at lag arrays.

to_kriging_specs(**kwargs)#: Return engine-ready specs for the active coupling model.

class krigekit.ObservationSet#

Observation data and search configuration for one variable.

property configured#: True once valid coordinates and values are present.

property nobs#: Number of observation rows (zero before set()).

property ndim#: Spatial dimension (zero before set()).

property ndrift#: Number of external-drift columns (zero when drift is absent).

set(coord, value, times=None, variance=None, sk_mean=None, drift=None)#

Set observation data, validating shapes.

coord/value are required. times, variance and drift are per-observation and replace any previous values (None clears them). sk_mean is updated only when provided, so it survives a data-only re-set().

set_search(nmax=None, maxdist=None, anis1=1.0, anis2=1.0, azimuth=0.0, dip=0.0, plunge=0.0, sector_search=False, time_at=None)#: Configure neighbour search (stored for transfer to the engine).

clear()#: Reset all data and search configuration to defaults.

validate(ndim=None)#: Validate configured data; optionally require a spatial dimension.

is_collocated_with(other)#: Return true when two variables share coordinates and times.

apply_to(kriging, ivar)#: Transfer observations (and drift) to a kriging object for ivar.

class krigekit.VariogramModel(structures=None)#

Bases: krigekit.variogram_base._VariogramModelBase

Python-side variogram model with Kriging.set_vgm-style specs.

The class is useful for fitting, plotting, and validating model choices before replaying the same structures into the Fortran solver. It supports additive nested structures and product structures using the same convention as krigekit.Kriging.set_vgm(): a structure with product=True is multiplied with the immediately preceding structure in covariance space.

Notes

covariance(h) and variogram(h) evaluate lag distances directly. Use calc_covariance() and calc_variogram() to evaluate between coordinates with each structure’s anisotropy parameters applied.

Create an empty marginal model or load set_vgm specifications.

calc_directional_average(rawvgm=None, store: bool = True, raw_kwargs=None, h_width=None, h_bins=15, cutoff=None, bandwidth=None, angle_tol=22.5, robust=False, include_minor2: bool = None, **kwargs)#

Average the raw cloud along the model’s fixed anisotropy axes.

The model orientation is taken from the stored structures. For 2D data, the output contains major and minor1 directions; for 3D data, minor2 is included unless include_minor2=False.

calc_params(fit_nugget: bool = True)#: Return the current flat (sill, range, ..., [nugget]) vector.

set_params(params=None, *, sills=None, ranges=None, sill=None, a_major=None, range_=None, nugget=None, fit_nugget: bool = True)#

Manually update fitted variogram parameters.

params uses the same flat convention returned by fit(): sill0, range0, sill1, range1, ..., [nugget]. For a single structure, the singular keywords sill, a_major/range_ and nugget are convenient alternatives. sills and ranges may be sequences with one value per stored structure.

Returns:: VariogramModel – self, so manual adjustments can be chained before plotting or applying to kriging.

fit(avgvgm=None, p0=None, x_col=('distance', 'mean'), y_col=('variogram', 'mean'), sigma_col=None, weight_col=None, weights=None, bounds=None, inplace: bool = False, makeplot: bool = False, fit_nugget: bool = True, raw_kwargs=None, avg_kwargs=None, **kwargs)#

Fit this model to an averaged variogram, returning a FitResult.

If avgvgm is omitted, the cached avg_variogram_ is used when available, otherwise it is computed from observations supplied by set_obs(). The theoretical fit delegates to VgmStructure.fit(); this wrapper adds the empirical-data convenience and updates the analysis fit-state caches.

Returns:: FitResult – .target is this model (inplace=True) or a new fitted model; .params/.cov/.metrics carry the fit outputs and .summary() gives a labelled table.

fit_aniso_angle(rawvgm=None, n_struct=None, set_ranges=True, raw_kwargs=None, **kwargs)#

Estimate the anisotropy orientation from the empirical cloud and apply it.

This is the first fitting step of the anisotropic workflow: run it before fit_anisotropy() so the directional binning and the sill/range fit use the correct axes. A 3-D cloud uses the multi-started model-based profile fit (krigekit.fit_aniso_angle()); a 2-D cloud falls back to the fast PCA azimuth (estimate_aniso_angle()). The fitted azimuth / dip / plunge are written into every structure, and with set_ranges=True the minor ranges are seeded from the fitted anisotropy ratios.

rawvgm defaults to the cached cloud, otherwise it is computed from observations. n_struct defaults to the number of structures. Returns self so it can precede the rest of the workflow.

calc_anisotropic_params(include_minor2: bool = False, fit_nugget: bool = True)#: Return the current anisotropic flat parameter vector.

set_anisotropic_params(params, include_minor2: bool = False, fit_nugget: bool = True)#

Manually update sills and anisotropic ranges from a flat vector.

params follows the fit_anisotropy() convention: sill, a_major, a_minor1 for each structure, optionally a_minor2, followed by a trailing nugget when fit_nugget=True.

fit_anisotropy(directional=None, p0=None, axis_col='direction', x_col='lag', y_col='variogram', sigma_col=None, weight_col=None, weights=None, bounds=None, inplace: bool = False, makeplot: bool = False, fit_nugget: bool = True, include_minor2: bool = None, raw_kwargs=None, directional_kwargs=None, maxfev=9999, ax=None, xlabel='Lag', ylabel='Semivariogram', **kwargs)#

Fit sills and anisotropic ranges with fixed orientation.

The fitted parameter vector is sill, a_major, a_minor1 for each structure, optionally a_minor2 for 3D fits, followed by a trailing nugget when fit_nugget=True. Directional data can be supplied directly, or it is computed with calc_directional_average() using the current azimuth, dip and plunge as fixed axes.

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

Add one nested variogram structure.

Parameters mirror krigekit.Kriging.set_vgm(), except ivar and jvar are omitted because this object represents one variable-pair model. Pass append=False to clear existing structures before adding the new one. Pass product=True to multiply this structure with the immediately preceding structure in covariance space.

The stored ranges and rotation angles define theoretical-model evaluation, directional axes, fitting, plotting overlays, and transfer to kriging. They do not implicitly transform a raw empirical cloud; pass anisotropy=... to calc_experimental() when desired. This separation is necessary because nested structures may have different anisotropy parameters.

Returns:: VariogramModel – self, so calls can be chained.

set_structure_params(index: int = 0, **params)#

Manually update fields on one stored variogram structure.

Parameters:

index (int, optional) – Zero-based structure index.
**params – Any VgmComponent field except append. Use this for edits that do not fit in the flat set_params vector, such as a_minor1, azimuth, dip or product.

set_anisotropy(structures=None, *, a_minor1=None, a_minor2=None, ratio_minor1=None, ratio_minor2=None, anis1=None, anis2=None, azimuth=None, dip=None, plunge=None)#

Apply anisotropy parameters to one or more structures.

Parameters:

structures (int or sequence of int, optional) – Zero-based structure indices to update. Defaults to all stored structures.
a_minor1 (float or sequence, optional) – Absolute minor-axis ranges. A sequence must have one value per selected structure.
a_minor2 (float or sequence, optional) – Absolute minor-axis ranges. A sequence must have one value per selected structure.
ratio_minor1 (float or sequence, optional) – Minor/major range ratios. anis1 and anis2 are accepted as aliases for compatibility with the kriging search terminology.
ratio_minor2 (float or sequence, optional) – Minor/major range ratios. anis1 and anis2 are accepted as aliases for compatibility with the kriging search terminology.
azimuth (float, optional) – Rotation angles in degrees.
dip (float, optional) – Rotation angles in degrees.
plunge (float, optional) – Rotation angles in degrees.

Returns:

VariogramModel – self.

covariance(h)#

Evaluate the nested/product covariance model at lag distance h.

Delegates to VgmStructure.covariance(); product groups are evaluated exactly like the Fortran engine.

property cov0#: Covariance at zero lag, including nugget and product groups.

variogram(h)#: Evaluate the semivariogram gamma(h) = C(0) - C(h).

calc_covariance(coord0, coord1, pairwise: bool = False)#

Evaluate covariance between coordinates, applying anisotropy.

Parameters:

coord0 (array-like) – Coordinates with shape (dim,) or (n, dim). By default, arrays with matching lengths are compared row-wise; if one side has one coordinate it is broadcast against the other side.
coord1 (array-like) – Coordinates with shape (dim,) or (n, dim). By default, arrays with matching lengths are compared row-wise; if one side has one coordinate it is broadcast against the other side.
pairwise (bool, optional) – If true, return the full (n0, n1) covariance matrix.

Returns:

numpy.ndarray or scalar-like – Covariance value(s) from the nested/product model.

calc_variogram(coord0, coord1, pairwise: bool = False)#: Evaluate semivariogram values between coordinates with anisotropy.

plot(avgvgm=None, ax=None, x_col=('distance', 'mean'), y_col=('variogram', 'mean'), h=None, plot_data: bool = True, plot_model: bool = True, annotate: bool = True, plotkws_data=None, plotkws_model=None, xlabel='Lag', ylabel='Semivariogram')#

Plot cached/explicit averaged data and the current model curve.

If avgvgm is omitted, the cached _avg table from calc_average() is used when available. The curve is evaluated with variogram(), so it represents the isotropic lag-distance model. Use calc_variogram() for anisotropy-aware coordinate evaluation.

plot_map(rawvgm=None, ax=None, angle_aniso='model', ellipse_aniso='model', estimate: bool = False, raw_kwargs=None, **kwargs)#

Plot a 2D variogram map from the cached raw cloud.

Parameters:

rawvgm (pandas.DataFrame, optional) – Raw variogram cloud. If omitted, the cached _raw table is used, or computed from stored observations.
angle_aniso ({"model", "estimate", None} or float, optional) – Angle overlay for maximum continuity. "model" uses the first structure’s azimuth, "estimate" estimates the angle from the raw cloud, and a float uses that azimuth directly.
ellipse_aniso ({"model", None} or tuple, optional) – Ellipse overlay. "model" uses the first structure’s major and first minor range.
estimate (bool, optional) – Shorthand for angle_aniso="estimate".
raw_kwargs (dict, optional) – Keyword arguments passed to experimental() if a new raw cloud must be computed.
**kwargs – Forwarded to plot_vgm_map().

Returns:

matplotlib.axes.Axes – Axis containing the variogram map.

plot_map3d(rawvgm=None, ax=None, angle_aniso='model', estimate: bool = False, raw_kwargs=None, **kwargs)#

Plot a 3D variogram map as orthogonal fence sections.

Calls plot_vgm_map3d(). By default (rotate_fences=False) fences align with the world X/Y/Z axes:

Fence A â€” horizontal XY plane (azimuth pattern).
Fence B (n_fences â‰¥ 2) â€” vertical XZ Eastâ€“West section (dip).
Fence C (n_fences â‰¥ 3) â€” vertical YZ Northâ€“South section.

A red line is projected onto each fence showing the major axis direction so the fitted orientation can be compared with the empirical map. Pass rotate_fences=True to rotate the fences to the model’s principal planes instead.

Parameters:

rawvgm (pandas.DataFrame, optional) – Raw 3D variogram cloud. If omitted the cached cloud is used, or computed from stored observations.
angle_aniso ({"model", "estimate", None} or float or tuple, optional) – Model orientation. "model" reads (azimuth, dip, plunge) from the first fitted structure; "estimate" estimates the major direction from the raw cloud; a float is azimuth only; a tuple is (azimuth[, dip[, plunge]]).
estimate (bool, optional) – Shorthand for angle_aniso="estimate".
raw_kwargs (dict, optional) – Forwarded to calc_experimental() when the raw cloud must be computed.
**kwargs – Forwarded to plot_vgm_map3d() (e.g. n_fences, rotate_fences, dx, dz, cutoff, vmax, fill_nan).

Returns:

matplotlib.axes.Axes – 3D axis containing the variogram map.

to_kriging_specs(replace: bool = False)#

Return structures as dictionaries accepted by Kriging.set_vgm.

Parameters:: replace (bool, optional) – If true, the first returned spec has append=False and later specs have append=True. This is convenient when applying a complete model to a reused krigekit.Kriging object.

apply_to(kriging, ivar: int, jvar: int, replace: bool = True)#

Apply this model to a krigekit.Kriging object.

The first structure clears any existing model for (ivar, jvar) when replace=True. Set replace=False to append all structures to an existing model.

to_temporal_specs()#

Return structures accepted by SpaceTimeKriging.set_vgm_temporal.

The one-dimensional a_major value is renamed to at_k. Spatial anisotropy fields are intentionally omitted because temporal marginal structures are one-dimensional.

apply_temporal_to(kriging, ivar: int, jvar: int)#

Append this model to a SpaceTimeKriging temporal marginal.

The target pair should not already contain temporal structures. The space-time API currently resets spatial and temporal marginals together, so this helper deliberately does not offer a replace mode.

class krigekit.VariogramSystem(nvar=None)#

Multivariable variogram system for cokriging workflows.

The system stores observations and variogram structures by 1-based variable pair (ivar, jvar). Each pair is a VgmStructure, while fit_lmc() fits all requested pairs together with positive-semidefinite coregionalization matrices.

Create an empty multivariable variogram system.

nvar=None selects dynamic mode: the variable count grows to the largest referenced 1-based index. An explicit nvar is a strict upper bound and access beyond it raises.

set_obs(ivar, coord, value, times=None, variance=None, sk_mean=None, drift=None)#

Store observations for variable ivar (ergonomic wrapper).

Delegates to self.obs[ivar].set(...) and invalidates cached empirical variograms for any pair containing ivar.

set_search(ivar, *, nmax=None, maxdist=None, anis1=1.0, anis2=1.0, azimuth=0.0, dip=0.0, plunge=0.0, sector_search=False, time_at=None)#

Configure per-variable neighbour search, transferred by apply().

Search settings are stored on the variable’s ObservationSet; the Fortran engine performs the actual neighbour selection.

set_vgm(ivar, jvar, vtype, **kwargs)#: Add one nested structure to the model for (ivar, jvar).

set_raw_vgm(ivar, jvar, rawvgm)#: Store an externally computed raw variogram cloud for a pair.

set_avg_vgm(ivar, jvar, avgvgm)#: Store an externally computed averaged variogram for a pair.

calc_experimental(ivar=None, jvar=None, cross='auto', store=True, **kwargs)#

Compute a raw empirical variogram cloud for one or all pairs.

With ivar=None (and no jvar), the cloud is computed for every pair that has a variogram model configured via set_vgm() (vgm.configured_pairs()), returning a {(ivar, jvar): cloud} dict. Otherwise a single pair is computed and returned.

Direct pairs use raw_vgm(). Cross pairs use the LMC cross-variogram estimator raw_cross_vgm() when the observations are collocated. Set cross="pseudo" to force cross_vgm(), or cross="lmc" to require collocated observations.

calc_empirical(*args, **kwargs)#: Alias for calc_experimental().

calc_average(ivar=None, jvar=None, rawvgm=None, store=True, raw_kwargs=None, **kwargs)#: Average one pair, or all cached raw variograms when no pair is given.

fit_pair(ivar, jvar=None, avgvgm=None, inplace=True, **kwargs)#

Fit one variable-pair model independently.

This is convenient for direct variograms, but fit_lmc() is the safer choice for cokriging because it constrains cross-pair sills.

fit(*, method='lmc', ivar=None, jvar=None, pairs=None, **kwargs)#

Fit the system, returning a FitResult (method facade).

method="lmc" (default) fits a joint linear model of coregionalization via fit_lmc(); method="pair" fits one variable pair independently via fit_pair() and requires ivar (with optional jvar). Subclasses extend this facade – IndicatorVariogramSystem adds method="closure".

fit_aniso_angle(ivar=1, jvar=None, pairs=None, n_struct=None, set_ranges=True, raw_kwargs=None, **kwargs)#

Fit the anisotropy orientation from one pair and apply it system-wide.

Run this before fit_lmc() / fit_pair(): the orientation is estimated from the (ivar, jvar) pair’s empirical cloud (the primary auto-variogram by default) with the model-based profile fit (krigekit.fit_aniso_angle(), 3-D) or the PCA azimuth (estimate_aniso_angle(), 2-D), then written into the structures selected by pairs (all configured pairs by default, since an LMC shares one orientation). set_ranges=True also seeds the minor ranges from the fitted ratios. Returns self.

fit_lmc(pairs=None, x_col=('distance', 'mean'), y_col=('variogram', 'mean'), sigma_col=None, weight_col=None, fit_ranges=True, fit_nugget=True, inplace=False, raw_kwargs=None, avg_kwargs=None, max_nfev=20000, **kwargs)#

Fit an additive linear model of coregionalization.

The sill matrix for each nested structure is parameterized as L @ L.T during optimization, so every fitted coregionalization matrix is positive semidefinite. Ranges are shared across all pairs. Use sigma_col for uncertainty-style residual scaling or weight_col for weighted least squares.

get_lmc_matrices(include_nugget=True)#: Return fitted/coregionalization matrices from current pair models.

set_markov_cross(primary, secondary, corr=None, structure='secondary', cross_nugget=0.0)#

Build the (primary, secondary) cross variogram by the Markov Model 1 (MM1) collocated-cokriging assumption.

For a sparsely sampled primary and a densely sampled secondary the cross-covariance cannot be fit from the primary’s (often nugget-dominated) structure – a joint fit_lmc() would drive it to zero. MM1 instead transfers it through the collocated correlation: the cross adopts the nested structure of one variable (structure="secondary" by default – the dense covariate that carries the spatial continuity), and each cross partial sill is

\[b_{ps}^{(k)} = \rho \, \sqrt{b_{pp}^{(k)} \, b_{ss}^{(k)}}\]

which is positive-semidefinite per structure for |rho| <= 1, so the coregionalization is valid by construction (no clamping needed). This is the appropriate model for sparse-hard + dense-soft cokriging (Almeida & Journel, 1994; Goovaerts, 1997). Markov Model 2 is not yet implemented.

Parameters:

primary (int) – 1-based indices; both auto-models must already be set via set_vgm() and share the same nested-structure count.
secondary (int) – 1-based indices; both auto-models must already be set via set_vgm() and share the same nested-structure count.
corr (float, optional) – Collocated cross-correlation in [-1, 1]. If None it is estimated from the collocated observations of the two variables (which must share coordinates; otherwise pass corr explicitly).
structure ({"secondary", "primary"}) – Which variable’s structure shapes/ranges the cross adopts.
cross_nugget (float) – Cross nugget partial sill (default 0).

validate_for(kriging, observations=True, variograms=True, pairs=None)#: Check that this system can be applied to kriging without leaving it partially configured. Raises on the first problem and mutates nothing; this is the validate-first half of apply().

apply_observations(kriging)#: Transfer configured observations (and drift) in ascending order.

apply_variograms(kriging, replace=True, pairs=None)#: Transfer configured pair structures in canonical upper-triangle order.

apply(kriging, observations=True, variograms=True, replace=True, pairs=None)#

Validate, then transfer observations and/or variograms to kriging.

Validation (validate_for()) completes before any mutation, so a failure leaves kriging unchanged. Observations are applied in ascending variable order (drift immediately after each base observation), then configured variograms in canonical upper-triangle order.

class krigekit.IndicatorVariogramSystem(categories=None, ncat=None)#

Bases: krigekit.variogram_system.VariogramSystem

Variogram system for K mutually exclusive categorical indicators.

Create an indicator system from category labels or a count.

Provide categories (an ordered label list) or ncat. The system has exactly K variables (ivar = 1..K); secondary co-variates are out of scope here – configure those on the kriging object directly.

property ncat#: Number of categories K (always the system’s variable count).

encode_indicators(categories)#: Return the one-hot indicator matrix for categories using labels.

set_categorical_obs(coord, categories, *, category_labels=None, variance=None)#

Encode raw categories into K indicator datasets and store them.

category_labels (or the constructor categories) fixes the label order; otherwise sorted unique values are used. Proportions are computed from the encoded indicators.

calc_proportions()#: Compute and store category proportions from configured indicators.

initial_indicator_covariance()#: Return the closed zero-lag covariance diag(p) - p p^T.

set_indicator_vgm(vtype='sph', nugget=0.0, sill=1.0, a_major=1.0, a_minor1=None, a_minor2=None, azimuth=0.0, dip=0.0, plunge=0.0, sill_strategy='theoretical', cross_strategy='closure', proportions=None)#

Configure all K^2 indicator pairs with one shared model shape.

sill_strategy sets the auto (diagonal) partial sills:

"theoretical" – p_k (1 - p_k) from category proportions;
"uniform" – the given sill for every category.

cross_strategy sets the off-diagonal partial sills:

"closure" (recommended) – -p_k p_l, the closed indicator covariance off-diagonal;
"independent" – zero (no co-kriging);
"proportional" – sqrt(s_k s_l) from the auto sills;
"uniform" – the given sill.

Legacy configurations map directly: old cross="same" is sill_strategy="uniform", cross_strategy="uniform"; old cross="independent"/"proportional" pair with sill_strategy="theoretical".

fit_indicator_lmc(pairs=None, shared_ranges=True, x_col=('distance', 'mean'), y_col=('variogram', 'mean'), sigma_col=None, weight_col=None, fit_nugget=True, inplace=False, raw_kwargs=None, avg_kwargs=None, max_nfev=20000, **kwargs)#

Fit a closed, PSD linear model of coregionalization for indicators.

Each nested coregionalization matrix is parameterized as B = Q L L^T Q^T with Q the contrast basis, guaranteeing positive semidefiniteness and closure (B 1 = 0) at every fitted lag. Ranges are shared across pairs when shared_ranges=True. Returns a FitResult whose target is the fitted indicator system.

validate_closure(tol=1e-06, include_nugget=True)#

Check that every fitted coregionalization matrix is closed.

Raises if any matrix (sill structures, plus the nugget when include_nugget) has a row sum exceeding tol in magnitude.

fit(*, method='closure', **kwargs)#

Fit the indicator system; method="closure" is the default.

"closure" runs fit_indicator_lmc(); "lmc" and "pair" fall through to the unconstrained VariogramSystem facade.

apply(kriging, observations=True, variograms=True, replace=True, pairs=None)#: Transfer indicators and structures, checking ncat agreement.

krigekit.calc_anisotropic_lag(lag, *, anis1=1.0, anis2=1.0, azimuth=0.0, dip=0.0, plunge=0.0)#

Calculate engine-compatible equivalent major-axis lag distance.

anis1 and anis2 are the first and second minor/major range ratios. Isotropic input therefore uses both ratios equal to one. The returned lag remains in coordinate units; dividing it by the major range gives the reduced lag used by the variogram kernels.

krigekit.calc_lag_vectors(coord0, coord1, pairwise=False)#

Return signed lag vectors coord1 - coord0.

Coordinates may be single points or (n, dim) arrays. By default, matching rows are compared and a single point is broadcast. With pairwise=True the result has shape (n0, n1, dim).

krigekit#

Public API#

Submodules#

Classes#

Functions#

Package Contents#

Array convention#

Typical workflow#

Fields#

Typical workflow#

Coordinate convention#

Typical workflow (single variable, sum-metric)#