diff --git a/README.md b/README.md
index b90ad129..fe5b206e 100644
--- a/README.md
+++ b/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()
+```
+
+
+
+
+
+
## 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 .
[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/
diff --git a/docs/source/conf.py b/docs/source/conf.py
index e89928fc..22349910 100644
--- a/docs/source/conf.py
+++ b/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",
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 7abb1e9e..1aff309b 100644
--- a/docs/source/index.rst
+++ b/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 `__
- `Spatio-Temporal Modelling `__
- `Normalizing Data `__
+- `Plurigaussian Field Generation (PGS) `__
- `Miscellaneous examples `__
@@ -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
==================
diff --git a/docs/source/pics/2d_pgs.png b/docs/source/pics/2d_pgs.png
new file mode 100644
index 00000000..61551681
Binary files /dev/null and b/docs/source/pics/2d_pgs.png differ
diff --git a/docs/source/pics/3d_pgs.png b/docs/source/pics/3d_pgs.png
new file mode 100644
index 00000000..5a00475a
Binary files /dev/null and b/docs/source/pics/3d_pgs.png differ
diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst
index 3c25f597..fa8b45c9 100644
--- a/docs/source/tutorials.rst
+++ b/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
diff --git a/examples/11_plurigaussian/00_simple.py b/examples/11_plurigaussian/00_simple.py
new file mode 100644
index 00000000..7e4b2dfc
--- /dev/null
+++ b/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)
diff --git a/examples/11_plurigaussian/01_pgs.py b/examples/11_plurigaussian/01_pgs.py
new file mode 100644
index 00000000..0b0a7f9c
--- /dev/null
+++ b/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.
diff --git a/examples/11_plurigaussian/02_spatial_relations.py b/examples/11_plurigaussian/02_spatial_relations.py
new file mode 100644
index 00000000..320e65ca
--- /dev/null
+++ b/examples/11_plurigaussian/02_spatial_relations.py
@@ -0,0 +1,109 @@
+"""
+Controlling Spatial Relations
+-----------------------------
+
+In this example we will try to understand how we can
+influence the spatial relationships of the different
+categories with the field `L`. For simplicity, we will
+start very similarly to the very 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])
+
+###############################################################################
+# Again, we will use the same parameters for both 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 prepare the `L` field, which will be a bit more complicated this
+# time. First, we will create a triangle. Next, we will create two rectangles
+# touching each other along one of their edges and both being directly above the
+# triangle, but without touching it directly. Finally, we will create a few
+# very narrow rectangles, which will not touch any other category shapes. The
+# implementation details are not very interesting, and can be skipped.
+
+# no. of grid cells of field L
+M = [60, 50]
+
+# size of the rectangles
+R = [10, 8]
+
+# positions of some of the shapes for concise indexing
+S1 = [1, -9]
+S2 = [-5, 3]
+S3 = [-5, -5]
+
+L = np.zeros(M)
+# a small upper triangular helper matrix to create the triangle
+tri = np.triu(np.ones((R[0], R[0])))
+# the triangle
+L[
+ M[0] // 2 + S1[0] : M[0] // 2 + S1[0] + R[0],
+ M[1] // 2 + S1[1] : M[1] // 2 + S1[1] + R[0],
+] = tri
+# the first rectangle
+L[
+ M[0] // 2 + S2[0] - R[0] // 2 : M[0] // 2 + S2[0] + R[0] // 2,
+ M[1] // 2 + S2[1] - R[1] // 2 : M[1] // 2 + S2[1] + R[1] // 2,
+] = 2
+# the second rectangle
+L[
+ M[0] // 2 + S3[0] - R[0] // 2 : M[0] // 2 + S3[0] + R[0] // 2,
+ M[1] // 2 + S3[1] - R[1] // 2 : M[1] // 2 + S3[1] + R[1] // 2,
+] = 3
+# some very narrow rectangles
+for i in range(4):
+ L[
+ M[0] // 2 + S1[0] : M[0] // 2 + S1[0] + R[0],
+ M[1] // 2
+ + S1[1]
+ + R[1]
+ + 3
+ + 2 * i : M[1] // 2
+ + S1[1]
+ + R[1]
+ + 4
+ + 2 * i,
+ ] = (
+ 4 + i
+ )
+
+###############################################################################
+# With the two SRFs and `L` ready, we can create the PGS.
+pgs = gs.PGS(dim, [field1, field2], L)
+
+###############################################################################
+# And now the plotting of the two Gaussian fields, `L`, and the PGS.
+
+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")
+
+###############################################################################
+# We can see that the two upper light and medium brown rectangles both fill up
+# large and rather smooth areas of the PGS. And they share very long common
+# borders due to the fact that these categories touch each other along one of
+# their edges. The next large area is the dark brown of the triangle. This
+# category is always very close to the light brown areas, but only sometimes
+# close to the medium brown areas, as they only share small parts in close
+# proximity to each other. Finally, we have the four stripes. They create
+# distorted stripes in the PGS. The lighter they get, the less area they fill.
+# This is due to the fact that their area is not only relatively small, but
+# also because they are increasingly further away from the center of `L`.
diff --git a/examples/11_plurigaussian/03_correlations.py b/examples/11_plurigaussian/03_correlations.py
new file mode 100644
index 00000000..98008a2d
--- /dev/null
+++ b/examples/11_plurigaussian/03_correlations.py
@@ -0,0 +1,69 @@
+"""
+Understanding the Influence of Variograms
+-----------------------------------------
+
+Up until now, we have only used very smooth Gaussian variograms for the
+underlying spatial random fields. Now, we will combine a smooth Gaussian
+field with a much rougher exponential field. This example should feel
+familiar, if you had a look at the previous examples.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+dim = 2
+# no. of cells in both dimensions
+N = [200, 200]
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+
+###############################################################################
+# Now, we generate fields with a Gaussian and an Exponential variogram.
+
+model1 = gs.Gaussian(dim=dim, var=1, len_scale=[50, 25])
+srf1 = gs.SRF(model1)
+field1 = srf1.structured([x, y], seed=20170519)
+model2 = gs.Exponential(dim=dim, var=1, len_scale=[40, 40])
+srf2 = gs.SRF(model2)
+field2 = srf2.structured([x, y], seed=19970221)
+
+###############################################################################
+# The `L` field will consist of a circle which contains one category and the
+# surrounding is the second category.
+
+# no. of grid cells of field L
+M = [200, 200]
+
+# radius of circle
+R = 25
+
+x_L = np.arange(M[0])
+y_L = np.arange(M[1])
+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
+
+###############################################################################
+# With the two SRFs and `L` ready, we can create the PGS.
+pgs = gs.PGS(dim, [field1, field2], L)
+
+###############################################################################
+# And now the plotting of the two Gaussian fields, `L`, and the PGS.
+
+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")
+plt.show()
+
+###############################################################################
+# In this PGS, we can see two different spatial structures combined. We see large
+# and rather smooth structures and shapes, which are surrounded by very rough and
+# unconnected patches.
diff --git a/examples/11_plurigaussian/04_3d_pgs.py b/examples/11_plurigaussian/04_3d_pgs.py
new file mode 100644
index 00000000..3e911176
--- /dev/null
+++ b/examples/11_plurigaussian/04_3d_pgs.py
@@ -0,0 +1,80 @@
+"""
+Creating a Three Dimensional PGS
+--------------------------------
+
+Let's create a 3d PGS! This will mostly feel very familiar, but the plotting
+will be a bit more involved.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import gstools as gs
+
+dim = 3
+# no. of cells in all dimensions
+N = [40] * dim
+
+x = np.arange(N[0])
+y = np.arange(N[1])
+z = np.arange(N[2])
+
+###############################################################################
+# Because we want to create a 3d PGS, we have to generate 3 SRFs. If we are
+# interested in even higher dimensions, we could solve this code repetition
+# by using a loop...
+
+model1 = gs.Gaussian(dim=dim, var=1, len_scale=[20, 10, 15])
+srf1 = gs.SRF(model1)
+field1 = srf1.structured([x, y, z], seed=20170519)
+model2 = gs.Exponential(dim=dim, var=1, len_scale=[5, 5, 5])
+srf2 = gs.SRF(model2)
+field2 = srf2.structured([x, y, z], seed=19970221)
+model3 = gs.Gaussian(dim=dim, var=1, len_scale=[7, 12, 18])
+srf3 = gs.SRF(model3)
+field3 = srf3.structured([x, y, z], seed=20011012)
+
+###############################################################################
+# The 3d `L` field will consist of a cube which contains one category and the
+# surrounding is the second category.
+
+# size of cube
+C = [18] * dim
+
+L = np.zeros(N)
+L[
+ N[0] // 2 - C[0] // 2 : N[0] // 2 + C[0] // 2,
+ N[1] // 2 - C[1] // 2 : N[1] // 2 + C[1] // 2,
+ N[2] // 2 - C[2] // 2 : N[2] // 2 + C[2] // 2,
+] = 1
+
+###############################################################################
+# With the three SRFs and `L` ready, we can create the 3d PGS.
+pgs = gs.PGS(dim, [field1, field2, field3], L)
+
+# ###############################################################################
+# For ploting the 3d PGS, we will use [PyVista](https://pyvista.org/) which works
+# nicely together with GSTools.
+
+import pyvista as pv
+
+grid = pv.ImageData(dimensions=N)
+
+# uncomment, if you want to see `L`, which is just a cube...
+# grid.point_data['L'] = np.meshgrid(L, indexing="ij")[0]
+# grid.plot(show_edges=True)
+
+grid.point_data["PGS"] = pgs.P.reshape(-1)
+grid.contour(isosurfaces=8).plot()
+
+###############################################################################
+# .. note::
+# PyVista does not work on readthedocs, but you can try it out yourself by
+# running the example yourself. You will get an interactive version of this
+# screenshot.
+
+###############################################################################
+#
+# .. image:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/3d_pgs.png
+# :width: 400px
+# :align: center
diff --git a/examples/11_plurigaussian/README.rst b/examples/11_plurigaussian/README.rst
new file mode 100644
index 00000000..66ccafa6
--- /dev/null
+++ b/examples/11_plurigaussian/README.rst
@@ -0,0 +1,14 @@
+Plurigaussian Field Generation
+==============================
+
+Plurigaussian field simulations (PGS) are used to simulate correlated fields
+of categorical data, e.g. lithofacies, hydrofacies, soil types, or
+cementitious materials.
+PGS uses one spatial random field (SRF) per dimension, e.g. two SRFs, when
+working with two dimensional data. Furthermore, PGS needs a field, which
+describes the categorical data and its spatial relations.
+This might sound more complicated then it is, as we will see in the following
+examples.
+
+Examples
+--------
diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py
index 11e63a2b..b25785f7 100644
--- a/src/gstools/__init__.py
+++ b/src/gstools/__init__.py
@@ -45,6 +45,7 @@
.. autosummary::
SRF
CondSRF
+ PGS
Covariance Base-Class
^^^^^^^^^^^^^^^^^^^^^
@@ -162,7 +163,7 @@
TPLSimple,
TPLStable,
)
-from gstools.field import SRF, CondSRF
+from gstools.field import PGS, SRF, CondSRF
from gstools.krige import Krige
from gstools.tools import (
DEGREE_SCALE,
diff --git a/src/gstools/field/__init__.py b/src/gstools/field/__init__.py
index 9d835238..9255c174 100644
--- a/src/gstools/field/__init__.py
+++ b/src/gstools/field/__init__.py
@@ -28,10 +28,19 @@
:toctree:
Field
+
+Plurigaussian Simulation
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autosummary::
+ :toctree:
+
+ PGS
"""
from gstools.field.base import Field
from gstools.field.cond_srf import CondSRF
+from gstools.field.pgs import PGS
from gstools.field.srf import SRF
-__all__ = ["SRF", "CondSRF", "Field"]
+__all__ = ["SRF", "CondSRF", "Field", "PGS"]
diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py
index d4cb0dea..5ffb2682 100755
--- a/src/gstools/field/generator.py
+++ b/src/gstools/field/generator.py
@@ -714,7 +714,7 @@ def update(self, model=None, seed=np.nan, period=None, mode_no=None):
the seed of the random number generator.
If :any:`None`, a random seed is used. If :any:`numpy.nan`,
the actual seed will be kept. Default: :any:`numpy.nan`
- period : :class:`list` or :any:`None, optional
+ period : :class:`list` or :any:`None`, optional
The spatial periodicity of the field, is often the domain size.
mode_no : :class:`list` or :any:`None`, optional
Number of Fourier modes per dimension.
diff --git a/src/gstools/field/pgs.py b/src/gstools/field/pgs.py
new file mode 100644
index 00000000..e768acb1
--- /dev/null
+++ b/src/gstools/field/pgs.py
@@ -0,0 +1,119 @@
+"""
+GStools subpackage providing plurigaussian simulation (PGS)
+
+.. currentmodule:: gstools.field.pgs
+
+The following classes are provided
+
+.. autosummary::
+ :toctree:
+
+ PGS
+"""
+
+# pylint: disable=C0103
+import numpy as np
+
+# very clunky way of supporting both np 1.x and 2.x exceptions
+try:
+ np.AxisError = np.exceptions.AxisError
+except AttributeError:
+ ...
+
+
+class PGS:
+ """A simple class to generate plurigaussian field simulations (PGS).
+
+ See [Ricketts2023]_ for more details.
+
+ Parameters
+ ----------
+ dim : :class:`int`
+ dimension of the field
+ fields : :class:`list` or :class:`numpy.ndarray`
+ For `dim > 1` a list of spatial random fields (SRFs), with
+ `len(fields) == dim`. For `dim == 1`, the SRF can be directly given,
+ instead of a list. This class supports structured and unstructured meshes.
+ All fields must have the same shapes.
+ facies : :class:`numpy.ndarray`
+ A `dim` dimensional structured field, whose values are mapped to the PGS.
+ It does not have to have the same shape as the `fields`, as the indices are
+ automatically scaled.
+
+ References
+ ----------
+ .. [Ricketts2023] Ricketts, E.J., Freeman, B.L., Cleall, P.J. et al.
+ A Statistical Finite Element Method Integrating a Plurigaussian Random
+ Field Generator for Multi-scale Modelling of Solute Transport in
+ Concrete. Transp Porous Med 148, 95–121 (2023)
+ https://doi.org/10.1007/s11242-023-01930-8
+ """
+
+ def __init__(self, dim, fields, facies):
+ # hard to test for 1d case
+ if dim > 1:
+ if dim != len(fields):
+ raise ValueError(
+ "PGS: Mismatch between dim. and no. of fields."
+ )
+ for d in range(1, dim):
+ if not fields[0].shape == fields[d].shape:
+ raise ValueError("PGS: Not all fields have the same shape.")
+ self._dim = dim
+ self._Zs = np.array(fields)
+ self._L = np.array(facies)
+ if len(self._L.shape) != dim:
+ raise ValueError("PGS: Mismatch between dim. and facies shape.")
+ self._P = self.calc_pgs()
+
+ def calc_pgs(self):
+ """Generate the plurigaussian field.
+
+ The PGS is saved as `self.P` and is also returned.
+
+ Returns
+ -------
+ pgs : :class:`numpy.ndarray`
+ the plurigaussian field
+ """
+ try:
+ mapping = np.stack(self._Zs, axis=1)
+ except np.AxisError:
+ # if dim==1, `fields` is prob. a raw field & not a 1-tuple or
+ # equivalent
+ if self._dim == 1:
+ self._Zs = np.array([self._Zs])
+ mapping = np.stack(self._Zs, axis=1)
+ else:
+ raise
+ pos_l = []
+ try:
+ # structured grid
+ centroid = self._Zs.mean(axis=tuple(range(1, self._dim + 1)))
+ except np.AxisError:
+ # unstructured grid
+ centroid = self._Zs.mean(axis=1)
+ for d in range(self._dim):
+ l = np.floor(self._Zs[d].min()) - 1
+ h = np.ceil(self._Zs[d].max()) + 1
+ m = (h + l) / 2.0
+ dist = max(np.abs(h - m), np.abs(l - m))
+ pos_l.append(
+ np.linspace(
+ centroid[d] - dist,
+ centroid[d] + dist,
+ self._L.shape[d],
+ )
+ )
+
+ P_dig = []
+ for d in range(self._dim):
+ P_dig.append(np.digitize(mapping[:, d], pos_l[d]))
+
+ # once Py3.11 has reached its EOL, we can drop the 1-tuple :-)
+ return self._L[(*P_dig,)]
+
+ @property
+ def P(self):
+ """:class:`numpy.ndarray`: The plurigaussian field"""
+ return self._P
diff --git a/tests/test_pgs.py b/tests/test_pgs.py
new file mode 100644
index 00000000..483d04e6
--- /dev/null
+++ b/tests/test_pgs.py
@@ -0,0 +1,217 @@
+"""
+This is the unittest of the PGS class.
+"""
+
+import unittest
+
+import numpy as np
+
+import gstools as gs
+
+
+class TestPGS(unittest.TestCase):
+ def test_struct_1d(self):
+ n = 100
+ x = np.arange(n)
+ model = gs.Gaussian(dim=1, var=2, len_scale=15)
+ srf = gs.SRF(model, seed=436239)
+ field = srf.structured((x,))
+
+ m = 10
+ L = np.zeros(n)
+ L[n // 3 - m // 2 : n // 3 + m // 2] = 1
+ L[n // 3 - 2 * m : n // 3 - m // 2] = 2
+ L[4 * n // 5 - m // 2 : 4 * n // 5 + m // 2] = 3
+
+ pgs = gs.PGS(1, field, L)
+ self.assertAlmostEqual(pgs.P[n // 2], 0.0)
+ self.assertAlmostEqual(pgs.P[0], 0.0)
+ self.assertAlmostEqual(pgs.P[-1], 1.0)
+ self.assertAlmostEqual(pgs.P[-20], 3.0)
+
+ def test_struct_2d(self):
+ n1 = 100
+ n2 = 100
+ pos = [np.arange(n1), np.arange(n2)]
+
+ model1 = gs.Gaussian(dim=2, var=1, len_scale=10)
+ srf1 = gs.SRF(model1, seed=20170519)
+ field1 = srf1.structured(pos)
+
+ model2 = gs.Gaussian(dim=2, var=5, len_scale=20)
+ srf2 = gs.SRF(model2, seed=20160519)
+ field2 = srf2.structured(pos)
+
+ # create rectangle in middle of L
+ m1 = 16
+ m2 = 16
+
+ L = np.zeros((n1, n2))
+ L[
+ n1 // 2 - m1 // 2 : n1 // 2 + m1 // 2,
+ n2 // 2 - m2 // 2 : n2 // 2 + m2 // 2,
+ ] = 1
+ L[
+ n1 // 2 - m1 // 2 + m1 : n1 // 2 + m1 // 2 + m1,
+ n2 // 2 - m2 // 2 : n2 // 2 + m2 // 2,
+ ] = 2
+ L[
+ n1 // 2 - m1 // 2 + m1 : n1 // 2 + m1 // 2 + m1,
+ n2 // 2 - m2 // 2 + m2 : n2 // 2 + m2 // 2 + m2,
+ ] = 3
+ L[
+ n1 // 3 - m1 // 2 : n1 // 3 + m1 // 2,
+ n2 // 3 - m2 // 2 : n2 // 3 + m2 // 2,
+ ] = 4
+
+ pgs = gs.PGS(2, [field1, field2], L)
+
+ self.assertAlmostEqual(pgs.P[n1 // 2, n2 // 2], 2.0)
+ self.assertAlmostEqual(pgs.P[0, 0], 1.0)
+ self.assertAlmostEqual(pgs.P[-1, -1], 1.0)
+ self.assertAlmostEqual(pgs.P[0, -1], 0.0)
+ self.assertAlmostEqual(pgs.P[-1, 0], 1.0)
+
+ def test_struct_3d(self):
+ n1 = 30
+ n2 = 30
+ n3 = 30
+ pos = [np.arange(n1), np.arange(n2), np.arange(n3)]
+
+ model1 = gs.Gaussian(dim=3, var=1, len_scale=10)
+ srf1 = gs.SRF(model1, seed=20170519)
+ field1 = srf1.structured(pos)
+
+ model2 = gs.Gaussian(dim=3, var=5, len_scale=20)
+ srf2 = gs.SRF(model2, seed=20160519)
+ field2 = srf2.structured(pos)
+
+ model3 = gs.Gaussian(dim=3, var=0.1, len_scale=5)
+ srf3 = gs.SRF(model3, seed=20191219)
+ field3 = srf3.structured(pos)
+
+ # create rectangle in middle of L
+ m1 = 10
+ m2 = 10
+ m3 = 10
+
+ L = np.zeros((n1, n2, n3))
+ L[
+ n1 // 2 - m1 // 2 : n1 // 2 + m1 // 2,
+ n2 // 2 - m2 // 2 : n2 // 2 + m2 // 2,
+ n3 // 2 - m3 // 2 : n3 // 2 + m3 // 2,
+ ] = 1
+
+ pgs = gs.PGS(3, [field1, field2, field3], L)
+
+ self.assertAlmostEqual(pgs.P[n1 // 2, n2 // 2, n3 // 2], 1.0)
+ self.assertAlmostEqual(pgs.P[n1 // 3, n2 // 3, n3 // 3], 1.0)
+ self.assertAlmostEqual(
+ pgs.P[2 * n1 // 3, 2 * n2 // 3, 2 * n3 // 3], 1.0
+ )
+ self.assertAlmostEqual(pgs.P[0, 0, 0], 1.0)
+ self.assertAlmostEqual(pgs.P[-1, -1, -1], 0.0)
+ self.assertAlmostEqual(pgs.P[-1, 0, 0], 1.0)
+ self.assertAlmostEqual(pgs.P[0, -1, 0], 1.0)
+ self.assertAlmostEqual(pgs.P[0, 0, -1], 1.0)
+ self.assertAlmostEqual(pgs.P[0, -1, -1], 1.0)
+ self.assertAlmostEqual(pgs.P[-1, 0, -1], 0.0)
+ self.assertAlmostEqual(pgs.P[-1, -1, 0], 1.0)
+
+ def test_struct_4d(self):
+ n1 = 20
+ n2 = 20
+ n3 = 20
+ n4 = 20
+ pos = [np.arange(n1), np.arange(n2), np.arange(n4), np.arange(n4)]
+
+ model1 = gs.Gaussian(dim=4, var=1, len_scale=10)
+ srf1 = gs.SRF(model1, seed=20170519)
+ field1 = srf1.structured(pos)
+
+ model2 = gs.Gaussian(dim=4, var=5, len_scale=20)
+ srf2 = gs.SRF(model2, seed=20160519)
+ field2 = srf2.structured(pos)
+
+ model3 = gs.Gaussian(dim=4, var=0.1, len_scale=5)
+ srf3 = gs.SRF(model3, seed=20191219)
+ field3 = srf3.structured(pos)
+
+ model4 = gs.Exponential(dim=4, var=0.5, len_scale=12)
+ srf3 = gs.SRF(model4, seed=20191219)
+ field4 = srf3.structured(pos)
+
+ # create rectangle in middle of L
+ m1 = 5
+ m2 = 5
+ m3 = 5
+ m4 = 5
+
+ L = np.zeros((n1, n2, n3, n4))
+ L[
+ n1 // 2 - m1 // 2 : n1 // 2 + m1 // 2,
+ n2 // 2 - m2 // 2 : n2 // 2 + m2 // 2,
+ n3 // 2 - m3 // 2 : n3 // 2 + m3 // 2,
+ n4 // 2 - m4 // 2 : n4 // 2 + m4 // 2,
+ ] = 1
+
+ pgs = gs.PGS(4, [field1, field2, field3, field4], L)
+
+ self.assertAlmostEqual(pgs.P[n1 // 2, n2 // 2, n3 // 2, n4 // 2], 1.0)
+ self.assertAlmostEqual(pgs.P[0, 0, 0, 0], 0.0)
+ self.assertAlmostEqual(pgs.P[-1, -1, -1, -1], 0.0)
+
+ def test_unstruct_2d(self):
+ n1 = 10
+ n2 = 8
+ rng = np.random.RandomState(seed=438430)
+ x_unstruct = rng.randint(0, n1, size=n1 * n2)
+ y_unstruct = rng.randint(0, n2, size=n1 * n2)
+
+ pos_unstruct = [x_unstruct, y_unstruct]
+
+ model1 = gs.Gaussian(dim=2, var=1, len_scale=10)
+ srf1 = gs.SRF(model1, seed=20170519)
+ field1_unstruct = srf1.unstructured(pos_unstruct)
+
+ model2 = gs.Gaussian(dim=2, var=5, len_scale=20)
+ srf2 = gs.SRF(model2, seed=20160519)
+ field2_unstruct = srf2.unstructured(pos_unstruct)
+
+ # create rectangle in middle of L
+ m1 = 4
+ m2 = 4
+
+ L_struct = np.zeros((n1, n2))
+ L_struct[
+ n1 // 2 - m1 // 2 : n1 // 2 + m1 // 2,
+ n2 // 2 - m2 // 2 : n2 // 2 + m2 // 2,
+ ] = 1
+ L_struct[
+ n1 // 2 - m1 // 2 + m1 : n1 // 2 + m1 // 2 + m1,
+ n2 // 2 - m2 // 2 : n2 // 2 + m2 // 2,
+ ] = 2
+ L_struct[
+ n1 // 2 - m1 // 2 + m1 : n1 // 2 + m1 // 2 + m1,
+ n2 // 2 - m2 // 2 + m2 : n2 // 2 + m2 // 2 + m2,
+ ] = 3
+ L_struct[
+ n1 // 3 - m1 // 2 : n1 // 3 + m1 // 2,
+ n2 // 3 - m2 // 2 : n2 // 3 + m2 // 2,
+ ] = 4
+
+ pgs = gs.PGS(2, [field1_unstruct, field2_unstruct], L_struct)
+
+ self.assertAlmostEqual(pgs.P[0], 4.0)
+ self.assertAlmostEqual(pgs.P[-1], 1.0)
+ self.assertAlmostEqual(pgs.P[n1 * n2 // 2], 1.0)
+
+ def test_assertions(self):
+ n = 30
+ pos = [np.arange(n), np.arange(n), np.arange(n)]
+ L_2d = np.empty((n, n))
+ field1 = np.empty((n, n))
+ field2 = np.empty((n - 1, n - 1))
+ self.assertRaises(ValueError, gs.PGS, 3, [0, 1], None)
+ self.assertRaises(ValueError, gs.PGS, 3, pos, L_2d)
+ self.assertRaises(ValueError, gs.PGS, 2, [field1, field2], L_2d)