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

        Click :ref:`here <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-27

.. code-block:: python3

    import gempy as gp

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

    # new for this
    from gempy.assets import kriging

    np.random.seed(5555)








.. GENERATED FROM PYTHON SOURCE LINES 28-30

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


.. GENERATED FROM PYTHON SOURCE LINES 32-37

.. code-block:: python3

    data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'
    geo_data = gp.create_data('kriging', extent=[0, 1000, 0, 50, 0, 1000], resolution=[50, 1, 50],
                              path_o=data_path + "/data/input_data/jan_models/model1_orientations.csv",
                              path_i=data_path + "/data/input_data/jan_models/model1_surface_points.csv")





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

 Out:

 .. code-block:: none

    Active grids: ['regular']




.. GENERATED FROM PYTHON SOURCE LINES 38-40

Setting and ordering the units and series:


.. GENERATED FROM PYTHON SOURCE LINES 42-45

.. code-block:: python3

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






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <style  type="text/css" >
    #T_da883_row0_col3{
                background-color:  #015482;
            }#T_da883_row1_col3{
                background-color:  #9f0052;
            }#T_da883_row2_col3{
                background-color:  #ffbe00;
            }</style><table id="T_da883_" ><thead>    <tr>        <th class="blank level0" ></th>        <th class="col_heading level0 col0" >surface</th>        <th class="col_heading level0 col1" >series</th>        <th class="col_heading level0 col2" >order_surfaces</th>        <th class="col_heading level0 col3" >color</th>        <th class="col_heading level0 col4" >id</th>    </tr></thead><tbody>
                    <tr>
                            <th id="T_da883_level0_row0" class="row_heading level0 row0" >0</th>
                            <td id="T_da883_row0_col0" class="data row0 col0" >rock2</td>
                            <td id="T_da883_row0_col1" class="data row0 col1" >Strat_Series</td>
                            <td id="T_da883_row0_col2" class="data row0 col2" >1</td>
                            <td id="T_da883_row0_col3" class="data row0 col3" >#015482</td>
                            <td id="T_da883_row0_col4" class="data row0 col4" >1</td>
                </tr>
                <tr>
                            <th id="T_da883_level0_row1" class="row_heading level0 row1" >1</th>
                            <td id="T_da883_row1_col0" class="data row1 col0" >rock1</td>
                            <td id="T_da883_row1_col1" class="data row1 col1" >Strat_Series</td>
                            <td id="T_da883_row1_col2" class="data row1 col2" >2</td>
                            <td id="T_da883_row1_col3" class="data row1 col3" >#9f0052</td>
                            <td id="T_da883_row1_col4" class="data row1 col4" >2</td>
                </tr>
                <tr>
                            <th id="T_da883_level0_row2" class="row_heading level0 row2" >2</th>
                            <td id="T_da883_row2_col0" class="data row2 col0" >basement</td>
                            <td id="T_da883_row2_col1" class="data row2 col1" >Basement_Series</td>
                            <td id="T_da883_row2_col2" class="data row2 col2" >1</td>
                            <td id="T_da883_row2_col3" class="data row2 col3" >#ffbe00</td>
                            <td id="T_da883_row2_col4" class="data row2 col4" >3</td>
                </tr>
        </tbody></table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 46-48

Calculating the model:


.. GENERATED FROM PYTHON SOURCE LINES 50-53

.. code-block:: python3

    interp_data = gp.set_interpolator(geo_data, compile_theano=True,
                                      theano_optimizer='fast_compile')





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

 Out:

 .. code-block:: none

    Setting kriging parameters to their default values.
    Compiling theano function...
    Level of Optimization:  fast_compile
    Device:  cpu
    Precision:  float64
    Number of faults:  0
    Compilation Done!
    Kriging values: 
                        values
    range              1415.1
    $C_o$            47678.57
    drift equations    [3, 3]




.. GENERATED FROM PYTHON SOURCE LINES 54-55

no mesh computed as basically 2D model

.. GENERATED FROM PYTHON SOURCE LINES 55-57

.. code-block:: python3

    sol = gp.compute_model(geo_data, compute_mesh=False)








.. GENERATED FROM PYTHON SOURCE LINES 58-60

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


.. GENERATED FROM PYTHON SOURCE LINES 62-64

.. code-block:: python3

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




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


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

 Out:

 .. code-block:: none


    <gempy.plot.visualization_2d.Plot2D object at 0x7fcc6baa02e0>



.. GENERATED FROM PYTHON SOURCE LINES 65-75

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 77-78

conditioning data (data measured at locations)

.. GENERATED FROM PYTHON SOURCE LINES 78-81

.. code-block:: python3

    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 82-83

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

.. GENERATED FROM PYTHON SOURCE LINES 83-85

.. code-block:: python3

    domain = kriging.domain(model=sol, domain=[2], data=cond_data)








.. GENERATED FROM PYTHON SOURCE LINES 86-89

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


.. GENERATED FROM PYTHON SOURCE LINES 91-94

.. code-block:: python3

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








.. GENERATED FROM PYTHON SOURCE LINES 95-98

.. code-block:: python3

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




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





.. GENERATED FROM PYTHON SOURCE LINES 99-102

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


.. GENERATED FROM PYTHON SOURCE LINES 105-108

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 110-112

.. code-block:: python3

    solution = kriging.create_kriged_field(domain, variogram_model)








.. GENERATED FROM PYTHON SOURCE LINES 113-117

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 119-121

.. code-block:: python3

    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.0</td>
          <td>25.0</td>
          <td>410.0</td>
          <td>2.24</td>
          <td>2.13</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>430.0</td>
          <td>2.20</td>
          <td>2.01</td>
        </tr>
        <tr>
          <th>2</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>450.0</td>
          <td>2.15</td>
          <td>1.90</td>
        </tr>
        <tr>
          <th>3</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>470.0</td>
          <td>2.10</td>
          <td>1.81</td>
        </tr>
        <tr>
          <th>4</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>490.0</td>
          <td>2.07</td>
          <td>1.77</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 122-125

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


.. GENERATED FROM PYTHON SOURCE LINES 127-132

.. code-block:: python3

    solution.plot_results(geo_data=geo_data, prop='val', contour=False,
                          direction='y', cell_number=0, alpha=0.7,
                          show_data=False, legend=True)
    plt.show()




.. image:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_003.png
    :alt: ch3 1 kriging interpolation and simulation
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none

    /WorkSSD/PythonProjects/gempy/gempy/assets/kriging.py:265: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
      cmap.set_bad(color='w', alpha=alpha) #define color and alpha for nan values




.. GENERATED FROM PYTHON SOURCE LINES 133-137

.. code-block:: python3

    solution.plot_results(geo_data=geo_data, prop='both', contour=False,
                          direction='y', cell_number=0, alpha=0,
                          interpolation='bilinear', show_data=False)
    plt.show()



.. image:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_004.png
    :alt: Estimated value, Variance
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none

    /WorkSSD/PythonProjects/gempy/gempy/assets/kriging.py:265: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
      cmap.set_bad(color='w', alpha=alpha) #define color and alpha for nan values




.. GENERATED FROM PYTHON SOURCE LINES 138-145

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 147-149

.. code-block:: python3

    solution_sim = kriging.create_gaussian_field(domain, variogram_model)








.. GENERATED FROM PYTHON SOURCE LINES 150-152

.. code-block:: python3

    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>16</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>410.0</td>
          <td>2.61</td>
          <td>1.48</td>
        </tr>
        <tr>
          <th>157</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>430.0</td>
          <td>2.38</td>
          <td>0.28</td>
        </tr>
        <tr>
          <th>276</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>450.0</td>
          <td>3.98</td>
          <td>0.24</td>
        </tr>
        <tr>
          <th>25</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>470.0</td>
          <td>2.54</td>
          <td>0.53</td>
        </tr>
        <tr>
          <th>267</th>
          <td>10.0</td>
          <td>25.0</td>
          <td>490.0</td>
          <td>0.76</td>
          <td>0.28</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 153-155

.. code-block:: python3

    solution_sim.results_df['estimated value']





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

 Out:

 .. code-block:: none


    16     2.61
    157    2.38
    276    3.98
    25     2.54
    267    0.76
           ... 
    55     3.76
    165    3.32
    427    2.01
    303    1.96
    187    2.04
    Name: estimated value, Length: 500, dtype: float64



.. GENERATED FROM PYTHON SOURCE LINES 156-157

sphinx_gallery_thumbnail_number = 3

.. GENERATED FROM PYTHON SOURCE LINES 157-161

.. code-block:: python3

    solution_sim.plot_results(geo_data=geo_data, prop='val', contour=False, direction='y', cell_number=0, alpha=0.7,
                              show_data=True, legend=True)
    plt.show()




.. image:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_005.png
    :alt: ch3 1 kriging interpolation and simulation
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none

    /WorkSSD/PythonProjects/gempy/gempy/assets/kriging.py:265: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
      cmap.set_bad(color='w', alpha=alpha) #define color and alpha for nan values




.. GENERATED FROM PYTHON SOURCE LINES 162-164

.. code-block:: python3

    solution_sim.plot_results(geo_data=geo_data, prop='both', contour=False, direction='y', cell_number=0, alpha=0,
                              interpolation='bilinear', show_data=False)
    plt.show()


.. image:: /tutorials/ch3-Interpolations/images/sphx_glr_ch3_1_kriging_interpolation_and_simulation_006.png
    :alt: Estimated value, Variance
    :class: sphx-glr-single-img


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

 Out:

 .. code-block:: none

    /WorkSSD/PythonProjects/gempy/gempy/assets/kriging.py:265: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
      cmap.set_bad(color='w', alpha=alpha) #define color and alpha for nan values





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

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


.. _sphx_glr_download_tutorials_ch3-Interpolations_ch3_1_kriging_interpolation_and_simulation.py:


.. only :: html

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



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



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


.. only:: html

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

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