.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/11_plurigaussian/05_conditioned.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_11_plurigaussian_05_conditioned.py: Creating conditioned PGS ------------------------ In case we have knowledge about values of the PGS in certain areas, we should incorporate that knowledge. We can do this by using conditional fields, which are created by combining kriged fields with SRFs. This can be done quite easily with GSTools. For more details, have a look at the kriging and conditional field examples. In this example, we assume that we know the categorical values of the PGS field to be 2 in the lower left 20 by 20 grid cells. Warning: Using PGS for conditioning fields is still a beta feature. .. GENERATED FROM PYTHON SOURCE LINES 16-30 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import gstools as gs dim = 2 # no. of cells in both dimensions N = [100, 80] # grid x = np.arange(N[0]) y = np.arange(N[1]) .. GENERATED FROM PYTHON SOURCE LINES 31-41 Now we want to read the known data in order to condition the PGS on them. Normally, we would probably get the data in a different format, but in order to keep the dependencies for this example to a minimum, we will simpy use a numpy format. After reading in the data, we will have a quick look at how the data looks like. We'll see the first few values, which are all 1. This value is not very important. However, it should be in the range of SRF values, to not mess with the kriging. The known value of the PGS, namely 2, will be set for the lithotypes field, which will map it the the conditioned SRF values of 1 at given positions. .. GENERATED FROM PYTHON SOURCE LINES 41-48 .. code-block:: Python cond_data = np.load("conditional_values.npz") cond_pos = cond_data["cond_pos"] cond_val = cond_data["cond_val"] print(f"first 5 conditional positions:\n{cond_pos[:, :5]}") print(f"first 5 conditional values:\n{cond_val[:5]}") .. rst-class:: sphx-glr-script-out .. code-block:: none first 5 conditional positions: [[0 1 2 3 4] [0 0 0 0 0]] first 5 conditional values: [1 1 1 1 1] .. GENERATED FROM PYTHON SOURCE LINES 49-53 With the conditional values ready, we can now set up the covariance model for the kriging. This knowledge has to normally be inferred, but here we just assume that we know the convariance structure of the underlying field. For better visualization, we use `Simple` kriging with a mean value of 0. .. GENERATED FROM PYTHON SOURCE LINES 53-59 .. code-block:: Python model = gs.Gaussian(dim=dim, var=1, len_scale=[10, 5], angles=np.pi / 8) krige = gs.krige.Simple(model, cond_pos=cond_pos, cond_val=cond_val, mean=0) cond_srf = gs.CondSRF(krige) cond_srf.set_pos([x, y], "structured") .. GENERATED FROM PYTHON SOURCE LINES 60-66 Now that the conditioned field class is set up, we can generate SRFs conditioned on our previous knowledge. We'll do that for the two SRFs needed for the PGS, and then we will also set up the PGS generator. Next, we'll use a little helper method, which can transform the coordinates from the SRFs to the lithotypes field. This helps us set up the area around the conditioned value `cond_val`. .. GENERATED FROM PYTHON SOURCE LINES 66-91 .. code-block:: Python field1 = cond_srf(seed=484739) field2 = cond_srf(seed=45755894) pgs = gs.PGS(dim, [field1, field2]) M = [100, 80] # size of the rectangle rect = [40, 32] lithotypes = np.zeros(M) # calculate grid axes of the lithotypes field pos_lith = pgs.calc_lithotype_axes(lithotypes.shape) # transform conditioned SRF value to lithotypes index pos_lith_ind = pgs.transform_coords( lithotypes.shape, [cond_val[0], cond_val[0]] ) # conditioned category of 2 around the conditioned values' positions lithotypes[ pos_lith_ind[0] - 5 : pos_lith_ind[0] + 5, pos_lith_ind[1] - 5 : pos_lith_ind[1] + 5, ] = 2 .. GENERATED FROM PYTHON SOURCE LINES 92-93 With the two SRFs and the lithotypes ready, we can create the actual PGS. .. GENERATED FROM PYTHON SOURCE LINES 93-96 .. code-block:: Python P = pgs(lithotypes) .. GENERATED FROM PYTHON SOURCE LINES 97-104 Finally, we can plot the PGS, but we will also show the lithotypes and the two original Gaussian SRFs. We will set the colours of the SRF correlation scatter plot to be the sum of their respective position tuples (x+y), to get a feeling for which point corresponds to which position. The more blue the points, the smaller the sum is. We can nicely see that many blue points gather in the highlighted rectangle of the lithotypes where the categorical value of 2 is set. .. GENERATED FROM PYTHON SOURCE LINES 104-123 .. code-block:: Python fig, axs = plt.subplots(2, 2) axs[0, 0].imshow(field1, cmap="copper", vmin=-3, vmax=2, origin="lower") axs[0, 1].imshow(field2, cmap="copper", vmin=-3, vmax=2, origin="lower") axs[1, 0].scatter( field1.flatten(), field2.flatten(), s=0.1, c=(x.reshape((len(x), 1)) + y.reshape((1, len(y))).flatten()), ) axs[1, 0].pcolormesh( pos_lith[0], pos_lith[1], lithotypes.T, alpha=0.3, cmap="copper" ) axs[1, 1].imshow(P, cmap="copper", origin="lower") plt.tight_layout() plt.show() .. image-sg:: /examples/11_plurigaussian/images/sphx_glr_05_conditioned_001.png :alt: 05 conditioned :srcset: /examples/11_plurigaussian/images/sphx_glr_05_conditioned_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 124-126 With all this set up, we can easily create an ensemble of PGS, which conform to the conditional values .. GENERATED FROM PYTHON SOURCE LINES 126-149 .. code-block:: Python seed = gs.random.MasterRNG(20170519) ens_no = 9 fields1 = [] fields2 = [] Ps = [] for i in range(ens_no): fields1.append(cond_srf(seed=seed())) fields2.append(cond_srf(seed=seed())) pgs = gs.PGS(dim, [fields1[-1], fields2[-1]]) Ps.append(pgs(lithotypes)) fig, axs = plt.subplots(3, 3) cnt = 0 for i in range(int(np.sqrt(ens_no))): for j in range(int(np.sqrt(ens_no))): axs[i, j].imshow(Ps[cnt], cmap="copper", origin="lower") cnt += 1 plt.tight_layout() plt.show() .. image-sg:: /examples/11_plurigaussian/images/sphx_glr_05_conditioned_002.png :alt: 05 conditioned :srcset: /examples/11_plurigaussian/images/sphx_glr_05_conditioned_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.584 seconds) .. _sphx_glr_download_examples_11_plurigaussian_05_conditioned.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 05_conditioned.ipynb <05_conditioned.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 05_conditioned.py <05_conditioned.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 05_conditioned.zip <05_conditioned.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_