pyinterp.OptimalInterpolation

class pyinterp.OptimalInterpolation(obs_coords, obs_values, obs_sigma2, *, covariance='gaussian', coordinate_system='cartesian', spheroid=None, time_scale=1.0)

Bases: object

Optimal Interpolation with anisotropic covariance and per-obs error.

The estimator indexes scattered observations and computes the BLUE (Best Linear Unbiased Estimator) at arbitrary query points.

Important

This is simple kriging: the BLUE about a known, zero mean. The weights are not constrained to sum to one, so observations must be anomalies relative to a zero background mean — subtract a mean / reference field before passing them in (e.g. use sea-level anomalies, not absolute sea-surface height). Where a query has no nearby, well-correlated observations the analysis relaxes toward 0 (the prior mean) and the formal error toward sigma, not toward the local data mean. Feeding non-centred values therefore produces a systematic bias toward zero in data-sparse regions.

The field standard deviation sigma is also evaluated only at the query point and assumed locally constant across that query’s neighbourhood; strong sigma gradients within one neighbourhood introduce a small modelling bias.

Parameters:

• obs_coords (NDArray2DFloat64) – Observation coordinates of shape (N, 3). In cartesian mode this is (x, y, t) in user units. In geographic mode this is (lon_deg, lat_deg, t_seconds).

• obs_values (NDArray1DFloat64) – Observed values, shape (N,).

• obs_sigma2 (NDArray1DFloat64) – Per-observation error variance, shape (N,). Must be strictly positive.

• covariance (CovarianceFunction) – Anisotropic covariance kernel name. Defaults to "gaussian". Note that "spherical" and "wendland" are positive-definite only in up to three dimensions; in geographic mode the analysis runs over three ECEF axes plus time, so for large decorrelation scales the covariance matrix is not guaranteed positive-definite with those two kernels — prefer "gaussian", "matern_12/32/52" or "cauchy" there.

• coordinate_system (CoordinateSystem) – "cartesian" (default) or "geographic". See module docstring for the trade-offs.

• spheroid (Spheroid | None) – Optional pyinterp.geometry.geographic.Spheroid used in geographic mode. Defaults to WGS-84.

• time_scale (float) – Scale factor applied to t when packing the internal R-tree, so the k-NN search (and any radius filter) treat space and time comparably. In cartesian mode it is in user-space-units per time-unit; in geographic mode it is in meters per second (time becomes a meters-equivalent before being mixed with the ECEF axes). It cancels out of the covariance — only neighbour selection is affected, so results with the default 1.0 are unchanged. Default 1.0.

Note

The internal R-tree uses the 4D Euclidean metric on the packed coordinates, mixing the three spatial axes with the time_scale-scaled time axis. In geographic mode the spatial part is the geodetic chord distance (very close to the great-circle distance for short ranges); set time_scale to a meters-per-second value so the k-NN search balances space and time. In cartesian mode the user is responsible for unit consistency between space and time.

Choosing time_scale:

The k-NN search ranks neighbours by the raw packed Euclidean distance, which mixes spatial units with time. With the default time_scale = 1 that ranking depends on the arbitrary ratio of spatial units to seconds, not on how correlated the neighbours actually are, so the retrieved set can be skewed toward one axis. A neighbour one decorrelation length away in space and one away in time are equally correlated with the query; to make the search treat them equally, scale time so those two steps have the same packed length:

time_scale ≈ L_spatial / L_t

(typical spatial decorrelation length over the temporal one). Then the k nearest neighbours are the k most correlated observations. Worked example (geographic): L_spatial = 150 km and L_t = 6 h give time_scale ≈ 150000 m / 21600 s ≈ 7 m/s. time_scale is baked into the tree at construction, so choose a value representative of the scales you query with; it never changes the covariance (hence the analysis is identical when all neighbours are used), only which neighbours are retrieved.

Index the observations and build the internal 4D R-tree.

Attributes

coordinate_system	Coordinate system in use.
covariance	Covariance kernel name.
n_observations	Number of indexed observations.

Special Methods

__call__(query_coords, *[, lx, ly, ...])

Run the OI analysis at a set of query points.