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

.. only:: html

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

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

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

.. _sphx_glr_tutorials_ch2-Geophysics_ch2_1_gravity.py:


2.1 Forward Gravity: Simple example
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 8-9

Importing gempy

.. GENERATED FROM PYTHON SOURCE LINES 9-22

.. code-block:: Python

    import gempy as gp
    import gempy_viewer as gpv

    # Aux imports
    import numpy as np
    import pandas as pd
    import os
    import matplotlib.pyplot as plt


    np.random.seed(1515)
    pd.set_option('display.precision', 2)








.. GENERATED FROM PYTHON SOURCE LINES 23-39

.. code-block:: Python

    data_path = os.path.abspath('../../data/input_data/tut_SandStone')

    # Importing the data from csv

    geo_model: gp.data.GeoModel = gp.create_geomodel(
        project_name='Greenstone',
        extent=[696000, 747000, 6863000, 6930000, -20000, 200],  # * Here we define the extent of the model
        refinement=5,  # * Here we define the number of octree levels. If octree levels are defined, the resolution is ignored.
        importer_helper=gp.data.ImporterHelper(
            path_to_orientations=data_path + "/SandStone_Foliations.csv",
            path_to_surface_points=data_path + "/SandStone_Points.csv",
            hash_surface_points=None,
            hash_orientations=None
        )
    )








.. GENERATED FROM PYTHON SOURCE LINES 40-49

.. code-block:: Python

    gp.map_stack_to_surfaces(
        gempy_model=geo_model,
        mapping_object={
            "EarlyGranite_Series": 'EarlyGranite',
            "BIF_Series": ('SimpleMafic2', 'SimpleBIF'),
            "SimpleMafic_Series": 'SimpleMafic1', 'Basement': '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>EarlyGranite_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 #015482;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>EarlyGranite</td></tr>
        </table>
            </td></tr>
        </table>
            <br>
        <table style="border-left:1.2px solid black;>
          <tr><th colspan="2"><b>StructuralGroup:</b></th></tr>
          <tr><td>Name:</td><td>BIF_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 #728f02;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>SimpleMafic2</td></tr>
        </table>
            <br>
        <table width="50%" style="border-left:15px solid #9f0052;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>SimpleBIF</td></tr>
        </table>
            </td></tr>
        </table>
            <br>
        <table style="border-left:1.2px solid black;>
          <tr><th colspan="2"><b>StructuralGroup:</b></th></tr>
          <tr><td>Name:</td><td>SimpleMafic_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 #ffbe00;">
          <tr><th colspan="2"><b>StructuralElement:</b></th></tr>
          <tr><td>Name:</td><td>SimpleMafic1</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;">EarlyGrani...</th><th style="transform: rotate(-35deg); height:150px; vertical-align: bottom; text-align: center;">BIF_Series</th><th style="transform: rotate(-35deg); height:150px; vertical-align: bottom; text-align: center;">SimpleMafi...</th><tr><th>EarlyGranite_Series</th><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td></tr><tr><th>BIF_Series</th><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td></tr><tr><th>SimpleMafic_Series</th><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td><td style="background-color: #FFB6C1; width: 20px; height: 20px; border: 1px solid black;"></td><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 50-52

.. code-block:: Python

    gpv.plot_2d(geo_model)




.. image-sg:: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_001.png
   :alt: Cell Number: mid Direction: y
   :srcset: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_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 0x7f919f1713c0>



.. GENERATED FROM PYTHON SOURCE LINES 55-58

Creating grid
~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 60-63

First we need to define the location of the devices. For this example we
can make a map:


.. GENERATED FROM PYTHON SOURCE LINES 65-73

.. code-block:: Python

    grav_res = 20
    X = np.linspace(7.050000e+05, 747000, grav_res)
    Y = np.linspace(6863000, 6925000, grav_res)
    Z = 300
    xyz = np.meshgrid(X, Y, Z)
    xy_ravel = np.vstack(list(map(np.ravel, xyz))).T
    xy_ravel





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

 .. code-block:: none


    array([[7.05000000e+05, 6.86300000e+06, 3.00000000e+02],
           [7.07210526e+05, 6.86300000e+06, 3.00000000e+02],
           [7.09421053e+05, 6.86300000e+06, 3.00000000e+02],
           ...,
           [7.42578947e+05, 6.92500000e+06, 3.00000000e+02],
           [7.44789474e+05, 6.92500000e+06, 3.00000000e+02],
           [7.47000000e+05, 6.92500000e+06, 3.00000000e+02]])



.. GENERATED FROM PYTHON SOURCE LINES 74-76

We can see the location of the devices relative to the model data:


.. GENERATED FROM PYTHON SOURCE LINES 78-83

.. code-block:: Python


    gpv.plot_2d(geo_model, direction='z', show=False)
    plt.scatter(xy_ravel[:, 0], xy_ravel[:, 1], s=1)
    plt.show()




.. image-sg:: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_002.png
   :alt: Cell Number: mid Direction: z
   :srcset: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 84-87

Now we need to create the grid centered on the devices (see:
https://github.com/cgre-aachen/gempy/blob/master/notebooks/tutorials/ch1-3-Grids.ipynb)


.. GENERATED FROM PYTHON SOURCE LINES 89-90

geo_model.set_centered_grid(xy_ravel, resolution=[10, 10, 15], radius=5000)

.. GENERATED FROM PYTHON SOURCE LINES 90-98

.. code-block:: Python


    gp.set_centered_grid(
        grid=geo_model.grid,
        centers=xy_ravel,
        resolution=np.array([10, 10, 15]),
        radius=np.array([5000, 5000, 5000])
    )





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

 .. code-block:: none

    Active grids: ['centered']

    CenteredGrid(centers=array([[7.05000000e+05, 6.86300000e+06, 3.00000000e+02],
           [7.07210526e+05, 6.86300000e+06, 3.00000000e+02],
           [7.09421053e+05, 6.86300000e+06, 3.00000000e+02],
           ...,
           [7.42578947e+05, 6.92500000e+06, 3.00000000e+02],
           [7.44789474e+05, 6.92500000e+06, 3.00000000e+02],
           [7.47000000e+05, 6.92500000e+06, 3.00000000e+02]]), resolution=array([10, 10, 15]), radius=array([5000, 5000, 5000]), kernel_grid_centers=array([[-5000.        , -5000.        ,  -300.        ],
           [-5000.        , -5000.        ,  -360.        ],
           [-5000.        , -5000.        ,  -383.36972966],
           ...,
           [ 5000.        ,  5000.        , -3407.68480754],
           [ 5000.        ,  5000.        , -4618.11403801],
           [ 5000.        ,  5000.        , -6300.        ]]), left_voxel_edges=array([[1709.43058496, 1709.43058496,  -30.        ],
           [1709.43058496, 1709.43058496,  -30.        ],
           [1709.43058496, 1709.43058496,  -11.68486483],
           ...,
           [1709.43058496, 1709.43058496, -435.56428767],
           [1709.43058496, 1709.43058496, -605.21461523],
           [1709.43058496, 1709.43058496, -840.942981  ]]), right_voxel_edges=array([[1709.43058496, 1709.43058496,  -30.        ],
           [1709.43058496, 1709.43058496,  -11.68486483],
           [1709.43058496, 1709.43058496,  -16.23606704],
           ...,
           [1709.43058496, 1709.43058496, -605.21461523],
           [1709.43058496, 1709.43058496, -840.942981  ],
           [1709.43058496, 1709.43058496, -840.942981  ]]))



.. GENERATED FROM PYTHON SOURCE LINES 99-101

.. code-block:: Python

    geo_model.grid.centered_grid.kernel_grid_centers





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

 .. code-block:: none


    array([[-5000.        , -5000.        ,  -300.        ],
           [-5000.        , -5000.        ,  -360.        ],
           [-5000.        , -5000.        ,  -383.36972966],
           ...,
           [ 5000.        ,  5000.        , -3407.68480754],
           [ 5000.        ,  5000.        , -4618.11403801],
           [ 5000.        ,  5000.        , -6300.        ]])



.. GENERATED FROM PYTHON SOURCE LINES 102-105

Now we need to compute the component tz (see
https://github.com/cgre-achen/gempy/blob/master/notebooks/tutorials/ch2-2-Cell_selection.ipynb)


.. GENERATED FROM PYTHON SOURCE LINES 107-111

.. code-block:: Python


    gravity_gradient = gp.calculate_gravity_gradient(geo_model.grid.centered_grid)
    gravity_gradient





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

 .. code-block:: none


    array([-0.00435884, -0.0035374 , -0.00260207, ..., -0.60455378,
           -0.888396  , -0.98280245])



.. GENERATED FROM PYTHON SOURCE LINES 112-118

.. code-block:: Python

    geo_model.geophysics_input = gp.data.GeophysicsInput(
        tz=gravity_gradient,
        densities=np.array([2.61, 2.92, 3.1, 2.92, 2.61]),
    )









.. GENERATED FROM PYTHON SOURCE LINES 119-122

Once we have created a gravity interpolator we can call it from compute
model as follows:


.. GENERATED FROM PYTHON SOURCE LINES 124-136

.. code-block:: Python

    geo_model.interpolation_options.mesh_extraction = False
    sol = gp.compute_model(
        gempy_model=geo_model,
        engine_config=gp.data.GemPyEngineConfig(
            backend=gp.data.AvailableBackends.numpy,
            dtype='float32'
        )
    )

    grav = sol.gravity






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

 .. code-block:: none

    Setting Backend To: AvailableBackends.numpy
    Chunking done: 16 chunks
    Chunking done: 16 chunks
    Chunking done: 16 chunks




.. GENERATED FROM PYTHON SOURCE LINES 137-139

.. code-block:: Python

    gpv.plot_2d(geo_model, cell_number=[-1], direction=['z'], show_data=False)




.. image-sg:: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_003.png
   :alt: Cell Number: -1 Direction: z
   :srcset: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_003.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 0x7f9191ce7250>



.. GENERATED FROM PYTHON SOURCE LINES 140-142

.. code-block:: Python

    gpv.plot_2d(geo_model, cell_number=['mid'], direction='x')




.. image-sg:: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_004.png
   :alt: Cell Number: mid Direction: x
   :srcset: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_004.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 0x7f918fbc49d0>



.. GENERATED FROM PYTHON SOURCE LINES 143-153

.. code-block:: Python

    gpv.plot_2d(geo_model, direction=['z'], height=7, show_results=False, show_data=True, show=False)
    plt.scatter(xy_ravel[:, 0], xy_ravel[:, 1], s=1)
    plt.imshow(sol.gravity.reshape(grav_res, grav_res),
               extent=(xy_ravel[:, 0].min() + (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                       xy_ravel[:, 0].max() - (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                       xy_ravel[:, 1].min() + (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2,
                       xy_ravel[:, 1].max() - (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2),
               cmap='viridis_r', origin='lower')
    plt.show()




.. image-sg:: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_005.png
   :alt: Cell Number: mid Direction: z
   :srcset: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 154-160

Plotting lithologies
^^^^^^^^^^^^^^^^^^^^

If we want to compute the lithologies we will need to create a normal
interpolator object as seen in the Chapter 1 of the tutorials


.. GENERATED FROM PYTHON SOURCE LINES 163-166

Now we can plot all together (change the alpha parameter to see the
gravity overlying):


.. GENERATED FROM PYTHON SOURCE LINES 168-183

.. code-block:: Python

    gpv.plot_2d(geo_model, cell_number=[-1], direction=['z'], show=False,
                kwargs_regular_grid={'alpha': .5})

    plt.scatter(xy_ravel[:, 0], xy_ravel[:, 1], s=1)
    plt.imshow(grav.reshape(grav_res, grav_res),
               extent=(xy_ravel[:, 0].min() + (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                       xy_ravel[:, 0].max() - (xy_ravel[0, 0] - xy_ravel[1, 0]) / 2,
                       xy_ravel[:, 1].min() + (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2,
                       xy_ravel[:, 1].max() - (xy_ravel[0, 1] - xy_ravel[30, 1]) / 2),
               cmap='viridis_r', origin='lower', alpha=.8)
    cbar = plt.colorbar()
    cbar.set_label(r'$\mu$gal')
    plt.show()

    # sphinx_gallery_thumbnail_number = -2



.. image-sg:: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_006.png
   :alt: Cell Number: -1 Direction: z
   :srcset: /tutorials/ch2-Geophysics/images/sphx_glr_ch2_1_gravity_006.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_tutorials_ch2-Geophysics_ch2_1_gravity.py:

.. only:: html

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

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

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

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

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


.. only:: html

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

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