Artifacts in generated random fields

[flydream0428]

On the online document, there is a simple example of random field generation.

import gstools as gs
# structured field with a size 100x100 and a grid-size of 1x1
x = y = range(100)
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model)
srf((x, y), mesh_type='structured')
srf.plot()

The code gives a perfect field like the one below,

However, if I change the model to the Exponential model using

model = gs.Exponential(dim=2, var=1, len_scale=10)

I will get very granular field,

Although, we expect different models will lead to different covariance, I don't think the difference can be so significant.

[MuellerSeb]

Hi there! Nice to see people question stuff like this. That's how it should be 😉

But I think you fell for a red herring here. To me, this looks absolutely fine. We could investigate two things here:

1. re-estimate the empirical variogram from the generated fields:

   import gstools as gs
   # structured field with a size 100x100 and a grid-size of 1x1
   x = y = range(100)
   model = gs.Gaussian(dim=2, var=1, len_scale=10)
   srf = gs.SRF(model, seed=3141592)
   srf((x, y), mesh_type="structured")
   
   bins, gamma = gs.vario_estimate((x, y), srf.field, mesh_type="structured")
   fit_model = gs.Gaussian(dim=2)
   fit_model.fit_variogram(bins, gamma)
   
   ax = fit_model.plot(x_max=max(bins), label="fitted Gaussian model")
   ax.scatter(bins, gamma, label="estimated variogram")
   model.plot(ax=ax, x_max=max(bins), label="original Gaussian model")
   
   model = gs.Exponential(dim=2, var=1, len_scale=10)
   srf.model = model
   srf((x, y), mesh_type="structured")
   
   fit_model = gs.Exponential(dim=2)
   bins, gamma = gs.vario_estimate((x, y), srf.field, mesh_type="structured")
   fit_model.fit_variogram(bins, gamma)
   
   ax = fit_model.plot(x_max=max(bins), label="fitted Exponential model")
   ax.scatter(bins, gamma, label="estimated variogram")
   model.plot(ax=ax, x_max=max(bins), label="original Exponential model")

   
   
   • as you see, the original variograms are pretty well matching the estimated ones
   • this means, that the fields indeed follow the input variogram models
2. A crucial information for the look of a field is the roughness information for a vairogram model. For the associated correlation function $\rho(r)$ the roughness information $\alpha$ is determined by:
   • $\rho(r)\approx 1 - k\cdot r^\alpha \quad \text{as} \quad r\to 0$
   • $\alpha$ is between 0 (very rough, pure nugget model) and 2 (very smooth like the gaussian)
   • This leads to the following observation:
     • Gaussian $\rho(r)=\exp(-r^2/\ell^2)$ has $\alpha=2$
     • Exponential $\rho(r)=\exp(-r/\ell)$ has $\alpha=1$
     • So these two models have very different roughness properties which is in line with the observation
   • you can play around with the Matern model, which has a shape parameter $\nu$ which is directly connected to the associated roughness information $\alpha=\min(\nu, 2)$
     • this is in line with the definition of the matern model, which is identical with the Exponential model for $\nu=1$ and converges to the Gaussian model for $\nu\to\infty$

I hope this helps a bit to understand the look of the generated fields and to get a better feeling for the models.

Cheers,
Sebastian

[flydream0428]

Hi Sebastian,

Thanks a lot for your help. I agree that the variogram you have used to verify the field confirms the general pattern is right.
However, I still believe there is a glitch here. I did more experiments to test this. I have compared the results from PyEMU and GsTools.
I have increased the range here to 10000 to make it more obvious.

x = np.linspace(minx, maxx, 150)
y = np.linspace(miny, maxy, 150)
model = gs.Exponential(dim=2, var=1, len_scale=10000)
srf = gs.SRF(model)
field = srf.structured([x, y])
srf.plot()

import pyemu 
v = pyemu.utils.geostats.ExpVario(a=10000, contribution=1.0)
pgs = pyemu.utils.geostats.GeoStruct(variograms=v, nugget=0.0)
delx = np.ones(150)*100
dely = delx
ss = pyemu.utils.geostats.SpecSim2d(delx,dely,pgs)
arr = np.squeeze(ss.draw_arrays(num_reals=1))
plt.imshow(arr)
plt.colorbar(shrink=1)

If you compare the two figures, the figure from GsTools has some strips in it. If we increase the range to 1e5, it's even more obvious.

Cheers,
Tao

[MuellerSeb]

That is an interesting picture 😉

What we are seeing here, are artifacts of the underlying Randomization Method (see Heße et al. 2014 for details).

The idea of the randomization method is to approximate the field in fourier space and transforming it back (there are lots of advantages in contrast to the sequential gaussian simulation approach and that is why we choose to use this method).

The approximation is cut of by a maximal number of sin and cos modes that are summed. If the number of modes is to low (1000 by default in GSTools), you can encounter the exact phenomena you just reproduced. That means, you actually see the superpositioned sin and cos waves.

To overcome these artifacts, you can increase the number of modes the generator uses directly in the SRF class (kwargs are forwarded to the generator used):

x = np.linspace(0, 15000, 150)
y = np.linspace(0, 15000, 150)
model = gs.Exponential(dim=2, var=1, len_scale=10000)
srf = gs.SRF(model, seed=3141592, mode_no=15000)
field = srf.structured([x, y])
srf.plot(contour_plot=False)
plt.show()

(in addition I disabled the contour plot to be closer to the imshow look and feel).

Although it is mentioned in the documentation of the RandMeth class, this feature may should get an example in the Random Field Generation tutorial section.

Hope that helps!
Cheers, Sebastian

[MuellerSeb]

The take-away message here is:

If the length scale is not at least a magnitude smaller than the field size (~factor 10), you need to increase the number of modes in the Randomization method.

[flydream0428]

Thanks a lot. This is very helpful.