Tutorials
Juyter Logo

Directional variogram estimation and fitting in 2D

%matplotlib widget
import matplotlib.pyplot as plt
plt.ioff()
# turn of warnings
import warnings
warnings.filterwarnings('ignore')

In this example, we demonstrate how to estimate a directional variogram by setting the direction angles in 2D.

Afterwards we will fit a model to this estimated variogram and show the result.

import numpy as np
import gstools as gs

Generating synthetic field with anisotropy and a rotation of 22.5 degree.

angle = np.pi / 8
model = gs.Exponential(dim=2, len_scale=[10, 5], angles=angle)
x = y = range(101)
srf = gs.SRF(model, seed=20220425)
field = srf((x, y), mesh_type="structured")

Now we are going to estimate a directional variogram with an angular tolerance of 11.25 degree and a bandwith of 8.

bins = range(0, 40, 2)
bin_center, dir_vario, counts = gs.vario_estimate(
    pos=(x, y),
    field=field,
    bin_edges=bins,
    direction=gs.rotated_main_axes(dim=2, angles=angle),
    angles_tol=np.pi / 16,
    bandwidth=8,
    mesh_type="structured",
    return_counts=True,
)

Afterwards we can use the estimated variogram to fit a model to it:

print("Original:")
print(model)
model.fit_variogram(bin_center, dir_vario)
print("Fitted:")
print(model)
Original:
Exponential(dim=2, var=1.0, len_scale=10.0, nugget=0.0, anis=[0.5], angles=[0.393])
Fitted:
Exponential(dim=2, var=0.921, len_scale=10.7, nugget=2.85e-24, anis=[0.541], angles=[0.393])

Plotting.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5])

ax1.scatter(bin_center, dir_vario[0], label="emp. vario: 1/8 $\pi$")
ax1.scatter(bin_center, dir_vario[1], label="emp. vario: 5/8 $\pi$")
ax1.legend(loc="lower right")

model.plot("vario_axis", axis=0, ax=ax1, x_max=max(bin_center), label="fit on axis 0")
model.plot("vario_axis", axis=1, ax=ax1, x_max=max(bin_center), label="fit on axis 1")
ax1.set_title("Fitting an anisotropic model")

srf.plot(ax=ax2)
ax2.set_aspect("equal")

Without fitting a model, we see that the correlation length in the main direction is greater than the transversal one.