-
Notifications
You must be signed in to change notification settings - Fork 69
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add a Fourier generator for periodic spatial random fields #302
base: main
Are you sure you want to change the base?
Conversation
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.
|
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 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() |
Please rebase with main. |
Good idea 😅 |
Now, the modes can be calculated internally from a rel. cutoff value and the periodicity.
|
Sorry, it was late last evening... I just fixed that problem. |
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: