Directional Variogram

[Bishama]

Hi,
I am following the tutorial for 2D directional variogram. I am generating a synthetic field using exponential model and then fitting the theoretical model manually using the same parameters used to define the field. However, the variogram doesn't fit even at shorter lag distances for 45 degree direction. If I increase the field size from 101 to 301 it gives a better fit.

Step 1: Generating synthetic field with anisotropy and a rotation of 45 degree.

angle = np.deg2rad(45) 
model = gs.Exponential(dim=2, len_scale=[25, 12.5], var=25, angles=angle)
x = y = range(101)
srf = gs.SRF(model, mean=50, seed=123456)
field = srf((x, y), mesh_type="structured")

Step2: Compute experimental variogram of the generated field

bins = range(0, 60, 2)
bin_center, dir_vario, counts = gs.vario_estimate(
    *((x, y), field, bins),
    direction=gs.rotated_main_axes(dim=2, angles=angle),
    angles_tol=np.deg2rad(15) ,
    bandwidth=8,
    mesh_type="structured",
    return_counts=True,
)

Step 3: Fit the manual theoretical variogram to the experimental variogram

fig1, ax1 = plt.subplots(figsize=[8, 8])
ax1.scatter(bin_center, dir_vario[0], label="Experimental 45°", color='blue')
ax1.scatter(bin_center, dir_vario[1], label="Experimental 135°", color='red')

# Plot the manually defined model 
manual_model = gs.Exponential(dim=2, var=25, len_scale=[25, 12.5], angles=angle)
manual_model.plot("vario_axis", axis=0, ax=ax1, x_max=max(bin_center), label="45°", linestyle='-', color='blue')
manual_model.plot("vario_axis", axis=1, ax=ax1, x_max=max(bin_center), label="135°", linestyle='-', color='red')

ax1.set_xlabel("Lag Distance")
ax1.set_ylabel("Semivariance")
ax1.legend()
ax1.set_title("Exponential Variogram")
plt.show()

Generated Field and Variogram Result:

Generated Field and Variogram by increasing the field size to 301

x = y = range(301)

My Question:

Why the variogram is not fitting for 45 degree direction for 100 by 100 area and increasing the size of area to 300 by 300 gives a better fit?

[LSchueler]

Hi Bishama, I'm sorry for replying so late!
I just reran your example with a different seed srf = gs.SRF(model, mean=50, seed=163456) and the fit looks much better. I guess with your seed of 123456 you created a freak field with an unusually long correlation length in the domain of 101 x 101. This area gets averaged out, if you increase your domain size to 301x301.
You should always keep in mind, that you are dealing with statistics and when working on only a single realization of a spatial random field you might just stumble upon one deviates from your expectations (at least on small domain sizes).

I hope that answers your question.

[Bishama]

Thank you for the response. That makes sense.