.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_basic/02_complex_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_02_complex_tharsis_model.py: Complex Geological Model - Tharsis Region ========================================= This example creates a complex 3D geological model of the Tharsis region using GemPy, demonstrating how to model erosive contacts between igneous intrusions and sedimentary sequences. **Overview** Building on the simple model from Example 01, this example demonstrates a more sophisticated geological scenario where Tournaisian Plutonites intrude through a conformable Devonian sedimentary sequence. The key challenge is modeling the **erosive contact** between the intrusive body and the host rocks. **Geological Context** The stratigraphy from youngest to oldest: 1. **Tournaisian Plutonites**: Late Carboniferous intrusive body with erosive contacts 2. **Visean Shales**: Marine shales deposited during the Visean stage 3. **Mid Devonian Siliciclastics**: Clastic sediments from the Middle Devonian 4. **Famennian Siliciclastics**: Upper Devonian basement rocks **Technical Approach** This example demonstrates a **two-model workflow**: 1. First, a stratigraphic model is created for the conformable sedimentary sequence 2. Then, a separate model is created for the plutonite intrusion 3. Finally, the models are merged by overwriting the sedimentary lithologies where the plutonite exists **Why two models instead of one?** Standard implicit modeling (using a single scalar field) can sometimes struggle with **erosive contacts** and complex topological relationships. By using two separate models, we gain several advantages: * **Implicit vs. Explicit Topology**: Instead of relying on the interpolator to correctly handle the cutting relationship (implicit), we explicitly define it during the merge step. This ensures the plutonite always "wins" and cuts cleanly through the host rock. * **Independent Constraints**: The sedimentary layers are constrained by their own orientations and surface points, while the plutonite is constrained by its unique geometry. This prevents data from one domain from "bleeding" into and distorting the other. * **Topological Robustness**: It ensures that no "floating" stratigraphic layers appear inside the intrusive body, which can happen in single-model workflows if the interpolated surfaces cross in unphysical ways. This approach allows each geological domain to be modeled with appropriate constraints while still producing a unified 3D model. .. note:: The coordinate system uses UTM Zone 29N projection (EPSG:32629) with elevations in meters above sea level. .. GENERATED FROM PYTHON SOURCE LINES 61-66 Import Libraries ---------------- We use **GemPy** for 3D implicit geological modeling and **gempy_viewer** for visualization. The ``helper_plotter`` module provides custom visualization functions for combined models. .. GENERATED FROM PYTHON SOURCE LINES 66-74 .. code-block:: Python import numpy as np import gempy as gp import gempy_viewer as gpv # 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 75-88 Define Consistent Color Scheme ------------------------------ **Consistent visualization** is crucial when comparing multiple models and plots. We define a color dictionary that will be used across all visualizations to ensure that each geological unit always appears in the same color. Color scheme rationale: * **Red** for Plutonites: Warm color representing igneous/intrusive rocks * **Blue** for Visean Shales: Cool color for marine sediments * **Green** for Mid Devonian: Intermediate color for older sediments * **Orange** for Famennian: Distinct color for basement rocks .. GENERATED FROM PYTHON SOURCE LINES 90-96 Import Paths and Helper Modules ------------------------------- The ``paths`` module provides centralized access to data file locations. The ``helper_plotter`` module contains custom visualization functions. The ``example_parameters`` module contains shared configuration like color schemes. .. GENERATED FROM PYTHON SOURCE LINES 96-103 .. code-block:: Python from mineye.config import paths from mineye.config.example_parameters import TharsisModelConfig from mineye.GeoModel import helper_plotter FORMATION_COLORS = TharsisModelConfig.TharsisDataProcessingConfig.FORMATION_COLORS .. GENERATED FROM PYTHON SOURCE LINES 104-116 Define Model Extent and Resolution ----------------------------------- **Model Extent**: The bounding box is tightly fit to the data coverage: * **X range**: 27.5 km (East-West direction) * **Y range**: 18.5 km (North-South direction) * **Z range**: 1.0 km (from -500m to +500m elevation) **Resolution**: We use a regular 64×64×64 grid instead of octrees for this model. This is necessary because we need to merge two separate models, which requires compatible regular grids. .. GENERATED FROM PYTHON SOURCE LINES 116-130 .. code-block:: Python min_x = -707500 max_x = -680000 min_y = 4530500 max_y = 4549000 max_z = 500 model_depth = -500 extent = [min_x, max_x, min_y, max_y, model_depth, max_z] # Use regular grid for model merging compatibility resolution = [64, 64, 64] refinement = 5 .. GENERATED FROM PYTHON SOURCE LINES 131-140 Load Structural Data Paths -------------------------- For the complex model, we have separate data files for: * **Sedimentary formations**: Contact points and orientations for the Devonian sequence * **Plutonite intrusion**: Contact points and orientations for the Tournaisian intrusion This separation allows independent quality control and adjustment of each dataset. .. GENERATED FROM PYTHON SOURCE LINES 140-147 .. code-block:: Python mod_or_sed_path = paths.get_orientations_path_sed_complex() mod_pts_sed_path = paths.get_points_path_sed_complex() mod_or_plut_path = paths.get_orientations_path_magmatic_complex() mod_pts_plut_path = paths.get_points_path_magmatic_complex() .. GENERATED FROM PYTHON SOURCE LINES 148-162 Create Stratigraphic Model -------------------------- **Step 1: Model the conformable sedimentary sequence** The sedimentary formations are modeled as a single stratigraphic series with conformable contacts. This means GemPy will interpolate smooth, parallel surfaces that honor the depositional geometry. The stratigraphic order (youngest to oldest): 1. Visean Shales 2. Mid Devonian Siliciclastics 3. Famennian Siliciclastics (basement) .. GENERATED FROM PYTHON SOURCE LINES 162-179 .. code-block:: Python stratigraphic_geo_model = gp.create_geomodel( project_name='stratigraphic_stack_model', extent=extent, refinement=refinement, resolution=resolution, importer_helper=gp.data.ImporterHelper( path_to_orientations=mod_or_sed_path, path_to_surface_points=mod_pts_sed_path, ) ) # Apply consistent colors to stratigraphic surfaces for element in stratigraphic_geo_model.structural_frame.structural_elements: if element.name in FORMATION_COLORS: element.color = FORMATION_COLORS[element.name] .. GENERATED FROM PYTHON SOURCE LINES 180-185 Define Stratigraphic Stack -------------------------- **Stratigraphic series** organize formations by their age relationships. All formations in ``Strat_Series1`` are conformable (parallel contacts). .. GENERATED FROM PYTHON SOURCE LINES 185-193 .. code-block:: Python gp.map_stack_to_surfaces( gempy_model=stratigraphic_geo_model, mapping_object={ "Strat_Series1": ("Visean Shales", "Mid Devonian Siliciclastics", "Famennian Siliciclastics") } ) .. raw:: html
Structural Groups: StructuralGroup:
Name:Strat_Series1
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Visean Shales

StructuralElement:
Name:Mid Devonian Siliciclastics

StructuralElement:
Name:Famennian Siliciclastics
Fault Relations:
Strat_Seri...
Strat_Series1
True
False


.. GENERATED FROM PYTHON SOURCE LINES 194-199 Add Topography and Compute Stratigraphic Model ----------------------------------------------- **Topography integration** ensures the model accurately represents surface geology. The DEM is cropped to match the model extent. .. GENERATED FROM PYTHON SOURCE LINES 199-215 .. code-block:: Python topography_path = paths.get_topography_path() gp.set_topography_from_file( grid=stratigraphic_geo_model.grid, filepath=topography_path, crop_to_extent=[ stratigraphic_geo_model.grid.extent[0], stratigraphic_geo_model.grid.extent[2], stratigraphic_geo_model.grid.extent[1], stratigraphic_geo_model.grid.extent[3] ] ) # Compute the stratigraphic model gp.compute_model(stratigraphic_geo_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Active grids: GridTypes.DENSE|TOPOGRAPHY|NONE Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 3 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, 3 DualContouringMeshes


.. GENERATED FROM PYTHON SOURCE LINES 216-225 Create Plutonite Model ---------------------- **Step 2: Model the intrusive body separately** The Tournaisian Plutonite has an **erosive contact** with the host rocks, meaning it cuts through the pre-existing stratigraphy. By modeling it separately, we can ensure the intrusion geometry is controlled by its own structural data without being influenced by the sedimentary contacts. .. GENERATED FROM PYTHON SOURCE LINES 225-255 .. code-block:: Python plutonite_id = 1 plutonite_geo_model = gp.create_geomodel( project_name='plutonite_model', extent=extent, refinement=refinement, resolution=resolution, importer_helper=gp.data.ImporterHelper( path_to_orientations=mod_or_plut_path, path_to_surface_points=mod_pts_plut_path, ) ) # Map the plutonite surface gp.map_stack_to_surfaces( gempy_model=plutonite_geo_model, mapping_object={ "Plutonite_Series": ["Tournaisian Plutonites"] } ) # Apply consistent color for element in plutonite_geo_model.structural_frame.structural_elements: if element.name in FORMATION_COLORS: element.color = FORMATION_COLORS[element.name] # Compute the plutonite model gp.compute_model(plutonite_geo_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH Using sequential processing for 1 surfaces .. raw:: html
Solutions: 5 Octree Levels, 1 DualContouringMeshes


.. GENERATED FROM PYTHON SOURCE LINES 256-269 Merge Models ------------ **Step 3: Combine the stratigraphic and plutonite models** The merging process: 1. Extract the lithology blocks from both models (3D arrays of formation IDs) 2. Create a mask identifying voxels where the plutonite exists 3. Overwrite the stratigraphic lithologies with the plutonite ID where the mask is True This approach implements the **erosive contact** - the plutonite "cuts through" the pre-existing sedimentary rocks. .. GENERATED FROM PYTHON SOURCE LINES 269-286 .. code-block:: Python # Get plutonite lithology block plut_lith_block = plutonite_geo_model.solutions.raw_arrays.lith_block plut_lith_block_reshaped = plut_lith_block.reshape(64, 64, 64) plutonite_mask = plut_lith_block_reshaped == plutonite_id # Get stratigraphic lithology block lith_block = stratigraphic_geo_model.solutions.raw_arrays.lith_block lith_block_reshaped = lith_block.reshape(64, 64, 64) # Get voxel coordinates for plotting voxel_coords = stratigraphic_geo_model.grid.regular_grid.values # Insert plutonite into stratigraphic model (erosive contact) lith_block_reshaped[plutonite_mask] = 6 lith_block_modified = lith_block_reshaped.flatten() .. GENERATED FROM PYTHON SOURCE LINES 287-299 Visualize Combined Model ------------------------ **Final 3D visualization** showing all geological units together. The plot respects topography by masking voxels above the surface. Formation ID mapping for visualization: * ID 1: Visean Shales * ID 2: Mid Devonian Siliciclastics * ID 3: Famennian Siliciclastics * ID 6: Tournaisian Plutonites (merged from separate model) .. GENERATED FROM PYTHON SOURCE LINES 299-329 .. code-block:: Python FORMATION_ID_MAP = { 1: 'Visean Shales', 2: 'Mid Devonian Siliciclastics', 3: 'Famennian Siliciclastics', 6: 'Tournaisian Plutonites', } # Get topography points for masking (X, Y, Z) topography_points = stratigraphic_geo_model.grid.topography.values p = gpv.plot_2d( stratigraphic_geo_model, section_names=['topography'], # this triggers the top-down geological map show_topography=True, show_lith=True, show_boundaries=True, show_data=False ) # Plot the combined model with topography masking helper_plotter.plot_combined_model( lith_block=lith_block_modified, voxel_coords=voxel_coords, formation_id_map=FORMATION_ID_MAP, formation_colors=FORMATION_COLORS, topography_points=topography_points, title='Combined Geological Model - Tharsis Region' ) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples_basic/images/sphx_glr_02_complex_tharsis_model_001.png :alt: Geological map :srcset: /examples_basic/images/sphx_glr_02_complex_tharsis_model_001.png :class: sphx-glr-multi-img * .. image-sg:: /examples_basic/images/sphx_glr_02_complex_tharsis_model_002.png :alt: Combined Geological Model - Tharsis Region :srcset: /examples_basic/images/sphx_glr_02_complex_tharsis_model_002.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/leguark/PycharmProjects/gempy_viewer/gempy_viewer/API/_plot_2d_sections_api.py:112: UserWarning: Section contacts not implemented yet. We need to pass scalar field for the sections grid warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 330-341 .. code-block:: Python gpv.plot_3d( model=stratigraphic_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_02_complex_tharsis_model_003.png :alt: 02 complex tharsis model :srcset: /examples_basic/images/sphx_glr_02_complex_tharsis_model_003.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 342-362 Summary and Next Steps ---------------------- **Key Takeaways**: * Complex geological relationships can be modeled by combining multiple GemPy models * Erosive contacts are implemented by overwriting lithology values * Consistent color schemes are essential for comparing multiple visualizations * Topography masking provides geologically realistic surface views **Next Steps**: * Example 03: Forward gravity modeling using this geological model * Example 04: Bayesian segmentation for lithological mapping * Probabilistic examples: Uncertainty quantification and joint inversion .. seealso:: * `GemPy Documentation `_ * `PyVista for 3D Visualization `_ sphinx_gallery_thumbnail_number = 1 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.502 seconds) .. _sphx_glr_download_examples_basic_02_complex_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: 02_complex_tharsis_model.ipynb <02_complex_tharsis_model.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02_complex_tharsis_model.py <02_complex_tharsis_model.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02_complex_tharsis_model.zip <02_complex_tharsis_model.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_