Plurigaussian fields

README.md

@@ -43,6 +43,7 @@ GeoStatTools provides geostatistical tools for various purposes:
 - data normalization and transformation
 - many readily provided and even user-defined covariance models
 - metric spatio-temporal modelling
+- plurigaussian field simulations (PGS)
 - plotting and exporting routines
 
 
@@ -110,6 +111,7 @@ The documentation also includes some [tutorials][tut_link], showing the most imp
 - [Geographic Coordinates][tut8_link]
 - [Spatio-Temporal Modelling][tut9_link]
 - [Normalizing Data][tut10_link]
+- [Plurigaussian Field Generation (PGS)][tut11_link]
 - [Miscellaneous examples][tut0_link]
 
 The associated python scripts are provided in the `examples` folder.
@@ -321,6 +323,52 @@ yielding
 [kraichnan_link]: https://doi.org/10.1063/1.1692799
 
 
+## Plurigaussian Field Simulation (PGS)
+
+With PGS, more complex categorical (or discrete) fields can be created.
+
+
+### Example
+
+```python
+import gstools as gs
+import numpy as np
+import matplotlib.pyplot as plt
+
+N = [180, 140]
+
+x, y = range(N[0]), range(N[1])
+
+# we need 2 SRFs
+model = gs.Gaussian(dim=2, var=1, len_scale=5)
+srf = gs.SRF(model)
+field1 = srf.structured([x, y], seed=20170519)
+field2 = srf.structured([x, y], seed=19970221)
+
+# with L, we prescribe the categorical data and its relations
+# here, we use 2 categories separated by a rectangle.
+R = [40, 32]
+L = np.zeros(N)
+L[
+    N[0] // 2 - R[0] // 2 : N[0] // 2 + R[0] // 2,
+    N[1] // 2 - R[1] // 2 : N[1] // 2 + R[1] // 2,
+] = 1
+
+pgs = gs.PGS(2, [field1, field2], L)
+
+fig, axs = plt.subplots(1, 2)
+axs[0].imshow(L, cmap="copper")
+axs[1].imshow(pgs.P, cmap="copper")
+plt.tight_layout()
+plt.savefig("2d_pgs.png")
+plt.show()
+```
+
+<p align="center">
+<img src="https://raw.githubusercontent.com/GeoStat-Framework/GSTools/main/docs/source/pics/2d_pgs.png" alt="PGS" width="600px"/>
+</p>
+
+
 ## VTK/PyVista Export
 
 After you have created a field, you may want to save it to file, so we provide
@@ -393,6 +441,7 @@ You can contact us via <[email protected]>.
 [tut8_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/08_geo_coordinates/index.html
 [tut9_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/09_spatio_temporal/index.html
 [tut10_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/10_normalizer/index.html
+[tut11_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/11_plurigaussian/index.html
 [tut0_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/00_misc/index.html
 [cor_link]: https://en.wikipedia.org/wiki/Autocovariance#Normalization
 [vtk_link]: https://www.vtk.org/

docs/source/conf.py

@@ -298,6 +298,7 @@ def setup(app):
         "../../examples/08_geo_coordinates/",
         "../../examples/09_spatio_temporal/",
         "../../examples/10_normalizer/",
+        "../../examples/11_plurigaussian/",
     ],
     # path where to save gallery generated examples
     "gallery_dirs": [
@@ -312,6 +313,7 @@ def setup(app):
         "examples/08_geo_coordinates/",
         "examples/09_spatio_temporal/",
         "examples/10_normalizer/",
+        "examples/11_plurigaussian/",
     ],
     # Pattern to search for example files
     "filename_pattern": r"\.py",

docs/source/index.rst

@@ -33,6 +33,7 @@ GeoStatTools provides geostatistical tools for various purposes:
 - data normalization and transformation
 - many readily provided and even user-defined covariance models
 - metric spatio-temporal modelling
+- plurigaussian field simulations (PGS)
 - plotting and exporting routines
 
 
@@ -164,6 +165,7 @@ showing the most important use cases of GSTools, which are
 - `Geographic Coordinates <examples/08_geo_coordinates/index.html>`__
 - `Spatio-Temporal Modelling <examples/09_spatio_temporal/index.html>`__
 - `Normalizing Data <examples/10_normalizer/index.html>`__
+- `Plurigaussian Field Generation (PGS) <examples/11_plurigaussian/index.html>`__
 - `Miscellaneous examples <examples/00_misc/index.html>`__
 
 
@@ -391,6 +393,54 @@ yielding
    :align: center
 
 
+Plurigaussian Field Simulation (PGS)
+====================================
+
+With PGS, more complex categorical (or discrete) fields can be created.
+
+
+Example
+-------
+
+.. code-block:: python
+
+   import gstools as gs
+   import numpy as np
+   import matplotlib.pyplot as plt
+
+   N = [180, 140]
+
+   x, y = range(N[0]), range(N[1])
+
+   # we need 2 SRFs
+   model = gs.Gaussian(dim=2, var=1, len_scale=5)
+   srf = gs.SRF(model)
+   field1 = srf.structured([x, y], seed=20170519)
+   field2 = srf.structured([x, y], seed=19970221)
+
+   # with L, we prescribe the categorical data and its relations
+   # here, we use 2 categories separated by a rectangle.
+   R = [40, 32]
+   L = np.zeros(N)
+   L[
+       N[0] // 2 - R[0] // 2 : N[0] // 2 + R[0] // 2,
+       N[1] // 2 - R[1] // 2 : N[1] // 2 + R[1] // 2,
+   ] = 1
+
+   pgs = gs.PGS(2, [field1, field2], L)
+
+   fig, axs = plt.subplots(1, 2)
+   axs[0].imshow(L, cmap="copper")
+   axs[1].imshow(pgs.P, cmap="copper")
+   plt.tight_layout()
+   plt.savefig("2d_pgs.png")
+   plt.show()
+
+.. image:: https://raw.githubusercontent.com/GeoStat-Framework/GSTools/main/docs/source/pics/2d_pgs.png
+   :width: 600px
+   :align: center
+
+
 VTK/PyVista Export
 ==================
 

docs/source/tutorials.rst

@@ -22,6 +22,7 @@ explore its whole beauty and power.
    examples/08_geo_coordinates/index
    examples/09_spatio_temporal/index
    examples/10_normalizer/index
+   examples/11_plurigaussian/index
    examples/00_misc/index
 
 .. only:: html

examples/11_plurigaussian/00_simple.py

@@ -0,0 +1,72 @@
+"""
+A First and Simple Example
+--------------------------
+
+As a first example, we will create a two dimensional plurigaussian field
+(PGS). Thus, we need two spatial random fields(SRF) and on top of that, we
+need a field describing the categorical data and its spatial relation.
+We will start off by creating the two SRFs with a Gaussian variogram, which
+makes the fields nice and smooth. But before that, we will import all
+necessary libraries and define a few variables, like the number of grid
+cells in each dimension.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+dim = 2
+# no. of cells in both dimensions
+N = [180, 140]
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+
+###############################################################################
+# In this first example we will use the same geostatistical parameters for
+# both fields for simplicity. Thus, we can use the same SRF instance for the
+# two fields.
+
+model = gs.Gaussian(dim=dim, var=1, len_scale=10)
+srf = gs.SRF(model)
+field1 = srf.structured([x, y], seed=20170519)
+field2 = srf.structured([x, y], seed=19970221)
+
+###############################################################################
+# Now, we will create the field describing the categorical data. For now, we
+# will only have two categories and we will address them by the integers 0 and 1.
+# We start off by creating a matrix of 0s from which we will select a rectangle
+# and fill that with 1s.
+# This field does not have to match the shape of the SRFs.
+
+M = [200, 160]
+
+# size of the rectangle
+R = [40, 32]
+
+L = np.zeros(M)
+L[
+    M[0] // 2 - R[0] // 2 : M[0] // 2 + R[0] // 2,
+    M[1] // 2 - R[1] // 2 : M[1] // 2 + R[1] // 2,
+] = 1
+
+###############################################################################
+# With the two SRFs and `L` ready, we can create our first PGS.
+
+pgs = gs.PGS(dim, [field1, field2], L)
+
+###############################################################################
+# Finally, we can plot the PGS, but we will also show the field `L` and the two
+# original Gaussian fields.
+
+fig, axs = plt.subplots(2, 2)
+
+axs[0, 0].imshow(field1, cmap="copper")
+axs[0, 1].imshow(field2, cmap="copper")
+axs[1, 0].imshow(L, cmap="copper")
+axs[1, 1].imshow(pgs.P, cmap="copper")
+
+# For more information on Plurigaussian fields and how they naturally extend
+# truncated Gaussian fields, we can recommend the book
+# [Plurigaussian Simulations in Geosciences](https://doi.org/10.1007/978-3-642-19607-2)

examples/11_plurigaussian/01_pgs.py

@@ -0,0 +1,118 @@
+"""
+Understanding PGS
+-----------------
+
+In this example we want to try to understand how exactly PGS are generated
+and how to influence them with the categorical field.
+First of all, we will set everything up very similar to the first example.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+dim = 2
+# no. of cells in both dimensions
+N = [100, 80]
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+
+###############################################################################
+# In this example we will use different geostatistical parameters for the
+# SRFs. We will create fields with a strong anisotropy, and on top of that they
+# will both be rotated.
+
+model1 = gs.Gaussian(dim=dim, var=1, len_scale=[20, 1], angles=np.pi / 8)
+srf1 = gs.SRF(model1, seed=20170519)
+field1 = srf1.structured([x, y])
+model2 = gs.Gaussian(dim=dim, var=1, len_scale=[1, 20], angles=np.pi / 4)
+srf2 = gs.SRF(model2, seed=19970221)
+field2 = srf2.structured([x, y])
+field1 += 5.0
+
+###############################################################################
+# Internally, each field's values are mapped along an axis, which can be nicely
+# visualized with a scatter plot. We can easily do that by flattening the 2d
+# field values and simply use matplotlib's scatter plotting functionality.
+# The x-axis shows field1's values and the y-axis shows field2's values.
+
+plt.scatter(field1.flatten(), field2.flatten(), s=0.1)
+
+###############################################################################
+# This mapping always has a multivariate Gaussian distribution and this is also
+# the field on which we define our categorical data `L` and their relations to each
+# other. Before providing further explanations, we will create `L`, which again
+# will have only two categories, but this time we will not prescribe a rectangle,
+# but a circle.
+
+# no. of grid cells of field L
+M = [51, 41]
+# we need the indices of L later
+x_L = np.arange(M[0])
+y_L = np.arange(M[1])
+
+# radius of circle
+R = 7
+
+L = np.zeros(M)
+mask = (x_L[:, np.newaxis] - M[0] // 2) ** 2 + (
+    y_L[np.newaxis, :] - M[1] // 2
+) ** 2 < R**2
+L[mask] = 1
+
+###############################################################################
+# Now, we look at every point in the scatter plot (which shows the field values)
+# shown above and map the categorical values of `L` to the positions (variables
+# `x` and `y` defined above) of these field values. This is probably much easier
+# to understand with a plot.
+# First, we calculate the indices of the L field, which we need for manually
+# plotting L together with the scatter plot. Normally they are computed internally.
+
+l = np.floor(field1.min()) - 1
+h = np.ceil(field1.max()) + 1
+m = (h + l) / 2.0
+dist = max(np.abs(h - m), np.abs(l - m))
+x_mean = field1.mean()
+x_l = np.linspace(
+    x_mean - dist,
+    x_mean + dist,
+    L.shape[0],
+)
+
+l = np.floor(field1.min()) - 1
+h = np.ceil(field1.max()) + 1
+m = (h + l) / 2.0
+dist = max(np.abs(h - m), np.abs(l - m))
+y_mean = field2.mean()
+y_l = np.linspace(
+    y_mean - dist,
+    y_mean + dist,
+    L.shape[1],
+)
+
+###############################################################################
+# We also compute the actual PGS now, to also plot that.
+
+pgs = gs.PGS(dim, [field1, field2], L)
+
+###############################################################################
+# And now to some plotting. Unfortunately, matplotlib likes to mess around with
+# the aspect ratios of the plots, so the left panel is a bit stretched.
+
+fig, axs = plt.subplots(2, 2)
+axs[0, 0].imshow(field1, cmap="copper")
+axs[0, 1].imshow(field2, cmap="copper")
+axs[1, 0].scatter(field1.flatten(), field2.flatten(), s=0.1, color="C0")
+axs[1, 0].pcolormesh(x_l, y_l, L.T, alpha=0.3)
+
+axs[1, 1].imshow(pgs.P, cmap="copper")
+plt.show()
+
+###############################################################################
+# The black areas show the category 0 and the orange areas show category 1. We
+# see that the majority of all points in the scatter plot are within the
+# yellowish circle, thus the orange areas are larger than the black ones. The
+# strong anisotropy and the rotation of the fields create these interesting
+# patterns which remind us of fractures in the subsurface.