.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "tutorials/ch3-Interpolations/ch3_1_kriging_interpolation_and_simulation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_tutorials_ch3-Interpolations_ch3_1_kriging_interpolation_and_simulation.py>`
        to download the full example code

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

.. _sphx_glr_tutorials_ch3-Interpolations_ch3_1_kriging_interpolation_and_simulation.py:


3.1: Simple example of kriging in gempy
=======================================

.. GENERATED FROM PYTHON SOURCE LINES 8-13

In this notebook it will be shown how to create a kriged or simulated
field in a simple geological model in gempy. We start by creating a
simple model with three horizontally layered units, as shown in the
gempy examples.


.. GENERATED FROM PYTHON SOURCE LINES 15-16

Importing GemPy

.. GENERATED FROM PYTHON SOURCE LINES 16-29

.. code-block:: Python

    import gempy as gp
    import gempy_viewer as gpv

    # Importing auxiliary libraries
    import numpy as np
    import matplotlib.pyplot as plt
    import os

    # new for this
    from gempy_plugins.kriging import kriging

    np.random.seed(5555)








.. GENERATED FROM PYTHON SOURCE LINES 30-32

Creating the model by importing the input data and displaying it:


.. GENERATED FROM PYTHON SOURCE LINES 34-47

.. code-block:: Python

    data_path = os.path.abspath('../../')

    geo_data: gp.data.GeoModel = gp.create_geomodel(
        project_name='kriging',
        extent=[0, 1000, 0, 50, 0, 1000],
        resolution=[50, 10, 50],
        refinement=1,
        importer_helper=gp.data.ImporterHelper(
            path_to_orientations=data_path + "/data/input_data/jan_models/model1_orientations.csv",
            path_to_surface_points=data_path + "/data/input_data/jan_models/model1_surface_points.csv",
        )
    )








.. GENERATED FROM PYTHON SOURCE LINES 48-50

Setting and ordering the units and series:


.. GENERATED FROM PYTHON SOURCE LINES 52-60

.. code-block:: Python

    gp.map_stack_to_surfaces(
        gempy_model=geo_data,
        mapping_object={
            "Strat_Series": ('rock2', 'rock1'),
            "Basement_Series": ('basement')
        }
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Could not find element 'basement' in any group.


.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">

            <table>
              <tr><td>Structural Groups:</td><td>
        <table style="border-left:1.2px solid black;>
          <tr><th colspan="2"><b>StructuralGroup:</b></th></tr>
          <tr><td>Name:</td><td>Strat_Series</td></tr>
          <tr><td>Structural Relation:</td><td>StackRelationType.ERODE</td></tr>
          <tr><td>Elements:</td><td>
        <table width="50%" style="border-left:15px solid #9f0052;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>rock2</td></tr>
        </table>
            <br>
        <table width="50%" style="border-left:15px solid #015482;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>rock1</td></tr>
        </table>
            </td></tr>
        </table>
            </td></tr>
              <tr><td>Fault Relations:</td><td><table style="border-collapse: collapse; table-layout: fixed;"><th></th><th style="transform: rotate(-35deg); height:150px; vertical-align: bottom; text-align: center;">Strat_Seri...</th><tr><th>Strat_Series</th><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td></tr></table></td></tr>
              <tr><td></td><td>
            <table>
              <tr>
                <td><div style="display: inline-block; background-color: #527682; width: 20px; height: 20px; border: 1px solid black;"></div> True</td>
                <td><div style="display: inline-block; background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></div> False</td>
              </tr>
            </table>
            </td></tr>
            </table>
        
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 61-63

Calculating the model:


.. GENERATED FROM PYTHON SOURCE LINES 65-66

no mesh computed as basically 2D model

.. GENERATED FROM PYTHON SOURCE LINES 66-68

.. code-block:: Python

    sol = gp.compute_model(geo_data)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Setting Backend To: AvailableBackends.numpy
    /home/leguark/gempy_engine/gempy_engine/modules/dual_contouring/dual_contouring_interface.py:38: RuntimeWarning: divide by zero encountered in divide
      weight_x = ((scalar_at_sp - scalar_8[:, :, 4:]) / scalar_dx).reshape(-1, 4, 1)
    /home/leguark/gempy_engine/gempy_engine/modules/dual_contouring/dual_contouring_interface.py:39: RuntimeWarning: divide by zero encountered in divide
      weight_y = ((scalar_at_sp - scalar_8[:, :, [2, 3, 6, 7]]) / scalar_d_y).reshape(-1, 4, 1)
    /home/leguark/gempy_engine/gempy_engine/modules/dual_contouring/dual_contouring_interface.py:48: RuntimeWarning: invalid value encountered in multiply
      intersect_dx = d_x[:, :, :] * weight_x[:, :, :]
    /home/leguark/gempy_engine/gempy_engine/modules/dual_contouring/dual_contouring_interface.py:49: RuntimeWarning: invalid value encountered in multiply
      intersect_dy = d_y[:, :, :] * weight_y[:, :, :]




.. GENERATED FROM PYTHON SOURCE LINES 69-71

So here is the very simple, basically 2D model that we created:


.. GENERATED FROM PYTHON SOURCE LINES 73-75

.. code-block:: Python

    gpv.plot_2d(geo_data, cell_number=0, show_data=False)




.. image-sg:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_001.png
   :alt: Cell Number: 0 Direction: y
   :srcset: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f490cb745e0>



.. GENERATED FROM PYTHON SOURCE LINES 76-86

1) Creating domain
------------------

Let us assume we have a couple of measurements in a domain of interest
within our model. In our case the unit of interest is the central rock
layer (rock1). In the kriging module we can define the domain by
handing over a number of surfaces by id - in this case the id of rock1
is 2. In addition we define four input data points in cond_data, each
defined by x,y,z coordinate and a measurement value.


.. GENERATED FROM PYTHON SOURCE LINES 88-89

conditioning data (data measured at locations)

.. GENERATED FROM PYTHON SOURCE LINES 89-92

.. code-block:: Python

    cond_data = np.array([[100, .5, 500, 2], [900, .5, 500, 1],
                          [500, .5, 550, 1], [300, .5, 400, 5]])








.. GENERATED FROM PYTHON SOURCE LINES 93-94

creating a domain object from the gempy solution, a defined domain conditioning data

.. GENERATED FROM PYTHON SOURCE LINES 94-101

.. code-block:: Python

    domain = kriging.Domain(
        model_solutions=sol,
        transform=geo_data.transform,
        domain=[2],
        data=cond_data
    )








.. GENERATED FROM PYTHON SOURCE LINES 102-105

2) Creating a variogram model
-----------------------------


.. GENERATED FROM PYTHON SOURCE LINES 107-113

.. code-block:: Python

    variogram_model = kriging.VariogramModel(
        theoretical_model='exponential',
        range_=200,
        sill=np.var(cond_data[:, 3])
    )








.. GENERATED FROM PYTHON SOURCE LINES 114-117

.. code-block:: Python

    variogram_model.plot(type_='both', show_parameters=True)
    plt.show()




.. image-sg:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_002.png
   :alt: Models of spatial correlation
   :srcset: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 118-121

3) Kriging interpolation
------------------------


.. GENERATED FROM PYTHON SOURCE LINES 124-127

In the following we define an object called kriging_model and set all
input parameters. Finally we generate the kriged field.


.. GENERATED FROM PYTHON SOURCE LINES 129-131

.. code-block:: Python

    kriging_solution = kriging.create_kriged_field(domain, variogram_model)








.. GENERATED FROM PYTHON SOURCE LINES 132-136

The result of our calculation is saved in the following dataframe,
containing an estimated value and the kriging variance for each point in
the grid:


.. GENERATED FROM PYTHON SOURCE LINES 138-141

.. code-block:: Python

    kriging_solution.results_df.head()







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>X</th>
          <th>Y</th>
          <th>Z</th>
          <th>estimated value</th>
          <th>estimation variance</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>410.0016</td>
          <td>2.241542</td>
          <td>2.103698</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>430.0016</td>
          <td>2.193555</td>
          <td>1.979810</td>
        </tr>
        <tr>
          <th>2</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>450.0016</td>
          <td>2.144345</td>
          <td>1.865440</td>
        </tr>
        <tr>
          <th>3</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>470.0016</td>
          <td>2.097883</td>
          <td>1.774970</td>
        </tr>
        <tr>
          <th>4</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>490.0016</td>
          <td>2.058535</td>
          <td>1.724997</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 142-145

It is also possible to plot the results in cross section similar to the
way gempy models are plotted.


.. GENERATED FROM PYTHON SOURCE LINES 148-165

.. code-block:: Python

    from gempy_viewer.modules.plot_2d.visualization_2d import Plot2D

    plot_2d: Plot2D = gpv.plot_2d(
        model=geo_data,
        cell_number=0,
        show_data=False,
        show=False,
        kwargs_lithology={ 'alpha': 0.5 }
    )
    kriging.plot_kriging_results(
        geo_data=geo_data,
        kriging_solution=kriging_solution,
        plot_2d=plot_2d,
        title='Kriging interpolation: Estimated values',
        result_column=['estimated value']
    )




.. image-sg:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_003.png
   :alt: Kriging interpolation: Estimated values
   :srcset: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 166-182

.. code-block:: Python

    plot_2d_both = gpv.plot_2d(
        model=geo_data,
        cell_number=[0, 0],
        show_data=False,
        show=False,
        kwargs_lithology={ 'alpha': 0.5 }
    )

    kriging.plot_kriging_results(
        geo_data=geo_data,
        kriging_solution=kriging_solution,
        plot_2d=plot_2d_both,
        title='Kriging interpolation: Estimated values',
        result_column=['estimated value', 'estimation variance']
    )




.. image-sg:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_004.png
   :alt: Kriging interpolation: Estimated values, Kriging interpolation: Estimated values
   :srcset: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 183-190

4) Simulated field
------------------

Based on the same objects (domain and varigoram model) also a simulated
field (stationary Gaussian Field) can be generated. A Sequential
Gaussian Simulation approach is applied in this module:


.. GENERATED FROM PYTHON SOURCE LINES 192-194

.. code-block:: Python

    solution_sim = kriging.create_gaussian_field(domain, variogram_model)








.. GENERATED FROM PYTHON SOURCE LINES 195-197

.. code-block:: Python

    solution_sim.results_df.head()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>X</th>
          <th>Y</th>
          <th>Z</th>
          <th>estimated value</th>
          <th>estimation variance</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>1774</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>410.0016</td>
          <td>4.420467</td>
          <td>0.122209</td>
        </tr>
        <tr>
          <th>4701</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>430.0016</td>
          <td>4.893653</td>
          <td>0.108040</td>
        </tr>
        <tr>
          <th>3463</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>450.0016</td>
          <td>4.874381</td>
          <td>0.177325</td>
        </tr>
        <tr>
          <th>3227</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>470.0016</td>
          <td>4.983143</td>
          <td>0.220953</td>
        </tr>
        <tr>
          <th>4528</th>
          <td>10.0016</td>
          <td>2.5016</td>
          <td>490.0016</td>
          <td>3.747620</td>
          <td>0.108041</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 198-200

.. code-block:: Python

    solution_sim.results_df['estimated value']





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    1774    4.420467
    4701    4.893653
    3463    4.874381
    3227    4.983143
    4528    3.747620
              ...   
    1712    3.540901
    390     3.863402
    1039    4.134937
    3379    3.144726
    4090    2.446597
    Name: estimated value, Length: 5000, dtype: float64



.. GENERATED FROM PYTHON SOURCE LINES 201-216

.. code-block:: Python

    plot_2d: Plot2D = gpv.plot_2d(
        model=geo_data,
        cell_number=0,
        show_data=False,
        show=False,
        kwargs_lithology={ 'alpha': 0.5 }
    )
    kriging.plot_kriging_results(
        geo_data=geo_data,
        kriging_solution=solution_sim,
        plot_2d=plot_2d,
        title='Kriging interpolation: Estimated values',
        result_column=['estimated value']
    )




.. image-sg:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_005.png
   :alt: Kriging interpolation: Estimated values
   :srcset: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 217-233

.. code-block:: Python

    plot_2d_both = gpv.plot_2d(
        model=geo_data,
        cell_number=[0, 0],
        show_data=False,
        show=False,
        kwargs_lithology={ 'alpha': 0.5 }
    )

    kriging.plot_kriging_results(
        geo_data=geo_data,
        kriging_solution=solution_sim,
        plot_2d=plot_2d_both,
        title='Kriging interpolation: Estimated values',
        result_column=['estimated value', 'estimation variance']
    )
    # sphinx_gallery_thumbnail_number = 3



.. image-sg:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_006.png
   :alt: Kriging interpolation: Estimated values, Kriging interpolation: Estimated values
   :srcset: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_006.png
   :class: sphx-glr-single-img






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

   **Total running time of the script:** (15 minutes 23.937 seconds)


.. _sphx_glr_download_tutorials_ch3-Interpolations_ch3_1_kriging_interpolation_and_simulation.py:

.. only:: html

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

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

      :download:`Download Jupyter notebook: ch3_1_kriging_interpolation_and_simulation.ipynb <ch3_1_kriging_interpolation_and_simulation.ipynb>`

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

      :download:`Download Python source code: ch3_1_kriging_interpolation_and_simulation.py <ch3_1_kriging_interpolation_and_simulation.py>`


.. only:: html

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

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