[SOLVED] Include the Matheron estimator in pykrige

Include the Matheron estimator in pykrige

Hi I have another question: how can I include the Matheron estimator in my program and still use Pykrige?

My program looks like this:

from pykrige.ok import OrdinaryKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt

x = np.array([-100, 280, -290, 23, 101, 110])
y = np.array([54, 100, 590, 470, 200, 25])
z = np.array([29.3, 21.0, 19.2, 29.1, 21.9, 23.1])

Ok = OrdinaryKriging(
    x,
    y,
    z,
    variogram_model='exponential ',
    verbose= True,
    enable_plotting= True
)

gridx =np.arange(-300, 300, 10, dtype='float64')
gridy =np.arange(0,600, 10, dtype='float64')
zstar, ss = Ok.execute("grid", gridx, gridy)

cax = plt.imshow(zstar, extent=(-300, 300, 0, 600), origin='lower')
plt.scatter(x, y, c=z, marker=".")
cbar = plt.colorbar(cax)
plt.title('Test Interpolation')
plt.show()

Solution

Sebastian here from the GeoStat-Framework. The Matheron estimator for the empirical variogram is used by default in PyKrige (and is the only option).

Since PyKrige and GSTools are working well together, you could use GSTools for the variogram analysis part, where you can also choose the "cressie" estimator (or "matheron" by default as well):

from matplotlib import pyplot as plt
from pykrige.ok import OrdinaryKriging
import gstools as gs
import numpy as np

x = np.array([-100, 280, -290, 23, 101, 110])
y = np.array([54, 100, 590, 470, 200, 25])
z = np.array([29.3, 21.0, 19.2, 29.1, 21.9, 23.1])

model = gs.Exponential(dim=2)
bin_center, gamma = gs.vario_estimate(
    pos=[x, y],
    field=z,
    bin_edges=range(0, 900, 150),
    estimator="cressie",
)
model.fit_variogram(bin_center, gamma)
ax = model.plot(x_max=max(bin_center))
ax.scatter(bin_center, gamma)

Ok = OrdinaryKriging(
    x,
    y,
    z,
    variogram_model=model,
    verbose=True,
    enable_plotting=True
)

gridx =np.arange(-300, 300, 10, dtype='float64')
gridy =np.arange(0,600, 10, dtype='float64')
zstar, ss = Ok.execute("grid", gridx, gridy)

cax = plt.imshow(zstar, extent=(-300, 300, 0, 600), origin='lower')
plt.scatter(x, y, c=z, marker=".")
cbar = plt.colorbar(cax)
plt.title('Test Interpolation')
plt.show()

The estimation in PyKrige has some slightly different assumptions (like a set number of lags), but in GSTools you have more fine control.

Hope this helps.

Cheers, Sebastian