pykriging#

Python wrapper for a Fortran kriging and SGSIM engine.

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

Public API#

Entry point

Kriging(st=False, **kwargs)

Factory that returns a _SpatialKriging when st=False (the default) or a SpaceTimeKriging when st=True.

Spatial usage:

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

Space-time usage:

k = Kriging(st=True, nvar=1)
k.set_st_model(...)

Classes (returned by Kriging factory)

_SpatialKriging — spatial kriging / co-kriging / SGSIM SpaceTimeKriging — 3-D + time kriging / SGSIM

Convenience functions

ordinary_kriging — one-shot point kriging cokriging — one-shot co-kriging sequential_gaussian_simulation — one-shot SGSIM spacetime_kriging — one-shot ST kriging spacetime_cokriging — one-shot ST co-kriging

Submodules#

Classes#

SpaceTimeKriging

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

Functions#

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

Package Contents#

pykriging.Kriging(st: bool = False, **kwargs)#

Create a kriging object — spatial or space-time.

This is the primary entry point for both backends. Pass st=True to obtain a SpaceTimeKriging instance; omit it (or st=False) for ordinary spatial kriging (_SpatialKriging).

Parameters:

st (bool, default False) – False — spatial kriging (_SpatialKriging). True — space-time kriging (SpaceTimeKriging).
**kwargs –
Forwarded verbatim to the chosen constructor.

Spatial-only parameters (ndim, varying_vgm, std_ck, pf_cache) are silently removed when st=True; a UserWarning is emitted for any that were supplied with a non-default value so they are not silently ignored.

Returns:

_SpatialKriging – When st=False. Full set of methods: set_obs(), set_vgm(), set_grid(), set_search(), solve(), get_results(), and more.
SpaceTimeKriging – When st=True. Same core methods plus set_st_model(), set_vgm_temporal(), and set_vgm_joint_sills(). set_grid() requires an extra time= argument. set_search() requires extra time_* arguments.

Notes

Factory, not a class: Kriging is a function that returns one of two concrete types. Use isinstance(k, _SpatialKriging) or isinstance(k, SpaceTimeKriging) for type checks.

Output shape: get_results() shapes are consistent within each backend but differ when nsim > 1:

spatial — estimate (nblocks, nvar, nsim) squeezed to (nblocks,) for nvar=1, nsim=1.
ST — estimate (nsim, nblocks) squeezed to (nblocks,) for nsim=1.

Examples

Spatial ordinary kriging:

k = Kriging(ndim=2, nvar=1)
k.set_obs(ivar=1, coord=coord, value=value, nmax=20)
k.set_grid(coord=grid_coord)
k.set_vgm(ivar=1, jvar=1, vtype="sph", sill=1.0, a_major=500)
k.set_search(ivar=1)
k.solve()
est, var = k.get_results()

Space-time ordinary kriging:

k = Kriging(st=True, nvar=1)
k.set_st_model(model='sum_metric', transform='bounded', at=5.0)
k.set_obs(ivar=1, coord=obs_coord4, value=obs_value, nmax=30)
k.set_vgm(ivar=1, jvar=1, vtype="sph", sill=0.8, a_major=1000)
k.set_vgm_temporal(ivar=1, jvar=1, vtype="exp", sill=0.6, at_k=10.0)
k.set_vgm_joint_sills(1, 1, 0.4)
k.set_grid(coord=grid_coord, time=grid_time)
k.set_search(ivar=1)
k.solve()
est, var = k.get_results()

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

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

Parameters:

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

Returns:

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

Example

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

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

One-shot ordinary co-kriging with multiple variables.

Parameters:

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

Returns:

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

Example

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

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

Sequential Gaussian Simulation.

Parameters:

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

Returns:

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

class pykriging.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, 4) shape — columns 0:3 are x,y,z and column 3 is time. Grid/block arrays still use (ngrid, 3) + time separately until set_grid is updated.

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_coord4, 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

nvar = 1#

ndrift = 0#

nsim = 0#

verbose = False#

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 ('sum_metric' or 'product_sum')
transform ('nug' | 'sph' | 'exp' | 'gau' | 'pow' | 'bsq' | 'cir' | 'lin') – Controls f(dt) used in the joint ST distance. Aliases: ‘linear’ -> ‘lin’, ‘bounded’ -> ‘exp’, ‘power’ -> ‘pow’.
at (joint temporal scale (same time units as observations))
time_nugget –

: f(dt) = nugget + sill * (1 - corefunc(|dt| / at)) for
nonzero dt; f(0) is always 0.
time_sill –

: f(dt) = nugget + sill * (1 - corefunc(|dt| / at)) for
nonzero dt; f(0) is always 0.
k_ps (product-sum coefficient k (model='product_sum' only))

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

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

Parameters:

ivar (variable index, 1-based)
coord ((nobs, 4) — columns 0:3 spatial (x,y,z), column 3 = time)
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, 4) with columns x, y, z, t. 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, rangescale: numpy.ndarray | None = None, localnugget: numpy.ndarray | None = None)#

Set point estimation targets.

coord : (ngrid, 3) spatial coordinates time : (ngrid,) prediction times

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(). If randpath/sample are None, they are generated internally.

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

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

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

Parameters:

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

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

Run the ST kriging or SGSIM loop.

Parameters:

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

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

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 (nsim, ngrid)
variance (ndarray, shape (ngrid,))

pykriging.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, 4) observation coordinates — columns 0:3 spatial, column 3 time)
obs_value ((nobs,) observed values)
grid_coord ((ngrid, 3) 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,))

pykriging.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, 4) arrays, one per variable — columns 0:3 spatial, column 3 time)
obs_values (list of (nobs_i,) arrays)
grid_coord ((ngrid, 3))
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,))