%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 gsGenerating 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")Loading...
Without fitting a model, we see that the correlation length in the main direction is greater than the transversal one.