.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_basic/01_simple_tharsis_model.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_basic_01_simple_tharsis_model.py: Simple Geological Model - Tharsis Region ========================================= This example creates a simple 3D geological model of the Tharsis region using GemPy. **Overview** The Tharsis mining district in southwestern Spain is part of the Iberian Pyrite Belt (IPB), one of the world's most important volcanogenic massive sulfide (VMS) provinces. This example demonstrates the fundamental workflow for creating a 3D geological model from structural data. **Geological Context** The Tharsis deposit consists primarily of Tournaisian-age plutonites that intrude into older Devonian basement rocks. The intrusion created significant hydrothermal alteration and mineralization, making it a target for geophysical exploration methods. **Technical Details** This model uses: * **Octree refinement (level 5)**: Adaptive mesh refinement that concentrates computational effort near geological interfaces while maintaining efficiency in homogeneous regions * **Real structural data**: Surface contact points and orientation measurements from field mapping * **Topographic integration**: High-resolution DEM to accurately represent surface geology .. note:: The coordinate system uses UTM Zone 29N projection (EPSG:32629) with elevations in meters above sea level. .. GENERATED FROM PYTHON SOURCE LINES 34-59 Import Libraries ---------------- We use **GemPy** for 3D implicit geological modeling. GemPy uses a universal co-kriging interpolation approach to construct continuous 3D geological surfaces from sparse data points. **Universal Co-Kriging** GemPy's core interpolation engine is based on Universal Co-Kriging, a geostatistical method that combines multiple types of information into a single consistent 3D model. .. math:: Z(x) = \sum_{i=1}^n \lambda_i Z(x_i) + \sum_{j=1}^m \beta_j f_j(x) Where: * :math:`Z(x)` is the value of the potential field at location :math:`x` * :math:`\sum \lambda_i Z(x_i)` represents the **kriging** part (honoring point data) * :math:`\sum \beta_j f_j(x)` represents the **universal** part (honoring orientations/gradients) This approach ensures that the resulting surfaces are both continuous and honor all structural constraints simultaneously. **NumPy** provides numerical computing capabilities for array operations. .. GENERATED FROM PYTHON SOURCE LINES 59-66 .. code-block:: Python import numpy as np import gempy as gp # Set random seed for reproducibility np.random.seed(1234) .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH .. GENERATED FROM PYTHON SOURCE LINES 67-68 Import paths configuration .. GENERATED FROM PYTHON SOURCE LINES 68-70 .. code-block:: Python from mineye.config import paths .. GENERATED FROM PYTHON SOURCE LINES 71-83 Define Model Extent and Resolution ----------------------------------- The model covers the Tharsis mining district in UTM coordinates (Zone 29N). **Model Extent**: The bounding box defines the 3D volume for interpolation: * **X range**: 33.96 km (East-West direction) * **Y range**: 25.12 km (North-South direction) * **Z range**: 1.005 km (from -500m to +505m elevation) The vertical extent captures both subsurface structures (down to 500m below sea level) and topography (up to 505m above sea level). .. GENERATED FROM PYTHON SOURCE LINES 83-103 .. code-block:: Python min_x = -707_521 # * Cropping the corrupted area of the geotiff max_x = -675558 min_y = 4526832 max_y = 4551949 max_z = 505 model_depth = -500 extent = [min_x, max_x, min_y, max_y, model_depth, max_z] # **Model resolution**: use octree with refinement level 5 # # Octree refinement adaptively subdivides the model space: # # * Level 5 creates up to 2^5 = 32 subdivisions per dimension # * Finer resolution near geological contacts # * Coarser resolution in homogeneous regions # * Significantly more efficient than regular grids for complex geometries resolution = None # Using octrees instead of regular grid refinement = 5 .. GENERATED FROM PYTHON SOURCE LINES 104-118 Get Data Paths -------------- Load paths to structural data **Structural geological data** consists of two types: 1. **Surface points**: 3D coordinates (X, Y, Z) marking where geological contacts were observed, typically from field mapping, drill holes, or cross-sections 2. **Orientations**: The 3D orientation (dip, azimuth) of geological surfaces at specific locations, providing information about how surfaces trend through space These data constrain the implicit interpolation algorithm to construct geologically plausible 3D surfaces. .. GENERATED FROM PYTHON SOURCE LINES 118-125 .. code-block:: Python mod_or_path = paths.get_orientations_path() mod_pts_path = paths.get_points_path() print(f"Orientations: {mod_or_path}") print(f"Points: {mod_pts_path}") .. rst-class:: sphx-glr-script-out .. code-block:: none Orientations: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/Model_Input_Data/Simple-Models/orientations_mod.csv Points: /home/leguark/PycharmProjects/Mineye-Terranigma/examples/Data/Model_Input_Data/Simple-Models/points_mod.csv .. GENERATED FROM PYTHON SOURCE LINES 126-140 Create GemPy Geological Model ------------------------------ Build the model structure with imported structural data **GemPy's implicit modeling approach**: * Uses a universal co-kriging interpolator to construct smooth 3D scalar fields * Each geological surface is represented as an isosurface of a scalar field * The gradient of the scalar field (orientation data) and the isovalue (surface points) constrain the interpolation * The result is a continuous 3D model that honors all input data The ``ImporterHelper`` automatically reads and formats CSV files containing the structural data. .. GENERATED FROM PYTHON SOURCE LINES 140-152 .. code-block:: Python simple_geo_model = gp.create_geomodel( project_name='simple_model', extent=extent, refinement=refinement, resolution=resolution, importer_helper=gp.data.ImporterHelper( path_to_orientations=mod_or_path, path_to_surface_points=mod_pts_path, ) ) .. GENERATED FROM PYTHON SOURCE LINES 153-167 Map Geological Units -------------------- Define the stratigraphic stack **Stratigraphic series** organize geological formations into age-related groups. In GemPy, each series can contain multiple surfaces (formations) that are interpolated together using the same scalar field. For this simple model, we have a single series containing the Tournaisian Plutonites. This intrusive body represents the youngest geological unit in the model, with all other rocks implicitly treated as the "basement" by GemPy. In more complex models, you would define multiple series representing different geological events (e.g., sedimentation, folding, faulting, intrusion). .. GENERATED FROM PYTHON SOURCE LINES 167-175 .. code-block:: Python gp.map_stack_to_surfaces( gempy_model=simple_geo_model, mapping_object={ "Tournaisian_Plutonites": ["Tournaisian Plutonites"], } ) .. raw:: html
Structural Groups: StructuralGroup:
Name:Tournaisian_Plutonites
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Tournaisian Plutonites
Fault Relations:
Tournaisia...
Tournaisian_Plutonites
True
False


.. GENERATED FROM PYTHON SOURCE LINES 176-193 Compute Model with Topography ------------------------------ Add topographic surface from DEM and compute the model **Digital Elevation Models (DEMs)** provide high-resolution surface topography that significantly improves the accuracy of near-surface geological modeling. **Importance of topography in geological modeling**: * Defines the true upper boundary of the model (not just a flat plane) * Critical for geophysical modeling (gravity, magnetics) where topography affects measurements * Enables accurate calculation of true thickness of geological units * Important for resource estimation and mine planning The DEM is loaded from a raster file and resampled to match the model grid. GemPy automatically clips the geological model at the topographic surface, ensuring that geological units only exist below ground. .. GENERATED FROM PYTHON SOURCE LINES 193-209 .. code-block:: Python topography_path = paths.get_topography_path() gp.set_topography_from_file( grid=simple_geo_model.grid, filepath=topography_path, crop_to_extent=[ simple_geo_model.grid.extent[0], simple_geo_model.grid.extent[2], simple_geo_model.grid.extent[1], simple_geo_model.grid.extent[3] ] ) # Compute the model with topography gp.compute_model(simple_geo_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Active grids: GridTypes.OCTREE|TOPOGRAPHY|NONE Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces /home/leguark/.venv/2025/lib/python3.13/site-packages/torch/utils/_device.py:103: UserWarning: Using torch.cross without specifying the dim arg is deprecated. Please either pass the dim explicitly or simply use torch.linalg.cross. The default value of dim will change to agree with that of linalg.cross in a future release. (Triggered internally at /pytorch/aten/src/ATen/native/Cross.cpp:63.) return func(*args, **kwargs) .. raw:: html
Solutions: 5 Octree Levels, 1 DualContouringMeshes


.. GENERATED FROM PYTHON SOURCE LINES 210-220 Visualize the Model ------------------- **2D Cross-Sections** Before diving into 3D, it's often helpful to inspect 2D cross-sections. This allows for a clear view of the internal structural relationships and the behavior of the interpolator along specific planes. We'll plot a cross-section along the Y-axis (North-South) at the center of the model. .. GENERATED FROM PYTHON SOURCE LINES 220-225 .. code-block:: Python import gempy_viewer as gpv gpv.plot_2d(simple_geo_model, direction='y', cell_number='mid', show_data=True) .. image-sg:: /examples_basic/images/sphx_glr_01_simple_tharsis_model_001.png :alt: Cell Number: mid Direction: y :srcset: /examples_basic/images/sphx_glr_01_simple_tharsis_model_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 226-242 **3D Visualization** Create 3D visualization of the geological model **GemPy Viewer** provides interactive 3D visualization capabilities using PyVista. Visualization features: * **3D lithology blocks**: Colored volumes representing different geological units * **Topographic surface**: The DEM draped over the model * **Input data**: Surface points and orientation measurements shown in 3D space * **Vertical exaggeration (ve=5)**: Stretches the Z-axis by a factor of 5 to enhance the visibility of subtle geological features The visualization is fully interactive, allowing rotation, zooming, and inspection of the geological model from any angle. .. GENERATED FROM PYTHON SOURCE LINES 242-253 .. code-block:: Python import gempy_viewer as gpv gpv.plot_3d(simple_geo_model, ve=5, image=False, kwargs_pyvista_bounds={ 'show_xlabels': False, 'show_ylabels': False, 'show_zlabels': False } ) .. image-sg:: /examples_basic/images/sphx_glr_01_simple_tharsis_model_002.png :alt: 01 simple tharsis model :srcset: /examples_basic/images/sphx_glr_01_simple_tharsis_model_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/leguark/PycharmProjects/gempy_viewer/gempy_viewer/modules/plot_3d/drawer_surfaces_3d.py:38: PyVistaDeprecationWarning: ../../../gempy_viewer/gempy_viewer/modules/plot_3d/drawer_surfaces_3d.py:38: Argument 'color' must be passed as a keyword argument to function 'BasePlotter.add_mesh'. From version 0.50, passing this as a positional argument will result in a TypeError. gempy_vista.surface_actors[element.name] = gempy_vista.p.add_mesh( .. GENERATED FROM PYTHON SOURCE LINES 254-285 Summary and Next Steps ---------------------- In this example, we demonstrated the complete workflow for creating a simple 3D geological model using GemPy: 1. ✓ Defined the model extent and resolution 2. ✓ Imported structural geological data (surface points and orientations) 3. ✓ Created a GemPy geomodel with stratigraphic organization 4. ✓ Computed the 3D geological model using implicit interpolation 5. ✓ Integrated high-resolution topography from a DEM 6. ✓ Visualized the resulting 3D model **Key Takeaways**: * GemPy uses implicit modeling to create continuous 3D geological surfaces from sparse data * Octree refinement provides efficient adaptive resolution * Topography integration is crucial for accurate geological modeling * The resulting model can be used for further analysis (geophysics, resource estimation, etc.) **Next Steps**: * Example 02: Forward gravity modeling using this geological model * Example 03: Bayesian segmentation for lithological mapping * Probabilistic examples: Uncertainty quantification and joint inversion .. seealso:: * `GemPy Documentation `_ * `PyVista for 3D Visualization `_ sphinx_gallery_thumbnail_number = 2 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 36.234 seconds) .. _sphx_glr_download_examples_basic_01_simple_tharsis_model.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01_simple_tharsis_model.ipynb <01_simple_tharsis_model.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01_simple_tharsis_model.py <01_simple_tharsis_model.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01_simple_tharsis_model.zip <01_simple_tharsis_model.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_