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

.. only:: html

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

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

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

.. _sphx_glr_examples_geometries_a01_horizontal_stratigraphic.py:


Model 1 - Horizontal Stratigraphy
===================================

This example demonstrates how to build a basic geological model with horizontally stacked layers using GemPy,
an open-source Python library for implicit geological modeling. In this script, we will:

- Set up and compute a simple geological model using input data (orientations and surface points).
- Visualize the model in 2D and 3D.
- Export model components (formations, orientations, surface points, and the regular grid) to VTK files.
- Extract mesh vertices and normals, save them to Excel files, and convert the vertex data to an XYZ file.

Each section includes detailed comments to explain its purpose and functionality.

.. GENERATED FROM PYTHON SOURCE LINES 17-18

Import necessary libraries

.. GENERATED FROM PYTHON SOURCE LINES 18-26

.. code-block:: Python

    import gempy as gp
    import gempy_viewer as gpv
    import pyvista as pv
    import vtk
    import pandas as pd

    # sphinx_gallery_thumbnail_number = 2








.. GENERATED FROM PYTHON SOURCE LINES 27-32

Generate the Geological Model
-----------------------------
In this section, we define the data source and create a GeoModel instance.
The model is defined with a specific spatial extent and resolution (refinement).
Input data (orientations and surface points) are loaded from a remote repository.

.. GENERATED FROM PYTHON SOURCE LINES 32-59

.. code-block:: Python


    # Define the path to the dataset hosted online
    data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'

    # Create a GeoModel instance with the given project name, extent, and refinement level.
    data = gp.create_geomodel(
        project_name='horizontal',
        extent=[0, 1000, 0, 1000, 0, 1000],
        refinement=6,
        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"
        )
    )

    # Map a geological series to the corresponding surfaces.
    # Here, the "Strat_Series" is associated with two formations ('rock2' and 'rock1'),
    # which defines the stacking order of horizontal layers.
    gp.map_stack_to_surfaces(
        gempy_model=data,
        mapping_object={"Strat_Series": ('rock2', 'rock1')}
    )

    # Compute the geological model using the provided data and mappings.
    gp.compute_model(data)
    geo_data = data





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

 .. code-block:: none

    Surface points hash:  6f1a39ed77e87a4057f03629c946b1876b87e24409cadfe0e1cf7ab1488f69e4
    Orientations hash:  04c307ae23f70252fe54144a2fb95ca7d96584a2d497ea539ed32dfd23e7cd5d
    Setting Backend To: AvailableBackends.numpy




.. GENERATED FROM PYTHON SOURCE LINES 60-66

2D Visualization of the Geological Model
-----------------------------------------
The following plots provide 2D views of the model.
- The first plot shows the initial model without computed results.
- The subsequent plots display the computed geological data in the x and y directions.
Note: The boundaries are hidden in the latter plots for a cleaner visualization.

.. GENERATED FROM PYTHON SOURCE LINES 66-74

.. code-block:: Python


    # Plot the initial model in the y direction (without displaying computed results)
    gpv.plot_2d(geo_data, direction=['y'], show_results=False)

    # Plot the computed model results in the x and y directions, including data but excluding boundaries
    gpv.plot_2d(geo_data, direction=['x'], show_data=True, show_boundaries=False)
    gpv.plot_2d(geo_data, direction=['y'], show_data=True, show_boundaries=False)




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_001.png
         :alt: Cell Number: mid Direction: y
         :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_002.png
         :alt: Cell Number: mid Direction: x
         :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_003.png
         :alt: Cell Number: mid Direction: y
         :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_003.png
         :class: sphx-glr-multi-img


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

 .. code-block:: none


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



.. GENERATED FROM PYTHON SOURCE LINES 75-79

3D Visualization and Export to VTK
----------------------------------
Visualize the model in 3D with data, results, and boundaries.
The generated 3D plot is also used to export various components of the model to VTK files.

.. GENERATED FROM PYTHON SOURCE LINES 79-91

.. code-block:: Python

    p = gpv.plot_3d(geo_data, show_data=True, show_results=True, show_boundaries=True, image=True)

    # Export different components of the model to VTK files:
    p.surface_poly['rock1'].save('rock1.vtk')  # Save formation 'rock1'
    p.surface_poly['rock2'].save('rock2.vtk')  # Save formation 'rock2'
    p.orientations_mesh.save('orientations.vtk')  # Save the orientations mesh
    p.surface_points_mesh.save('surface_points.vtk')  # Save the surface points mesh

    # Retrieve and export the regular grid (volume representation) of the model
    box = p.regular_grid_actor.GetMapper().GetInput()
    box.save('box.vtk')




.. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_004.png
   :alt: a01 horizontal stratigraphic
   :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_004.png
   :class: sphx-glr-single-img

.. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_005.png
   :alt: a01 horizontal stratigraphic
   :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 92-95

Display the Exported VTK Files with PyVista
-------------------------------------------
Load and plot each exported VTK file to verify the export.

.. GENERATED FROM PYTHON SOURCE LINES 95-102

.. code-block:: Python

    pv.read('rock1.vtk').plot(show_edges=False)
    pv.read('rock2.vtk').plot(show_edges=False)
    pv.read('orientations.vtk').plot(show_edges=False)
    pv.read('surface_points.vtk').plot(show_edges=False)
    pv.read('box.vtk').plot(show_edges=False)





.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_006.png
          :alt: a01 horizontal stratigraphic
          :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_006.png
          :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_007.png
          :alt: a01 horizontal stratigraphic
          :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_007.png
          :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_008.png
          :alt: a01 horizontal stratigraphic
          :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_008.png
          :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_009.png
          :alt: a01 horizontal stratigraphic
          :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_009.png
          :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_010.png
          :alt: a01 horizontal stratigraphic
          :srcset: /examples/geometries/images/sphx_glr_a01_horizontal_stratigraphic_010.png
          :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 103-107

Extract and Save Mesh Vertices and Normals
------------------------------------------
The following functions extract vertex coordinates and corresponding normal vectors
from a mesh (vtkPolyData), and then save this data into Excel files for further analysis.

.. GENERATED FROM PYTHON SOURCE LINES 107-193

.. code-block:: Python


    def generate_normals(polydata):
        """
        Generate normals for the given polydata if they are not already computed.

        Parameters:
            polydata (vtk.vtkPolyData): Input polydata for which normals are to be computed.

        Returns:
            vtk.vtkPolyData: The polydata with computed point normals.
        """
        normal_generator = vtk.vtkPolyDataNormals()
        normal_generator.SetInputData(polydata)
        normal_generator.ComputePointNormalsOn()
        normal_generator.ComputeCellNormalsOff()
        normal_generator.Update()
        return normal_generator.GetOutput()


    def get_vertices_and_normals(mesh):
        """
        Extract vertices and normals from a mesh.

        Parameters:
            mesh: The mesh object (typically from a formation) to extract surface vertices and normals.

        Returns:
            tuple: Two lists containing vertices and normals respectively.
        """
        # Extract the surface of the mesh
        surface_mesh = mesh.extract_surface()
        polydata = surface_mesh

        # Ensure normals are computed for the polydata
        polydata_with_normals = generate_normals(polydata)

        # Extract points (vertices)
        points = polydata_with_normals.GetPoints()
        vertices = [points.GetPoint(i) for i in range(points.GetNumberOfPoints())]

        # Extract normals associated with the points
        normals_array = polydata_with_normals.GetPointData().GetNormals()
        normals = [normals_array.GetTuple(i) for i in range(normals_array.GetNumberOfTuples())]

        return vertices, normals


    def save_to_excel(vertices, normals, vertices_file, normals_file):
        """
        Save the vertices and normals to separate Excel files.

        Parameters:
            vertices (list): List of vertex coordinates.
            normals (list): List of normal vectors.
            vertices_file (str): File name for saving vertices.
            normals_file (str): File name for saving normals.
        """
        # Create DataFrames from the lists
        vertices_df = pd.DataFrame(vertices, columns=['X', 'Y', 'Z'])
        normals_df = pd.DataFrame(normals, columns=['x', 'y', 'z'])

        # Write the DataFrames to Excel files
        vertices_df.to_excel(vertices_file, index=False)
        normals_df.to_excel(normals_file, index=False)


    # Extract vertices and normals from the mesh of formation 'rock1'
    mesh = p.surface_poly['rock1']
    vertices, normals = get_vertices_and_normals(mesh)

    # Create DataFrames for later processing or verification
    vertices_df = pd.DataFrame(vertices, columns=['X', 'Y', 'Z'])
    normals_df = pd.DataFrame(normals, columns=['x', 'y', 'z'])

    # Define file names for the Excel outputs
    vertices_file = "rock1_vertices.xlsx"
    normals_file = "rock1_norms.xlsx"

    # Save the extracted data to Excel files
    save_to_excel(vertices, normals, vertices_file, normals_file)

    # Optionally, read back the Excel files to verify their content
    pd.read_excel(vertices_file)
    pd.read_excel(normals_file)







.. 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>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>...</th>
          <td>...</td>
          <td>...</td>
          <td>...</td>
        </tr>
        <tr>
          <th>251</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>252</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>253</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>254</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
        <tr>
          <th>255</th>
          <td>0</td>
          <td>0</td>
          <td>1</td>
        </tr>
      </tbody>
    </table>
    <p>256 rows × 3 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 194-198

Convert Vertices DataFrame to an XYZ File
-----------------------------------------
This function converts a DataFrame containing vertex coordinates to an XYZ format file,
which is a simple text file with one point per line.

.. GENERATED FROM PYTHON SOURCE LINES 198-217

.. code-block:: Python


    def dataframe_to_xyz(df, filename):
        """
        Write vertex coordinates from a DataFrame to an XYZ file.

        Parameters:
            df (pandas.DataFrame): DataFrame containing 'X', 'Y', and 'Z' columns.
            filename (str): Output filename for the XYZ file.
        """
        with open(filename, 'w') as f:
            for index, row in df.iterrows():
                f.write(f"{row['X']} {row['Y']} {row['Z']}\n")


    # Specify the output file name for the XYZ file
    filename = "output.xyz"

    # Convert the vertices DataFrame to an XYZ file
    dataframe_to_xyz(vertices_df, filename)








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

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


.. _sphx_glr_download_examples_geometries_a01_horizontal_stratigraphic.py:

.. only:: html

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

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

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

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

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

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

      :download:`Download zipped: a01_horizontal_stratigraphic.zip <a01_horizontal_stratigraphic.zip>`


.. only:: html

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

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