Include the Matheron estimator in pykrige

Question

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()

score 1 · Accepted Answer · answered Jun 13 '22 at 10:36

1

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

answered Jun 13 '22 at 10:36

MuellerSeb

796
6
11

Hi Sebastian, thanks for your answer I have another question about Pykrige and hope you can answer it for me: https://stackoverflow.com/questions/72574394/advance-a-basemap-plot – Weiss Jun 13 '22 at 10:54
hi I want to plot Isotopic data which are dependend from height. So I would use more likly universal kriging. However I wounder If I schould use for this propose 3D UK or normal UK? and which drift function would you recomend for Interpolating precipitation data? – Weiss Jun 14 '22 at 09:27
3D UK is for 3D data (e.g. depending on lat-lon-altitude). But I guess precipitation is mostly only given for lat-lon coordinates. We often use the DEM as an external drift (called specified drift in PyKrige and needs to be provided at the output grid as well). – MuellerSeb Jun 14 '22 at 13:55
Thanks for your answer! However in fact, the isotope data depends on the altitude, so 3DUK seems to be a good choice, doesn't it? Or should I just take the UK and use the altitude as a driver? – Weiss Jun 14 '22 at 14:01
You mean as a "drift"? You could try that ;-) – MuellerSeb Jun 14 '22 at 14:03
Since I have not been programming for very long, I would like to know if there is an example somewhere that explains how the drift is composed! and yes I mean drift ;) – Weiss Jun 14 '22 at 14:09
Here is an example with specified drift (3 drifts, you only want 1) in 3D: https://geostat-framework.readthedocs.io/projects/pykrige/en/stable/examples/02_kriging3D.html#sphx-glr-examples-02-kriging3d-py – MuellerSeb Jun 15 '22 at 08:51
That means: k3d3, ss3d = uk3d.execute( "grid", gridx, gridy, gridz, specified_drift_arrays=[ zg] or? – Weiss Jun 15 '22 at 08:55
maby can this be right? https://stackoverflow.com/questions/72627986/plotting-3d-data-in-2d-map (but all in all thanks for your patience:)) – Weiss Jun 15 '22 at 09:37
Hi @greeeeeeen do you know if there have ever been errors when trying to display UK3D in Basemap as a 2D image? – Weiss Jun 15 '22 at 10:54
hi, sry that I ask you so often, hovever I think you are the only one who can help me in this question: https://stackoverflow.com/questions/72631305/advance-a-interpolation I would be very happy if you could have a look at this :) All the best – Weiss Jun 18 '22 at 10:48

Include the Matheron estimator in pykrige

1 Answers1