Juyter Logo

Finding the best fitting variogram model

%matplotlib widget
import matplotlib.pyplot as plt
# turn of warnings
import warnings
import numpy as np
import gstools as gs

Generate a synthetic field with an exponential model.

x = np.random.RandomState(20220425).rand(1000) * 100.0
y = np.random.RandomState(20220426).rand(1000) * 100.0

model = gs.Exponential(dim=2, var=2, len_scale=8)
srf = gs.SRF(model, mean=0, seed=20220425)
field = srf((x, y))

Estimate the variogram of the field with 40 bins and plot the result.

bin_center, gamma = gs.vario_estimate((x, y), field)

Define a set of models to test.

models = {
    "Gaussian": gs.Gaussian,
    "Exponential": gs.Exponential,
    "Matern": gs.Matern,
    "Stable": gs.Stable,
    "Rational": gs.Rational,
    "Circular": gs.Circular,
    "Spherical": gs.Spherical,
    "SuperSpherical": gs.SuperSpherical,
    "JBessel": gs.JBessel,
scores = {}

Iterate over all models, fit their variogram and calculate the r2 score.

# plot the estimated variogram
plt.scatter(bin_center, gamma, color="k", label="data")
ax = plt.gca()

# fit all models to the estimated variogram
for name, model in models.items():
    fit_model = model(dim=2)
    para, pcov, r2 = fit_model.fit_variogram(bin_center, gamma, return_r2=True)
    fit_model.plot(x_max=max(bin_center), ax=ax)
    scores[name] = r2

Create a ranking based on the score and determine the best models

ranking = sorted(scores.items(), key=lambda item: item[1], reverse=True)
print("RANKING by Pseudo-r2 score")
for i, (model, score) in enumerate(ranking, 1):
    print(f"{i:>6}. {model:>15}: {score:.5}")
RANKING by Pseudo-r2 score
     1.          Stable: 0.99567
     2.  SuperSpherical: 0.99548
     3.          Matern: 0.9952
     4.     Exponential: 0.99251
     5.        Rational: 0.99205
     6.       Spherical: 0.98856
     7.        Circular: 0.98365
     8.        Gaussian: 0.98064
     9.         JBessel: 0.97989