Add a Fourier generator for periodic spatial random fields

[LSchueler]

This is a first attempt to include a new generator based on the Fourier method. One of its desirable properties is the periodicity of the generated fields, which can be a necessity for lower boundary conditions in atmospherical physics.

Tests are still completely missing and I think the length scales of the fields are currently not working as expected, which is probably due to the rescaling of the Fourier modes. Before commiting any more to this, I'll wait for the comments from an atmospherical physicist.

Here is a simple script to show you how to currently use the new generator:

import numpy as np
import gstools as gs
import matplotlib.pyplot as pt

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

model = gs.Gaussian(dim=2, var=1, len_scale=30.0)
srf = gs.SRF(model, generator='Fourier', modes_truncation=[1.0, 1.0], modes_no=[80, 80], seed=3684903)

field = srf((x, y), mesh_type='structured')

srf.plot()
pt.show()

[mschlutow]

Something's not quite right. The shape of the histogram is not Gaussian and its variance not unity. Becomes even more crooked when more modes are added.

import numpy as np
import gstools as gs
import matplotlib.pyplot as pt

Nx=512
Ny=256

x = np.arange(Nx)
y = np.arange(Ny)

model = gs.Gaussian(dim=2, var=1, len_scale=100.0)
srf = gs.SRF(model, generator='Fourier', modes_truncation=[1.0, 1.0], modes_no=[80, 80], seed=3684903)

field = srf((x, y), mesh_type='structured')

srf.plot()
pt.show()

pt.hist(field.flatten('C'),bins=64)
pt.show()

[LSchueler]

Thanks for providing the script :+1: I guess the reason is that you only sample very few correlation lengths. I just ran your script, but with len_scale=1.0 and this is the result:

[mschlutow]

Thanks for the quick reply. I am afraid that it is not only the length scale that leads to uneven distribution.

model = gs.Gaussian(dim=2, var=1.0, len_scale=200.0)
srf = gs.SRF(model, generator='Fourier', modes_truncation=[1.0, 1.0], modes_no=[512, 256])

Why is it symmetrical?

[LSchueler]

Yeah, there is definitely something wrong with the length scales. Here is a script to estimate the variograms, which also gives you the actual length scale of a field. If you want to use your own field generator, simply ignore everything until line 21 and give your field to the variogram estimator as the field keyword in line 29. I have to see, if I find some time this evening to check the rescaling of the Fourier modes.

import numpy as np
import gstools as gs
import matplotlib.pyplot as pt

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

dim = 2
model = gs.Gaussian(dim=dim, var=1, len_scale=[5.0, 3.0])
# fourier method
srf = gs.SRF(
    model,
    generator='Fourier',
    modes_truncation=[1.0, 1.0],
    modes_no=[100, 60],
    seed=3684903,
)
# randomization method
# srf = gs.SRF(model, seed=3684903)

field = srf((x, y), mesh_type='structured')

x_max = 40.0
angle = 0.0

bins = np.arange(0, x_max, 2)

bin_center, dir_vario, counts = gs.vario_estimate(
    *((x, y), field, bins),
    direction=gs.rotated_main_axes(dim, angle),
    angles_tol=np.pi / 16,
    bandwidth=8,
    mesh_type="structured",
    return_counts=True,
)

model.fit_variogram(bin_center, dir_vario)
print(f'Fitted model: {model}')

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

ax1.scatter(bin_center, dir_vario[0], label="emp. vario in x-dir")
ax1.scatter(bin_center, dir_vario[1], label="emp. vario: in y-dir")
ax1.legend(loc="lower right")

model.plot("vario_axis", axis=0, ax=ax1, x_max=x_max, label="fit in x-dir")
model.plot("vario_axis", axis=1, ax=ax1, x_max=x_max, label="fit in y-dir")
ax1.set_title("Fitting an anisotropic model")

srf.plot(ax=ax2)
pt.show()

[LSchueler]

Yeah, I think the final problem to solve is the scaling of the Fourier modes. This is an example, where I manually figured out the rescaling for one specific set of parameters and the (anisotropic) length scales are estimated correctly.

[MuellerSeb]

Please rebase with main.

[LSchueler]

Please rebase with main.

Good idea :sweat_smile:

[mschlutow]

GSTools/examples/00_misc/06_fourier.py throws error

/home/mschlutow/GSTools/src/gstools/field/plot.py:407: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
  fig.show()
Traceback (most recent call last):

  File ~/anaconda3/envs/nansen/lib/python3.9/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File ~/GSTools/examples/00_misc/06_fourier.py:55
    modes = [np.arange(0, modes_cutoff[d], modes_delta[d]) for d in 2]

TypeError: 'int' object is not iterable```

[LSchueler]

Sorry, it was late last evening... I just fixed that problem.

[LSchueler]

Wow, this was a looong journey. But I think we are finally there. We did a lot of testing in the background and the Fourier method seems to be working robustly now. @MuellerSeb Do you have time for a review?

[MuellerSeb]

Thanks @LSchueler for the great work. Will have a look as soon as I find the time! Could you update the branch with main?

[LSchueler]

I merged the main branch and I also implemented the Fourier method in GSTools-Core, so included that import in GSTools, if GSTools-Core is available.

[LSchueler]

This is ready for the next review round.

[LSchueler]

I can't believe it, after nearly 1 1/2 years, this is ready! Thanks for the reviews.