|
| 1 | +""" |
| 2 | +A First and Simple Example |
| 3 | +-------------------------- |
| 4 | +
|
| 5 | +As a first example, we will create a two dimensional plurigaussian field |
| 6 | +(PGS). Thus, we need two spatial random fields(SRF) and on top of that, we |
| 7 | +need a field describing the categorical data and its spatial relation. |
| 8 | +We will start off by creating the two SRFs with a Gaussian variogram, which |
| 9 | +makes the fields nice and smooth. But before that, we will import all |
| 10 | +necessary libraries and define a few variables, like the number of grid |
| 11 | +cells in each dimension. |
| 12 | +""" |
| 13 | + |
| 14 | +import matplotlib.pyplot as plt |
| 15 | +import numpy as np |
| 16 | + |
| 17 | +import gstools as gs |
| 18 | + |
| 19 | +dim = 2 |
| 20 | +# no. of cells in both dimensions |
| 21 | +N = [180, 140] |
| 22 | + |
| 23 | +x = np.arange(N[0]) |
| 24 | +y = np.arange(N[1]) |
| 25 | + |
| 26 | +############################################################################### |
| 27 | +# In this first example we will use the same geostatistical parameters for |
| 28 | +# both fields for simplicity. Thus, we can use the same SRF instance for the |
| 29 | +# two fields. |
| 30 | + |
| 31 | +model = gs.Gaussian(dim=dim, var=1, len_scale=10) |
| 32 | +srf = gs.SRF(model) |
| 33 | +field1 = srf.structured([x, y], seed=20170519) |
| 34 | +field2 = srf.structured([x, y], seed=19970221) |
| 35 | + |
| 36 | +############################################################################### |
| 37 | +# Now, we will create the lithotypes field describing the categorical data. For |
| 38 | +# now, we will only have two categories and we will address them by the |
| 39 | +# integers 0 and 1. We start off by creating a matrix of 0s from which we will |
| 40 | +# select a rectangle and fill that with 1s. This field does not have to match |
| 41 | +# the shape of the SRFs. |
| 42 | + |
| 43 | +centroid = [200, 160] |
| 44 | + |
| 45 | +# size of the rectangle |
| 46 | +rect = [40, 32] |
| 47 | + |
| 48 | +lithotypes = np.zeros(centroid) |
| 49 | +lithotypes[ |
| 50 | + centroid[0] // 2 - rect[0] // 2 : centroid[0] // 2 + rect[0] // 2, |
| 51 | + centroid[1] // 2 - rect[1] // 2 : centroid[1] // 2 + rect[1] // 2, |
| 52 | +] = 1 |
| 53 | + |
| 54 | +############################################################################### |
| 55 | +# With the two SRFs and the L-field ready, we can create our first PGS. First, we |
| 56 | +# set up an instance of the PGS class and then we are ready to calculate the |
| 57 | +# field by calling the instance and handing over the L-field. |
| 58 | + |
| 59 | +pgs = gs.PGS(dim, [field1, field2]) |
| 60 | +P = pgs(lithotypes) |
| 61 | + |
| 62 | +############################################################################### |
| 63 | +# Finally, we can plot the PGS, but we will also show the L-field and the two |
| 64 | +# original Gaussian fields. |
| 65 | + |
| 66 | +fig, axs = plt.subplots(2, 2) |
| 67 | + |
| 68 | +axs[0, 0].imshow(field1, cmap="copper", origin="lower") |
| 69 | +axs[0, 1].imshow(field2, cmap="copper", origin="lower") |
| 70 | +axs[1, 0].imshow(lithotypes, cmap="copper", origin="lower") |
| 71 | +axs[1, 1].imshow(P, cmap="copper", origin="lower") |
| 72 | + |
| 73 | +# For more information on Plurigaussian fields and how they naturally extend |
| 74 | +# truncated Gaussian fields, we can recommend the book |
| 75 | +# [Plurigaussian Simulations in Geosciences](https://doi.org/10.1007/978-3-642-19607-2) |
0 commit comments