.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/11_plurigaussian/01_pgs.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_11_plurigaussian_01_pgs.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_11_plurigaussian_01_pgs.py:


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.

.. GENERATED FROM PYTHON SOURCE LINES 9-22

.. 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]

    x = np.arange(N[0])
    y = np.arange(N[1])








.. GENERATED FROM PYTHON SOURCE LINES 23-26

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.

.. GENERATED FROM PYTHON SOURCE LINES 26-35

.. code-block:: Python


    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








.. GENERATED FROM PYTHON SOURCE LINES 36-40

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.

.. GENERATED FROM PYTHON SOURCE LINES 40-43

.. code-block:: Python


    plt.scatter(field1.flatten(), field2.flatten(), s=0.1)




.. image-sg:: /examples/11_plurigaussian/images/sphx_glr_01_pgs_001.png
   :alt: 01 pgs
   :srcset: /examples/11_plurigaussian/images/sphx_glr_01_pgs_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 44-49

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.

.. GENERATED FROM PYTHON SOURCE LINES 49-65

.. code-block:: Python


    # 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








.. GENERATED FROM PYTHON SOURCE LINES 66-70

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.

.. GENERATED FROM PYTHON SOURCE LINES 70-76

.. code-block:: Python


    pgs = gs.PGS(dim, [field1, field2])
    P = pgs(lithotypes)

    x_lith, y_lith = pgs.calc_lithotype_axes(lithotypes.shape)








.. GENERATED FROM PYTHON SOURCE LINES 77-79

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.

.. GENERATED FROM PYTHON SOURCE LINES 79-88

.. code-block:: Python


    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")




.. image-sg:: /examples/11_plurigaussian/images/sphx_glr_01_pgs_002.png
   :alt: 01 pgs
   :srcset: /examples/11_plurigaussian/images/sphx_glr_01_pgs_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 89-94

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.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.682 seconds)


.. _sphx_glr_download_examples_11_plurigaussian_01_pgs.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: 01_pgs.ipynb <01_pgs.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 01_pgs.py <01_pgs.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: 01_pgs.zip <01_pgs.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_