krigekit.kriging_st#

Python wrapper for the space-time kriging C API (krige_st_* entry points).

Mirrors the structure of _kriging.py but exposes:

  • SpaceTimeKriging class — full control over the ST kriging workflow

  • spacetime_kriging() — one-shot ordinary ST kriging

  • spacetime_cokriging() — one-shot ordinary ST co-kriging

Coordinate convention (same as base):

All spatial coord arrays are (nobs, 3) — rows are points, columns are x, y, z. Time arrays are 1-D, shape (nobs,), in any consistent unit (e.g. decimal years).

Variogram spec formats:

Spatial (9 values): “vtype nugget sill a_major a_minor1 a_minor2 azimuth dip plunge” Temporal (4 values): “vtype nugget sill at_k”

ST model parameters (set once via set_st_model):

model : ‘sum_metric’ or ‘product_sum’ transform : ‘nug’, ‘sph’, ‘exp’, ‘gau’, ‘pow’, ‘bsq’, ‘cir’, or ‘lin’ at : joint temporal scale (same time units as input) time_nugget, time_sill

: f(dt) scale for the joint temporal distance

k_ps : product-sum coefficient k (only for model=’product_sum’)

Classes#

SpaceTimeKriging

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

Functions#

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.

Module Contents#

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

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

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

Coordinate convention#

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

Typical workflow (single variable, sum-metric)#

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

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

Parameters:

  • model (str) – 'sum_metric' or 'product_sum'.

  • transform (str) – Variogram type used for f(dt): 'nug' | 'sph' | 'exp' | 'gau' | 'pow' | 'bsq' | 'cir' | 'lin'. Aliases: 'linear' → 'lin', 'bounded' → 'exp', 'power' → 'pow'.

  • at (float) – Joint temporal scale (same time units as observations).

  • time_nugget (float) – Nugget jump in f(dt) for the sum-metric model (applied for dt ≠ 0).

  • time_sill (float) – Upper scale in f(dt): f(dt) = time_nugget + time_sill * (1 - corefunc(|dt| / at)) for dt ≠ 0; f(0) = 0 always.

  • k_ps (float) – Product-sum coefficient k (model='product_sum' only).

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

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

Accepts two coordinate formats — pick whichever matches your workflow:

Combined format (default):

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

Split format (explicit time array):

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

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

Parameters:

  • value ((nobs,) observed values)

  • variance ((nobs,) measurement error variance (default: zeros))

  • nmax (max neighbours)

  • maxdist (max search radius in km-equivalent space (same units as h_ST))

  • sk_mean (global mean for simple kriging (unbias=0); default 0)

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

Replace observation values for variable ivar in-place.

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

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

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

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

Set time-aware ST gradient observation pairs.

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

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

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

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

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

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

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

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

Set joint sills for the sum-metric model.

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

Example:

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

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

Set point estimation targets.

Accepts two coordinate formats — pick whichever matches your workflow:

Combined format (consistent with set_obs()):

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

Split format (explicit time array):

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

  • rangescale ((ngrid,) local range scale factors (default: ones))

  • localnugget ((ngrid,) local nugget additions (default: zeros))

set_grid_cv()#

Cross-validation mode: predict at observation locations.

set_grid_drift(drift: numpy.ndarray)#

Drift values at estimation grid. drift shape: (ngrid, ndrift).

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

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

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

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

Parameters:

  • ivar (int) – Variable index (1-based).

  • time_at (float) – Temporal scale factor (default 1.0) to convert the time axis into km-equivalent search units: the search-tree time coordinate is t * time_at. This ensures that L2 distance in the 4D search space matches the sum-metric space-time distance: h_ST = sqrt(h_S^2 + (time_at * dt)^2). Normally, you should pass the same value as at in set_st_model().

  • anis1 (float) – Spatial minor/major anisotropy ratio (default 1.0).

  • anis2 (float) – Spatial vertical/major anisotropy ratio (default 1.0).

  • azimuth (float) – Azimuth of the spatial major axis in degrees (default 0.0, clockwise from North).

  • dip (float) – Dip angle of the spatial major axis in degrees (default 0.0, positive downward).

  • plunge (float) – Plunge angle of the spatial major axis in degrees (default 0.0).

  • sector_search (bool) – Enable sector (octant) search limiting candidates per sector. If True, candidate neighbours are partitioned into 8 spatial octants centered on the prediction location. At most nmax (from set_obs()) candidates are selected per octant. This ensures a balanced spatial distribution of neighbours and prevents clustering artifacts. The maximum total neighbours selected is 8 * nmax. If search anisotropy is enabled, spatial coordinates are rotated/scaled according to the anisotropy parameters before sector assignment.

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

Run the ST kriging or SGSIM loop.

Parameters:

  • nthread (int, optional) – Maximum number of OpenMP threads. 0 (default) lets the OpenMP runtime choose (respects OMP_NUM_THREADS).

  • ncache (int, optional) – Number of per-thread multi-slot hcache entries for this solve. None keeps the compiled/object default, 0 disables hcache, and 1 builds a one-slot hcache for overhead comparisons.

get_factor() → dict#

Return cached factor matrices and the assembled linear system.

property solver_stats: dict#

Solver statistics from the most recent solve() call.

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

chol_ok

Blocks solved by Cholesky factorization (either a fresh factorize or a cache hit that reused a previously computed Cholesky factor).

ssytrf_fact

Number of SSYTRF (Bunch-Kaufman LDL^T) factorizations performed. Each one is O(n³) but occurs only once per unique neighbourhood. A non-zero value means Cholesky failed for at least one neighbourhood; a value equal to 1 with global search means the factorization was done once and cached for all blocks.

ssytrf_reuse

Blocks solved by a cached SSYTRF factorization using SSYTRS, which is O(n²). When this is large relative to ssytrf_fact the SSYTRF caching is working effectively.

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

k.solve()
s = k.solver_stats
# Expected: chol_ok=0, ssytrf_fact=1, ssytrf_reuse=9999
get_results(copy: bool = False, squeeze: bool = True) → tuple[np.ndarray, np.ndarray]#

Retrieve kriging estimate and variance.

Parameters:

  • copy (bool, default False) – If True, return C-contiguous copies for downstream NumPy/Pandas use. If False, return views / Fortran-order arrays when possible.

  • squeeze (bool, default True) – If True, return a 1-D estimate when nsim == 1.

Returns:

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

  • variance (ndarray, shape (ngrid,))

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

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

Parameters:

  • obs_coord ((nobs, ndim+1) observation coordinates — first ndim cols spatial, last col time)

  • obs_value ((nobs,) observed values)

  • grid_coord ((ngrid, ndim) prediction spatial coordinates)

  • grid_time ((ngrid,) prediction times)

  • spatial_spec (dict or list[dict] spatial variogram structure(s))

  • temporal_spec (dict or list[dict] temporal variogram structure(s))

  • joint_sills (list[float] joint sills (sum-metric only))

  • model ('sum_metric' or 'product_sum')

  • transform ('nug', 'sph', 'exp', 'gau', 'pow', 'bsq', 'cir', or 'lin')

  • at (joint temporal scale (also used as time_at for the KD-tree))

  • time_nugget (temporal variogram nugget/sill for set_vgm_temporal)

  • time_sill (temporal variogram nugget/sill for set_vgm_temporal)

  • nmax (max neighbours)

  • maxdist (max search radius in km-equivalent space (h_ST units))

  • nthread (max OMP threads for this call (0 = OMP default))

  • ncache (per-thread hcache slots for this solve; None uses default)

Returns:

  • estimate ((ngrid,))

  • variance ((ngrid,))

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

One-shot ordinary space-time co-kriging.

Parameters:

  • obs_coords (list of (nobs_i, ndim+1) arrays, one per variable — first ndim cols spatial, last col time)

  • obs_values (list of (nobs_i,) arrays)

  • grid_coord ((ngrid, ndim))

  • grid_time ((ngrid,))

  • spatial_specs (dict (ivar,jvar) -> dict or list[dict])

  • temporal_specs (dict (ivar,jvar) -> dict or list[dict])

  • joint_sills (dict (ivar,jvar) -> list[float])

  • nthread (max OMP threads for this call (0 = OMP default))

  • ncache (per-thread hcache slots for this solve; None uses default)

Returns:

  • estimate ((ngrid,))

  • variance ((ngrid,))