Collocated Cokriging

[n0228a]

This PR adds collocated cokriging methods to GSTools for multivariate geostatistical estimation, allowing users to improve sparse primary variable estimates by leveraging densely-sampled secondary variable data.

Collocated Cokriging Features

• Two cokriging variants: Simple Collocated (SCCK) and Intrinsic Collocated (ICCK)
• Inherits from Krige base class for full API compatibility
• Implements Markov Model I (MM1): C_YZ(h) = ρ_YZ(0)·√(C_Z(h)·C_Y(h))
• Accessible via gstools.cokriging submodule or top-level: from gstools import SimpleCollocated, IntrinsicCollocated
• Full variance estimation with return_var=True
• Parameter validation: cross_corr ∈ [-1, 1], secondary_var > 0
• Boundary conditions: ρ=0 recovers Simple Kriging, ρ=±1 gives zero ICCK variance
• SCCK can exhibit variance inflation (σ²_SCCK ≥ σ²_SK) in unstable configurations
• ICCK provides stable variance: σ²_ICCK = (1-ρ₀²)·σ²_SK ≤ σ²_SK

Simple Collocated Cokriging

Uses only collocated secondary data at the estimation point:

import gstools as gs
model = gs.Gaussian(dim=1, var=0.5, len_scale=2)
scck = gs.SimpleCollocated(
    model, cond_pos, cond_val,
    cross_corr=0.8, secondary_var=1.5,
    mean=1.0, secondary_mean=0.5
)
field, var = scck(pos, secondary_data=sec_data, return_var=True)

Intrinsic Collocated Cokriging

Uses collocated secondary data plus secondary values at all primary locations for more stable variance:

icck = gs.IntrinsicCollocated(
    model, cond_pos, cond_val,
    sec_cond_pos, sec_cond_val,
    cross_corr=0.8, secondary_var=1.5,
    mean=1.0, secondary_mean=0.5
)
field, var = icck(pos, secondary_data=sec_data, return_var=True)

Implementation Details

• Located in new gstools.cokriging submodule with base class and methods
• CollocatedCokriging base class handles MM1 formulation and variance computation
• SimpleCollocated implements SCCK estimator: Z*_SCCK = Z*_SK·(1-k·λ_Y0) + λ_Y0·(Y(u0)-m_Y) + k·λ_Y0·m_Z
• IntrinsicCollocated implements ICCK estimator using secondary weights
• Both methods require secondary_data parameter on call
• Cross-covariance computed as: C_YZ(0) = ρ_YZ(0)·√(C_Z(0)·C_Y(0))
• Inherits all Krige functionality: grids, models, dimensions, anisotropy, drift, normalizers

Tests

11 tests in test_cokriging.py testing only cokriging logic:

• Parameter validation (secondary_data required, cross_corr bounds, secondary_var positive, ICCK data length)
• Boundary conditions (ρ=0 equals SK, ρ=±1 near-zero ICCK variance)
• Variance formulas (SCCK and ICCK mathematical correctness)
• Edge cases (SCCK variance inflation, SCCK vs ICCK stability comparison)
• Exact Interpolation at conditioning points

All tests passing with no warnings.

Examples

Two minimal examples following GSTools conventions (76 & 78 lines):

• 10_simple_collocated_cokriging.py - SCCK demonstration
• 11_intrinsic_collocated_cokriging.py - ICCK demonstration
• Same data with 2 spatial gap features demonstrating cokriging value
• Two-panel plots showing raw data (left) and SK vs cokriging comparison (right)
• Synthetic correlated secondary data for reproducibility

References

• Samson, M., & Deutsch, C. V. (2020). Collocated Cokriging. Geostatistics Lessons. http://geostatisticslessons.com/lessons/collocatedcokriging
• Wackernagel, H. (2003). Multivariate Geostatistics. Springer, Berlin.

[MuellerSeb]

This is great! I love it. I wanted this for a long time but never found the time to work on it.

I have some general questions/remarks:

1. What about creating a Correlogram (abstract)-base-class with the MarkovModel1 as one child that holds all the logic currently added with several keywords in CollocatedCokriging.
   I think it could look like this:

import gstools as gs

correlo_model = gs.MarkovModel1(
    model=gs.Gaussian(dim=1, var=0.5, len_scale=2),
    cross_corr=0.8,
    secondary_var=1.5,
)
scck = gs.SimpleCollocated(
    correlo_model, cond_pos, cond_val,
    mean=1.0, secondary_mean=0.5
)

Maybe the means could be combined in a list then, or we could also move them into the MarkovModel1 (not sure what makes more sense here), or keep them as is.
We would then also need to think about which functionality we move into the Correlogram class. (e.g. _compute_covariances)
As always my motivation for this would be, to be future prove (for MM2 and maybe other possible models).

2. How does/should Cokriging work with normalizer, fit_normalizer and fit_variogram?

3. Do we also need an estimator for the empirical correlogram?

4. please use math environments in the documentation for formuals.

5. There should be a _apply_intrinsic_collocated method like the _apply_simple_collocated method.

[n0228a]

1: I'm not sure if this will work very well with the approach I chose. As you maybe have seen, after some linear algebra reformulations of the problem, SCCK and ICCK can be shown as an extension to simple Kriging, making it not necessary to construct new matrices, that would need to be solved for every estimation point in collocated cokriging. That is also where the normalizer goes currently., it applies in the way it would apply for Simple Kriging.
3. As for the Markov Model 1 so far, I don't think a cmpirical correlogramm would be necessary, for mm2 also not as I understand the theory (there you would need to fit to the primary variogram again based on your secondary variable and the correlation coefficient).
4. I will update that.
5. I will update that too.

[LSchueler]

Hi, I also really like your work, thanks a lot!

I suggest that we have a chat about @MuellerSeb's suggestions. The Correlogram base class he is talking about would ensure that we can later add more cokriging methods without breaking the interface.

The source checks fail. You have to run a linter on your code, see the instructions at the end of this file.

I also think that cokriging is missing after line 148 of src/gstools/__init__.py.

Depending on the time you still have for this work, I think an example with 2d data would be really valuable, as I guess most applications will be 2d and many of our users like to have very specific examples.

[n0228a]

Hi,
I updated the documentation and added the apply_intrinsic method. The linter wants to sort by alphabet, but that raises import errors, the only fix I found is to import Cokrige after Krige, because it inherits from Krige. Maybe you know something.

[LSchueler]

I think there are no import problems, the linter is just too stupid.
You can change line 148 in src/gstools/__init__.py to

from gstools import (  # noqa: I001

put cokriging right after kriging and the linter should stop complaining.

[n0228a]

I tried to trace the circular import error and it occurs because krige needs field, it then initializes field with the submodule cond_srf and this needs again krige. I still don't understand why this doesn't happened before.
krige → field → cond_srf → krige
I'm not sure what the preferred solution here would be or just keep as is and put cokrige in different position?