Video Tutorial “code-along”: Onlap relations

This tutorial demonstrates step-by-step how to incorporate onlap relations to our geological models created with gempy. It follows the Video tutorial series available on the [gempy YouTube channel](https://www.youtube.com/@GemPy3D). Please follow the first and second part of the tutorials to learn the basics of modeling with gempy before diving into this tutorial.

Video tutorial 11: Basic onlap scenario

# Required imports
import gempy as gp
import gempy_viewer as gpv
import numpy as np

# Path to input data
data_path = "https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/"
path_to_data = data_path + "/data/input_data/video_tutorials_v3/"

# Create instance of geomodel
geo_model_onlap = gp.create_geomodel(
    project_name = 'tutorial_model_onlap_1',
    extent=[0,2000,0,1000,0,1000],
    resolution=[100,50,50],
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=path_to_data+"tutorial_model_onlap_1_orientations.csv",
        path_to_surface_points=path_to_data+"tutorial_model_onlap_1_surface_points.csv"
    )
)

Surface points hash:  2de2814bb46c6cbaa69572685aff54886e4ce3e3620384a0c70e3c52b71c0bd1
Orientations hash:  059fe82183b03d9c50b4eb8249e2ea99297a704ce85138678a7d93e7d19b7b12

# Map geological series to surfaces
gp.map_stack_to_surfaces(
    gempy_model=geo_model_onlap,
    mapping_object={
        "Young_Series": ("basin_fill_2", "basin_fill_1"),
        "Old_Series": ("basin_top", "basin_bottom")
    }
)

# Alternative way of mapping geological series to surfaces
# gp.map_stack_to_surfaces(
#     gempy_model=geo_model_onlap,
#     mapping_object={
#         "Young_Series": ("basin_fill_2", "basin_fill_1"),
#         "Onlap_Series": ("basin_top"),
#         "Old_Series": ("basin_bottom")
#     }
# )

Structural Groups:	StructuralGroup: Name: Young_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: basin_fill_2 StructuralElement: Name: basin_fill_1 StructuralGroup: Name: Old_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: basin_top StructuralElement: Name: basin_bottom
Fault Relations:	Young_Seri... Old_Series Young_Series Old_Series
	True False

# Display a basic cross section of input data
gpv.plot_2d(geo_model_onlap, show_data=True)

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

# Compute a solution for the model
gp.compute_model(geo_model_onlap)

Setting Backend To: AvailableBackends.numpy
Chunking done: 10 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
Chunking done: 19 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))

Solutions: 4 Octree Levels, 4 DualContouringMeshes

# Display the result in 2d section
gpv.plot_2d(geo_model_onlap, show_boundaries=False)

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

# Set the relation of the youngest group to Onlap
from gempy_engine.core.data.stack_relation_type import StackRelationType
geo_model_onlap.structural_frame.structural_groups[0].structural_relation = StackRelationType.ONLAP

# Display updated strucutral frame
geo_model_onlap.structural_frame

Structural Groups:	StructuralGroup: Name: Young_Series Structural Relation: StackRelationType.ONLAP Elements: StructuralElement: Name: basin_fill_2 StructuralElement: Name: basin_fill_1 StructuralGroup: Name: Old_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: basin_top StructuralElement: Name: basin_bottom
Fault Relations:	Young_Seri... Old_Series Young_Series Old_Series
	True False

# Compute a solution for the model
gp.compute_model(geo_model_onlap)

Setting Backend To: AvailableBackends.numpy
Chunking done: 10 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
Chunking done: 19 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))

Solutions: 4 Octree Levels, 4 DualContouringMeshes

# Display the result in 2d section
gpv.plot_2d(geo_model_onlap, show_boundaries=False)

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

Video tutorial 12: Advanced onlap - Subduction zone

# Create instance of geomodel
geo_model_subduction = gp.create_geomodel(
    project_name = 'tutorial_model_onlap_2',
    extent=[0,2000,0,1000,0,1000],
    resolution=[100,50,50],
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=path_to_data+"tutorial_model_onlap_2_orientations.csv?",
        path_to_surface_points=path_to_data+"tutorial_model_onlap_2_surface_points.csv?"
    )
)

Surface points hash:  d0999adb1b7cfca54abd040f517aeba042d7f43a69b507c69ac6ad8cd6682e13
Orientations hash:  ebada610b5a68730bf1d03323594e755fbe986cd07857ee99d8b758716e3dfaa

# Display a basic cross section of input data
gpv.plot_2d(geo_model_subduction)

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

# Map geological series to surfaces
gp.map_stack_to_surfaces(
    gempy_model=geo_model_subduction,
    mapping_object={
        "Top": ("continental_top"),
        "Continental_Series": ("continental_shallow", "continental_deep"),
        "Oceanic_Series": ("oceanic_top", "oceanic_bottom")
    }
)

Structural Groups:	StructuralGroup: Name: Top Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: continental_top StructuralGroup: Name: Continental_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: continental_shallow StructuralElement: Name: continental_deep StructuralGroup: Name: Oceanic_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: oceanic_top StructuralElement: Name: oceanic_bottom
Fault Relations:	Top Continenta... Oceanic_Se... Top Continental_Series Oceanic_Series
	True False

# Set the realtion of the youngest and second youngest group to Onlap
geo_model_subduction.structural_frame.structural_groups[0].structural_relation = StackRelationType.ONLAP
geo_model_subduction.structural_frame.structural_groups[1].structural_relation = StackRelationType.ONLAP

# Display updated structural frame
geo_model_subduction.structural_frame

Structural Groups:	StructuralGroup: Name: Top Structural Relation: StackRelationType.ONLAP Elements: StructuralElement: Name: continental_top StructuralGroup: Name: Continental_Series Structural Relation: StackRelationType.ONLAP Elements: StructuralElement: Name: continental_shallow StructuralElement: Name: continental_deep StructuralGroup: Name: Oceanic_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: oceanic_top StructuralElement: Name: oceanic_bottom
Fault Relations:	Top Continenta... Oceanic_Se... Top Continental_Series Oceanic_Series
	True False

# Create a simple topography using numpy

# Define grid spacing
spacing = 20

# Generate grid
x = np.arange(geo_model_subduction.grid.regular_grid.extent[0], geo_model_subduction.grid.regular_grid.extent[1] + spacing, spacing)
y = np.arange(geo_model_subduction.grid.regular_grid.extent[2], geo_model_subduction.grid.regular_grid.extent[3] + spacing, spacing)
X, Y = np.meshgrid(x, y)

# Define elevation (z) based on x, creating a simple mountain range
Z = np.ones_like(X) * 590  # Default elevation
Z[(X >= 570) & (X < 1000)] = 590 + (200 * (X[(X >= 570) & (X < 1000)] - 600) / 400)
Z[(X >= 1000) & (X < 1300)] = 810 - (250 * (X[(X >= 1000) & (X < 1300)] - 1000) / 300)
Z[X >= 1300] = 540

# Flatten the data into (N,3) shape
topography_points = np.vstack((X.ravel(), Y.ravel(), Z.ravel())).T

# Set topography from numpy array
gp.set_topography_from_arrays(grid=geo_model_subduction.grid, xyz_vertices=topography_points)

Active grids: GridTypes.NONE|TOPOGRAPHY|DENSE

<gempy.core.data.grid_modules.topography.Topography object at 0x7fbcb6c87550>

# Compute a solution for the model
gp.compute_model(geo_model_subduction)

Setting Backend To: AvailableBackends.numpy
Chunking done: 14 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
Chunking done: 10 chunks
Chunking done: 12 chunks
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))
/home/leguark/TeamCity/work/3a8738c25f60c3c9/venv/lib/python3.10/site-packages/gempy_engine/modules/activator/_soft_segment.py:95: RuntimeWarning: overflow encountered in exp
  return 1.0 / (1.0 + bt.t.exp(x))

Solutions: 4 Octree Levels, 5 DualContouringMeshes

# Display the result in 2d section
gpv.plot_2d(geo_model_subduction, show_topography=True, show_boundaries=False)

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

# Display 3d plot of final model
gpv.plot_3d(geo_model_subduction, show_topography=True, image=True)

# sphinx_gallery_thumbnail_number = -1

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7fbcbca71ab0>