Note
Go to the end to download the full example code.
Video Tutorial “code-along”: Modeling step by step¶
This tutorial demonstrates step-by-step geological modeling using the gempy and gempy_viewer libraries. It follows the Video tutorial series available on the gempy YouTube channel.
Video tutorial 1: Introduction¶
The first video is an introduction to GemPy - please view online before starting the tutorial.
Video tutorial 2: Input data¶
# Required imports
import gempy as gp
import gempy_viewer as gpv
# 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 = gp.create_geomodel(
project_name = 'tutorial_model',
extent=[0,2500,0,1000,0,1000],
resolution=[100,40,40],
importer_helper=gp.data.ImporterHelper(
path_to_orientations=path_to_data+"tutorial_model_orientations.csv",
path_to_surface_points=path_to_data+"tutorial_model_surface_points.csv"
)
)
Surface points hash: 7c6d3e04ab03a4b8324d9c91d56c30f9e6a7cb6c22c6f2ee69a5dd001c63337a
Orientations hash: 63e42d294dec66b4db2f175bc7b58553ee89d68f3072d36402963c90b0ef5262
# Display a basic cross section of input data
gpv.plot_2d(geo_model)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3eec88cec0>
# Manually add a surface point
gp.add_surface_points(
geo_model=geo_model,
x=[2250],
y=[500],
z=[750],
elements_names=['rock1']
)
# Show added point in cross section
gpv.plot_2d(geo_model)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3eea7065d0>
Video tutorial 3: Structural frame¶
# View structural frame
geo_model.structural_frame
# View structural elements
geo_model.structural_frame.structural_elements
[Element(
name=rock1,
color=#015482,
is_active=True
), Element(
name=rock2,
color=#9f0052,
is_active=True
), Element(
name=rock3,
color=#ffbe00,
is_active=True
), Element(
name=basement,
color=#728f02,
is_active=True
)]
# Define structural groups and age/stratigraphic relationship
gp.map_stack_to_surfaces(
gempy_model=geo_model,
mapping_object={
"Strat_Series2": ("rock3"),
"Strat_Series1": ("rock2", "rock1")
}
)
Video tutorial 4: Computation and results¶
# View interpolation options
geo_model.interpolation_options
InterpolationOptions(kernel_options=KernelOptions(range=1.7, c_o=10.0, uni_degree=1, i_res=4.0, gi_res=2.0, number_dimensions=3, kernel_function=AvailableKernelFunctions.cubic, kernel_solver=Solvers.DEFAULT, compute_condition_number=False, optimizing_condition_number=False, condition_number=None), evaluation_options=EvaluationOptions(_number_octree_levels=1, _number_octree_levels_surface=4, octree_curvature_threshold=-1.0, octree_error_threshold=1.0, octree_min_level=2, mesh_extraction=True, mesh_extraction_masking_options=<MeshExtractionMaskingOptions.INTERSECT: 3>, mesh_extraction_fancy=True, evaluation_chunk_size=500000, compute_scalar=True, compute_scalar_gradient=False, verbose=False), debug=False, cache_mode=<CacheMode.IN_MEMORY_CACHE: 3>, cache_model_name='tutorial_model', block_solutions_type=<BlockSolutionType.DENSE_GRID: 2>, sigmoid_slope=5000000)
# Compute a solution for the model
gp.compute_model(geo_model)
Setting Backend To: AvailableBackends.PYTORCH
GPU enabled. Using device: cuda
GPU device count: 1
Current GPU device: 0
Chunking done: 7 chunks
Chunking done: 30 chunks
# Display the result in 2d section
gpv.plot_2d(geo_model, cell_number=20)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3eea788190>
# Some examples of how to access results
print(geo_model.solutions.raw_arrays.lith_block)
print(geo_model.grid.dense_grid.values)
[4 4 4 ... 1 1 1]
[[ 12.5 12.5 12.5]
[ 12.5 12.5 37.5]
[ 12.5 12.5 62.5]
...
[2487.5 987.5 937.5]
[2487.5 987.5 962.5]
[2487.5 987.5 987.5]]
Video tutorial 5: 2D visualization¶
# 2d plotting options
gpv.plot_2d(geo_model, show_value=True, show_lith=False, show_scalar=True, series_n=1, cell_number=25)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3eec8c22c0>
# Create custom section lines
gp.set_section_grid(
grid=geo_model.grid,
section_dict={
'section1': ([0, 0], [2500, 1000], [100, 50]),
'section2': ([1000, 1000], [1500, 0], [100, 100]),
}
)
Active grids: GridTypes.DENSE|SECTIONS|NONE
# Show custom cross-section traces
gpv.plot_section_traces(geo_model)
<function plot_section_traces at 0x7f3efa82a980>
# Recompute model as a new grid was added
gp.compute_model(geo_model)
Setting Backend To: AvailableBackends.PYTORCH
GPU enabled. Using device: cuda
GPU device count: 1
Current GPU device: 0
Chunking done: 8 chunks
Chunking done: 33 chunks
# Display custom cross-sections
gpv.plot_2d(geo_model, section_names=['section1', 'section2'], show_data=False)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3eec9b68d0>
Video tutorial 6: 3D visualization¶
# Display the result in 3d
gpv.plot_3d(geo_model, show_lith=True, show_boundaries=True, ve=None)
<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f3eea7d2e40>
# How to access DC meshes
geo_model.solutions.dc_meshes[0].dc_data
# transform mesh vertices to original coordinate system
back_transformed_vertices = geo_model.input_transform.apply_inverse(geo_model.solutions.dc_meshes[0].vertices)
Video tutorial 7: Topography¶
# Setting a randomly generated topography
import numpy as np
gp.set_topography_from_random(
grid=geo_model.grid,
fractal_dimension=1.2,
d_z=np.array([700, 900]),
topography_resolution=np.array([250, 100])
)
Active grids: GridTypes.DENSE|TOPOGRAPHY|SECTIONS|NONE
Topography(_regular_grid=RegularGrid(resolution=array([100, 40, 40]), extent=array([ 0., 2500., 0., 1000., 0., 1000.]), values=array([[ 12.5, 12.5, 12.5],
[ 12.5, 12.5, 37.5],
[ 12.5, 12.5, 62.5],
...,
[2487.5, 987.5, 937.5],
[2487.5, 987.5, 962.5],
[2487.5, 987.5, 987.5]], shape=(160000, 3)), mask_topo=array([], shape=(0, 3), dtype=bool), _transform=None, _base_resolution=array([2, 2, 2])), values_2d=array([[[ 0. , 0. , 846.48200887],
[ 0. , 10.1010101 , 845.7885738 ],
[ 0. , 20.2020202 , 845.13500395],
...,
[ 0. , 979.7979798 , 858.17962427],
[ 0. , 989.8989899 , 858.56778738],
[ 0. , 1000. , 859.02124602]],
[[ 10.04016064, 0. , 847.36106954],
[ 10.04016064, 10.1010101 , 846.65906498],
[ 10.04016064, 20.2020202 , 845.99218121],
...,
[ 10.04016064, 979.7979798 , 858.79922685],
[ 10.04016064, 989.8989899 , 859.15630107],
[ 10.04016064, 1000. , 859.59405788]],
[[ 20.08032129, 0. , 848.23831087],
[ 20.08032129, 10.1010101 , 847.54698227],
[ 20.08032129, 20.2020202 , 846.88378095],
...,
[ 20.08032129, 979.7979798 , 859.36667512],
[ 20.08032129, 989.8989899 , 859.73843469],
[ 20.08032129, 1000. , 860.20142948]],
...,
[[2479.91967871, 0. , 843.59900934],
[2479.91967871, 10.1010101 , 842.92517567],
[2479.91967871, 20.2020202 , 842.26117265],
...,
[2479.91967871, 979.7979798 , 856.39660819],
[2479.91967871, 989.8989899 , 856.81354684],
[2479.91967871, 1000. , 857.26684359]],
[[2489.95983936, 0. , 844.67814891],
[2489.95983936, 10.1010101 , 843.99656816],
[2489.95983936, 20.2020202 , 843.31958929],
...,
[2489.95983936, 979.7979798 , 857.0343389 ],
[2489.95983936, 989.8989899 , 857.43709416],
[2489.95983936, 1000. , 857.87090113]],
[[2500. , 0. , 845.62858323],
[2500. , 10.1010101 , 844.92906332],
[2500. , 20.2020202 , 844.25939646],
...,
[2500. , 979.7979798 , 857.60270315],
[2500. , 989.8989899 , 858.02933643],
[2500. , 1000. , 858.47392974]]],
shape=(250, 100, 3)), source=None, values=array([[ 0. , 0. , 846.48200887],
[ 0. , 10.1010101 , 845.7885738 ],
[ 0. , 20.2020202 , 845.13500395],
...,
[2500. , 979.7979798 , 857.60270315],
[2500. , 989.8989899 , 858.02933643],
[2500. , 1000. , 858.47392974]], shape=(25000, 3)), resolution=(250, 100), raster_shape=())
# Recompute model as a new grid was added
gp.compute_model(geo_model)
Setting Backend To: AvailableBackends.PYTORCH
GPU enabled. Using device: cuda
GPU device count: 1
Current GPU device: 0
Chunking done: 9 chunks
Chunking done: 37 chunks
# Display a cross-section with topography
gpv.plot_2d(geo_model, show_topography=True)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3eec8a6be0>
# Displaying a geological map
gpv.plot_2d(geo_model, show_topography=True, section_names=['topography'], show_boundaries=False, show_data=False)
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f3e5468df20>
# Display the 3d model with topography
gpv.plot_3d(geo_model, show_lith=True, show_topography=True)
# sphinx_gallery_thumbnail_number = -1
<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f3e541e5950>
Total running time of the script: (0 minutes 15.259 seconds)
