Merge pull request #302 from GeoStat-Framework/fourier-gen

Add a Fourier generator for periodic spatial random fields

Author: LSchueler

.github/workflows/main.yml

@@ -48,9 +48,9 @@ jobs:
         run: |
           python -m pylint src/gstools/
 
-      - name: cython-lint check
-        run: |
-          cython-lint src/gstools/
+      #- name: cython-lint check
+        #run: |
+          #cython-lint src/gstools/
 
   build_wheels:
     name: wheels for ${{ matrix.os }}

README.md

@@ -34,7 +34,7 @@
 <img align="right" width="450" src="https://raw.githubusercontent.com/GeoStat-Framework/GSTools/main/docs/source/pics/demonstrator.png" alt="">
 
 GeoStatTools provides geostatistical tools for various purposes:
-- random field generation
+- random field generation, including periodic boundaries
 - simple, ordinary, universal and external drift kriging
 - conditioned field generation
 - incompressible random vector field generation

docs/source/index.rst

@@ -24,7 +24,7 @@ Purpose
 
 GeoStatTools provides geostatistical tools for various purposes:
 
-- random field generation
+- random field generation, including periodic boundaries
 - simple, ordinary, universal and external drift kriging
 - conditioned field generation
 - incompressible random vector field generation

examples/01_random_field/08_fourier.py

@@ -0,0 +1,44 @@
+"""
+Generating a Simple Periodic Random Field
+-----------------------------------------
+
+In this simple example we are going to learn how to generate periodic spatial
+random fields. The Fourier method comes naturally with the property of
+periodicity, so we'll use it to create the random field.
+"""
+
+import numpy as np
+
+import gstools as gs
+
+# We start off by defining the spatial grid. For the sake of simplicity, we
+# use a square domain. We set the optional argument `endpoint` to `False`, to
+# not make the domain in each dimension one grid cell larger than the
+# periodicity.
+L = 500.0
+x = np.linspace(0, L, 256, endpoint=False)
+y = np.linspace(0, L, 128, endpoint=False)
+
+# Now, we create a Gaussian covariance model with a correlation length which is
+# roughly half the size of the grid.
+model = gs.Gaussian(dim=2, var=1, len_scale=200)
+
+# Next, we hand the cov. model to the spatial random field class `SRF`
+# and set the generator to `"Fourier"`. The argument `period` is set to the
+# domain size. If only a single number is given, the same periodicity is
+# applied in each dimension, as shown in this example. The `mode_no` argument
+# sets the number of Fourier modes. If only an integer is given, that number
+# of modes is used for all dimensions.
+srf = gs.SRF(
+    model,
+    generator="Fourier",
+    period=L,
+    mode_no=32,
+    seed=1681903,
+)
+
+# Now, we can calculate the field with the given parameters.
+srf((x, y), mesh_type="structured")
+
+# GSTools has a few simple visualization methods built in.
+srf.plot()

examples/01_random_field/09_fourier_trans.py

@@ -0,0 +1,49 @@
+"""
+Generating a Transformed Periodic Random Field
+----------------------------------------------
+
+Building on the precious example, we are now going to generate periodic
+spatial random fields with a transformation applied, resulting in a level set.
+"""
+
+import numpy as np
+
+import gstools as gs
+
+# We start off by defining the spatial grid. As in the previous example, we do
+# not want to include the endpoints.
+L = np.array((500, 400))
+x = np.linspace(0, L[0], 300, endpoint=False)
+y = np.linspace(0, L[1], 200, endpoint=False)
+
+# Instead of using a Gaussian covariance model, we will use the much rougher
+# exponential model and we will introduce an anisotropy by using two different
+# length scales in the x- and y-directions
+model = gs.Exponential(dim=2, var=2, len_scale=[80, 20])
+
+# Same as before, we set up the spatial random field. But this time, we will
+#  use a periodicity which is equal to the domain size in x-direction, but
+# half the domain size in y-direction. And we will use different `mode_no` for
+# the different dimensions.
+srf = gs.SRF(
+    model,
+    generator="Fourier",
+    period=[L[0], L[1] / 2],
+    mode_no=[30, 20],
+    seed=1681903,
+)
+# and compute it on our spatial domain
+srf((x, y), mesh_type="structured")
+
+# With the field generated, we can now apply transformations starting with a
+# discretization of the field into 4 different values
+thresholds = np.linspace(np.min(srf.field), np.max(srf.field), 4)
+srf.transform("discrete", store="transform_discrete", values=thresholds)
+srf.plot("transform_discrete")
+
+# This is already a nice result, but we want to pronounce the peaks of the
+# field. We can do this by applying a log-normal transformation on top
+srf.transform(
+    "lognormal", field="transform_discrete", store="transform_lognormal"
+)
+srf.plot("transform_lognormal")

examples/01_random_field/README.rst

@@ -11,6 +11,12 @@ semi-variogram. This is done by using the so-called randomization method.
 The spatial random field is represented by a stochastic Fourier integral
 and its discretised modes are evaluated at random frequencies.
 
+In case you want to generate spatial random fields with periodic boundaries,
+you can use the so-called Fourier method. See the corresponding examples for
+how to do that. The spatial random field is represented by a stochastic
+Fourier integral and its discretised modes are evaluated at equidistant
+frequencies.
+
 GSTools supports arbitrary and non-isotropic covariance models.
 
 Examples