Merge pull request #370 from GeoStat-Framework/pgs

Plurigaussian fields

Author: LSchueler

README.md

@@ -43,6 +43,7 @@ GSTools 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,51 @@ 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 `lithotypes`, we prescribe the categorical data and its relations
+# here, we use 2 categories separated by a rectangle.
+rect = [40, 32]
+lithotypes = np.zeros(N)
+lithotypes[
+    N[0] // 2 - rect[0] // 2 : N[0] // 2 + rect[0] // 2,
+    N[1] // 2 - rect[1] // 2 : N[1] // 2 + rect[1] // 2,
+] = 1
+
+pgs = gs.PGS(2, [field1, field2])
+P = pgs(lithotypes)
+
+fig, axs = plt.subplots(1, 2)
+axs[0].imshow(lithotypes, cmap="copper")
+axs[1].imshow(P, cmap="copper")
+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 +440,7 @@ You can contact us via <info@geostat-framework.org>.
 [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

@@ -300,6 +300,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": [
@@ -314,6 +315,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 @@ GSTools 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,53 @@ 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 `lithotypes`, we prescribe the categorical data and its relations
+   # here, we use 2 categories separated by a rectangle.
+   rect = [40, 32]
+   lithotypes = np.zeros(N)
+   lithotypes[
+       N[0] // 2 - rect[0] // 2 : N[0] // 2 + rect[0] // 2,
+       N[1] // 2 - rect[1] // 2 : N[1] // 2 + rect[1] // 2,
+   ] = 1
+
+   pgs = gs.PGS(2, [field1, field2])
+   P = pgs(lithotypes)
+
+   fig, axs = plt.subplots(1, 2)
+   axs[0].imshow(lithotypes, cmap="copper")
+   axs[1].imshow(P, cmap="copper")
+   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/pics/2d_pgs.png

@@ -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

docs/source/pics/3d_pgs.png

@@ -20,6 +20,6 @@
 ###############################################################################
 # The result displayed with Paraview:
 #
-# .. image:: https://raw.githubusercontent.com/GeoStat-Framework/GeoStat-Framework.github.io/master/img/paraview.png
+# .. image:: ../../pics/paraview.png
 #    :width: 400px
 #    :align: center

docs/source/tutorials.rst

@@ -53,6 +53,6 @@
 ###############################################################################
 # The result should look like this:
 #
-# .. image:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_pyvista_cut.png
+# .. image:: ../../pics/GS_pyvista_cut.png
 #    :width: 400px
 #    :align: center

examples/00_misc/01_export.py

@@ -59,6 +59,6 @@
 ###############################################################################
 # The result should look like this:
 #
-# .. image:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_3d_vector_field.png
+# .. image:: ../../pics/GS_3d_vector_field.png
 #    :width: 400px
 #    :align: center

examples/01_random_field/06_pyvista_support.py

@@ -0,0 +1,75 @@
+"""
+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 lithotypes 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.
+
+centroid = [200, 160]
+
+# size of the rectangle
+rect = [40, 32]
+
+lithotypes = np.zeros(centroid)
+lithotypes[
+    centroid[0] // 2 - rect[0] // 2 : centroid[0] // 2 + rect[0] // 2,
+    centroid[1] // 2 - rect[1] // 2 : centroid[1] // 2 + rect[1] // 2,
+] = 1
+
+###############################################################################
+# With the two SRFs and the L-field ready, we can create our first PGS. First, we
+# set up an instance of the PGS class and then we are ready to calculate the
+# field by calling the instance and handing over the L-field.
+
+pgs = gs.PGS(dim, [field1, field2])
+P = pgs(lithotypes)
+
+###############################################################################
+# Finally, we can plot the PGS, but we will also show the L-field and the two
+# original Gaussian fields.
+
+fig, axs = plt.subplots(2, 2)
+
+axs[0, 0].imshow(field1, cmap="copper", origin="lower")
+axs[0, 1].imshow(field2, cmap="copper", origin="lower")
+axs[1, 0].imshow(lithotypes, cmap="copper", origin="lower")
+axs[1, 1].imshow(P, cmap="copper", origin="lower")
+
+# 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/04_vector_field/01_3d_vector_field.py

@@ -0,0 +1,93 @@
+"""
+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 `lithotypes` and their
+# relations to each other. Before providing further explanations, we will
+# create the lithotypes field, which again will have only two categories, but
+# this time we will not prescribe a rectangle, but a circle.
+
+# no. of grid cells of L-field
+M = [51, 41]
+# we need the indices of `lithotypes` later
+x_lith = np.arange(M[0])
+y_lith = np.arange(M[1])
+
+# radius of circle
+radius = 7
+
+lithotypes = np.zeros(M)
+mask = (x_lith[:, np.newaxis] - M[0] // 2) ** 2 + (
+    y_lith[np.newaxis, :] - M[1] // 2
+) ** 2 < radius**2
+lithotypes[mask] = 1
+
+###############################################################################
+# We can compute the actual PGS now. As a second step, we use a helper function
+# to recalculate the axes on which the lithotypes are defined. Normally, this
+# is handled internally. But in order to show the scatter plot together with
+# the lithotypes, we need the axes here.
+
+pgs = gs.PGS(dim, [field1, field2])
+P = pgs(lithotypes)
+
+x_lith, y_lith = pgs.calc_lithotype_axes(lithotypes.shape)
+
+###############################################################################
+# 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", origin="lower")
+axs[0, 1].imshow(field2, cmap="copper", origin="lower")
+axs[1, 0].scatter(field1.flatten(), field2.flatten(), s=0.1, color="C0")
+axs[1, 0].pcolormesh(x_lith, y_lith, lithotypes.T, alpha=0.3, cmap="copper")
+
+axs[1, 1].imshow(P, cmap="copper", origin="lower")
+
+###############################################################################
+# 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.