Describe the bug
kriging() builds the ordinary-kriging matrix in _build_kriging_matrix (xrspatial/interpolate/_kriging.py) with:
K[:n, :n] = vario_func(D)
D is the pairwise distance matrix, so its diagonal is 0. Evaluating the variogram model at distance 0 returns the nugget c0 (for example _spherical(0, c0, c, a) == c0), so the nugget ends up on the matrix diagonal.
By definition the semivariogram has gamma(0) = 0: gamma(h) = 0.5 * Var(Z(x+h) - Z(x)), and that difference is identically zero at h = 0. The nugget is the one-sided limit as h approaches 0 from above, not the value at 0. Standard ordinary-kriging implementations such as PyKrige zero the diagonal of the variogram block.
Effect
When the fitted nugget is zero the diagonal is already 0 and nothing changes, which is why the current tests pass: their trend-dominated data fits a near-zero nugget. When the nugget is non-zero (real data with measurement noise or microscale variation), putting it on the diagonal:
- forces the predictor to reproduce the data exactly instead of smoothing through the nugget, which is wrong for noisy data
- biases the kriging variance downward, since the diagonal term feeds into the variance via k0^T K_inv k0
Expected behavior
The variogram block diagonal should be 0. With a non-zero nugget the predictor should smooth rather than hit data points exactly, and the kriging variance should reflect the nugget.
Fix
Zero the diagonal of the variogram block after evaluating the model in _build_kriging_matrix.
Found during an accuracy sweep of xrspatial/interpolate/_kriging.py.
Describe the bug
kriging()builds the ordinary-kriging matrix in_build_kriging_matrix(xrspatial/interpolate/_kriging.py) with:Dis the pairwise distance matrix, so its diagonal is 0. Evaluating the variogram model at distance 0 returns the nuggetc0(for example_spherical(0, c0, c, a) == c0), so the nugget ends up on the matrix diagonal.By definition the semivariogram has gamma(0) = 0: gamma(h) = 0.5 * Var(Z(x+h) - Z(x)), and that difference is identically zero at h = 0. The nugget is the one-sided limit as h approaches 0 from above, not the value at 0. Standard ordinary-kriging implementations such as PyKrige zero the diagonal of the variogram block.
Effect
When the fitted nugget is zero the diagonal is already 0 and nothing changes, which is why the current tests pass: their trend-dominated data fits a near-zero nugget. When the nugget is non-zero (real data with measurement noise or microscale variation), putting it on the diagonal:
Expected behavior
The variogram block diagonal should be 0. With a non-zero nugget the predictor should smooth rather than hit data points exactly, and the kriging variance should reflect the nugget.
Fix
Zero the diagonal of the variogram block after evaluating the model in
_build_kriging_matrix.Found during an accuracy sweep of
xrspatial/interpolate/_kriging.py.