Artefacts when using Spherical model for SRF generation

[rhugonnet]

Hi there,
First, thanks for the great work!

I use a lot of Spherical models in my work, and I can't seem to create a smooth random field. Might be related to #126.
Version is 1.3.3.

Reproducible example with both short and long-range covariance models:

import gstools as gs
import numpy as np

x = np.arange(0, 100)
y = np.arange(0, 100)

model_s = gs.Spherical(dim=2, var=1, len_scale=5)
model_l = gs.Spherical(dim=2, var=1, len_scale=80)

srf_s = gs.SRF(model_s, mode_no=1000)
srf_l = gs.SRF(model_l, mode_no=1000)

field_s = srf_s((x, y), mesh_type='structured')
field_l = srf_l((x, y), mesh_type='structured')

srf_s.plot()
srf_l.plot()

The short-range looks OK, although there some grains:

The long-range doesn't look so good. We see the long-range correlation, but there's plenty of small grainy artefact:

It's nothing like the Gaussian-based fields, for sure!

I played with the mode_no (between 100 and 10,000) and upscaling (trying out "coarse_graining") arguments, but unfortunately without real success... 😞

[LSchueler]

Hi!

The spherical model is much rougher than the Gaussian model. That's what you are seeing in your plots. The spatial resoultion you are using cannot resolve the roughness. Try increasing the spatial resolution of your domain, e.g.

x = np.arange(0, 100, 0.1)
y = np.arange(0, 100, 0.1)

You can also try using a different kind of plot, which does a better job at interpolating the colours for rougher fields. We provide some convenience plotting functionality, but nothing too fleshed out.

[MuellerSeb]

This perfectly demonstrates, that the roughness-information is a=1 for the spherical model and that means, it is not differentiable (only for a=2) and that means, you can zoom in as much as you want, it will never get smooth.

[rhugonnet]

@LSchueler @MuellerSeb Thanks for the quick answers! I think issue can be closed then. Maybe a note in the documentation on less differentiables models would help first users grasp this concept?

Do you know if there any other method than Hesse et al. (2014) that might perform better at creating a random field for less differentiable models, maybe unconditional Gaussian simulation for example?

[LSchueler]

What exactly do you mean by "performing better"? I hope this doesn't sound too much like marketing BS, but it doesn't get much better than this 😉
As I suggested, you could have a look at different visualisation functions or even libraries. Maybe contour plots already improve the visualisation enough.
Here is a quick contour plot for your second example, where you can still see the roughness, but the ugly patches from the colour mesh interpolation is gone:

[MuellerSeb]

Maybe a note in the documentation on less differentiables models would help first users grasp this concept?

I was already thinking about adding an attribute to the covmodel class holding the roghness-information. But calculating it numerically is quite hard. Need to think about it.

[rhugonnet]

Ahah it doesn't, but I'm not sure my question was clear 😅! By "performing better", I meant entirely other approaches to the simulation of correlated fields. If I'm not mistaken, GSTools currently includes only that of Hesse et al. (2014) based on a Fourier method (i.e. cos/sin) to simulate the random field.
How will the Hesse et al. (2014) method compare to, for example, an unconditional/conditional Gaussian simulation based on kriging to simulate the correlated random field (usually using a sequential path of propagation)? Intuitively, I assume kriging equations might be less affected by a covariance model that is rougher than the Fourier method is. But maybe I'm missing something here!

[MuellerSeb]

From my point of view, a sequential gaussian simulation (SGS) should provide similar results but at much higher computational cost. And if the resulting field from SGS doesn't reflect the roughness of the variogram model used, this would be a bad thing.

[rhugonnet]

Thanks for the clarification!
I expected the roughness to be a numerical artefact of generation method of the correlation field (as it is not differentiable, it is hard to reach a Fourier representation), rather than something inherent to the variogram (i.e. a correlation field that would perfectly match the variogram).
I might do some Gaussian simulation tests at some point to better understand this 🤔. I think all my questions have been answered, thanks again!