Using OrdinaryKriging vs GSTools for variogram fit

[moffat]

Hi,

I'm new to GSTools/PyKrige, and I'm trying to go through some of the website's examples to understand what the code is doing. One thing I was trying to do is to compare the variogram fit using the OrdinaryKriging example here versus using the model output of fit_variogram in GSTools. Using the same "data" in the OrdinaryKriging example, I wanted to compare the output of:

import numpy as np
import gstools as gs
import pykrige.kriging_tools as kt
from pykrige.ok import OrdinaryKriging
# generate a synthetic field with an exponential model
data = np.array(
    [
        [0.3, 1.2, 0.47],
        [1.9, 0.6, 0.56],
        [1.1, 3.2, 0.74],
        [3.3, 4.4, 1.47],
        [4.7, 3.8, 1.74],
    ]
)

gridx = np.arange(0.0, 5.5, 0.5)
gridy = np.arange(0.0, 5.5, 0.5)

# estimate the variogram of the field
bin_center, gamma = gs.vario_estimate((data[:,0], data[:,1]), data[:,2])
# fit the variogram with a stable model. (no nugget fitted)
fit_model = gs.Linear(dim=2)
fit_model.fit_variogram(bin_center, gamma)
# output
ax = fit_model.plot(x_max=max(bin_center))
ax.scatter(bin_center, gamma)
print(fit_model)

OK = OrdinaryKriging(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    variogram_model="linear",
    verbose=True,
    enable_plotting=True,
)


OK_gs = OrdinaryKriging(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    fit_model,
    verbose=True,
    enable_plotting=True,
)

z, ss = OK.execute("grid", gridx, gridy)
z_gs, ss_gs = OK_gs.execute("grid", gridx, gridy)

This did not go well because the bings chosen by gs.vario_estimate were in a narrow range (max~1.75). Forcing something more reasonable like

bin_center, gamma = gs.vario_estimate((data[:,0], data[:,1]), data[:,2],bin_edges=([1.5 , 2, 3, 4, 5, 6]))

gave more reasonable results. I guess I was surprised that OrdinaryKriging was choosing what appeared to be the right bins while vario_estimate did not.

A follow-up question is that I tried the above with a Gaussian model, and could not match the results of OrdinaryKriging using vario_estimate. I'm probably missing something. Any thoughts are appreciated.

[LSchueler]

Hi moffat,

GSTools uses Sturges' rule to determine the number of bins if it's not given. This rule is known to sometimes perform poorly for small sample sizes. PyKrige on the other hand uses a fixed value of bins if none is given and I think you just got lucky with your values.
Unfortunately, the responsibilities of GSTools and PyKrige are not very clearly separated at the moment. But as long as you don't want to do very fancy kriging, you could simply stick with GSTools also for the ordinary kriging.

Does that answer your first question?

Regarding the Gaussian model, my initial guess is that the difference comes from the different bins?

[moffat]

Hi @LSchueler

Thanks, this is very helpful re: 1st question.

The Gaussian model answers are very different for each toolbox. This might the result of using these tiny test datasets. I'll try to put together some example code and figures to illustrate the differences.

Is there an obvious way to use the same bins with both toolboxes? I guess I could use the output of OrdinaryKriging and use that in vario_estimate, but is there a way to set custom bin edges with OrdinaryKriging?

[LSchueler]

I'm sorry for the late reply.

You mean you want to set the bin edges for the OrdinaryKriging from PyKrige? You can at least set the argument nlags, e.g.

nlags = 10

OK = OrdinaryKriging(
    data[:, 0],
    data[:, 1],
    data[:, 2],
    variogram_model="linear",
    nlags=nlags,
    verbose=True,
    enable_plotting=True,
)

If you really want to go down that path, you could edit this part of PyKrige to set custom bin edges.
The easier solution would be to basically copy the linked code and use that calculation to set the custom bin edges in the GSTools code.

[moffat]

I'm sorry myself for the late reply @LSchueler :)

The answer makes sense. Thank you!