.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/b_fundamentals/a01_basics.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_tutorials_b_fundamentals_a01_basics.py: Basics: Data Classes and Modeling API ======================================== A closer look at gempy's data classes and modeling API This tutorial is a written companion to the video tutorials, covering the same simple fault model but going into more technical depth: the data classes gempy is built on, constructing structural elements and groups directly rather than through a CSV import, inspecting a computed model's solutions and meshes, and saving a model to disk. .. GENERATED FROM PYTHON SOURCE LINES 14-19 .. code-block:: Python import numpy as np import gempy as gp import gempy_viewer as gpv .. GENERATED FROM PYTHON SOURCE LINES 20-37 gempy's data classes ----------------------- gempy uses a small set of Python classes to store everything that goes into a model: - :obj:`gempy.core.data.GeoModel` - :obj:`gempy.core.data.StructuralFrame` - :obj:`gempy.core.data.StructuralGroup` - :obj:`gempy.core.data.StructuralElement` - :obj:`gempy.core.data.SurfacePointsTable` - :obj:`gempy.core.data.OrientationsTable` - :obj:`gempy.core.data.Grid` A ``GeoModel`` holds one ``StructuralFrame``, which is an ordered list of ``StructuralGroup`` objects (also called series or stacks), each containing one or more ``StructuralElement`` objects -- a lithological unit or a fault surface, defined by a ``SurfacePointsTable`` and an ``OrientationsTable``. The rest of this tutorial builds up a model through these classes and looks at each one along the way. .. GENERATED FROM PYTHON SOURCE LINES 39-50 Model setup ------------- Surface points mark the **bottom** of a layer (if you need the top of a formation -- modeling an intrusion, say -- use an inverted orientation instead). Data can be supplied from CSV files, as here, or built up point by point in code, which the next tutorial covers. The model's ``extent`` defines the volume used for interpolation and plotting, and should enclose all the input data. ``refinement`` sets the number of octree levels used to extract smooth surfaces (see the Grids tutorial for the full explanation of how this interacts with ``resolution``). .. GENERATED FROM PYTHON SOURCE LINES 52-66 .. code-block:: Python data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' geo_model = gp.create_geomodel( project_name='Tutorial_Basics', extent=[0, 2000, 0, 2000, 0, 750], refinement=6, importer_helper=gp.data.ImporterHelper( path_to_orientations=data_path + "/data/input_data/getting_started/simple_fault_model_orientations.csv", path_to_surface_points=data_path + "/data/input_data/getting_started/simple_fault_model_points.csv", hash_surface_points="4cdd54cd510cf345a583610585f2206a2936a05faaae05595b61febfc0191563", hash_orientations="7ba1de060fc8df668d411d0207a326bc94a6cdca9f5fe2ed511fd4db6b3f3526" ) ) .. rst-class:: sphx-glr-script-out .. code-block:: none Surface points hash: 4cdd54cd510cf345a583610585f2206a2936a05faaae05595b61febfc0191563 Orientations hash: 7ba1de060fc8df668d411d0207a326bc94a6cdca9f5fe2ed511fd4db6b3f3526 .. GENERATED FROM PYTHON SOURCE LINES 67-75 ``ImporterHelper`` bundles everything needed to import data from various sources -- here, CSV files fetched over HTTP and verified against a known hash, matching every other tutorial in this documentation. Reviewing the imported data ------------------------------ The raw imported points and orientations are available as ``surface_points_copy`` and ``orientations_copy``: .. GENERATED FROM PYTHON SOURCE LINES 77-79 .. code-block:: Python geo_model.surface_points_copy .. raw:: html
XYZidnugget
700.001000.00300.00497450340.00
600.001000.00200.00497450340.00
500.001000.00100.00497450340.00
1000.001000.00600.00497450340.00
1100.001000.00700.00497450340.00
1000.0050.00350.001829037370.00
1000.00150.00333.331829037370.00
1000.00300.00333.331829037370.00
1000.00500.00366.671829037370.00
1000.001000.00433.331829037370.00
...............
1100.001700.00473.334262594380.00
1100.001950.00490.004262594380.00
0.001000.00433.334262594380.00
300.001000.00400.004262594380.00
600.001000.00366.674262594380.00
1300.001000.00566.674262594380.00
1600.001000.00550.004262594380.00
1900.001000.00566.674262594380.00
1700.00500.00533.334262594380.00
1700.001500.00516.674262594380.00


.. GENERATED FROM PYTHON SOURCE LINES 80-82 .. code-block:: Python geo_model.orientations_copy .. raw:: html
XYZG_xG_yG_zidnugget
500.001000.00300.00-0.95-0.000.32497450340.01
400.001000.00420.000.320.000.952470168680.01
1000.001000.00300.000.320.000.953323137030.01


.. GENERATED FROM PYTHON SOURCE LINES 83-87 Each structural element is internally tracked by a numeric ID. Note these aren't the small, sequential IDs used to color the lithology block in plots -- they're derived directly from each element's name and used for tracking identity regardless of reordering. ``element_id_name_map`` looks up which ID corresponds to which element: .. GENERATED FROM PYTHON SOURCE LINES 89-91 .. code-block:: Python geo_model.structural_frame.element_id_name_map .. rst-class:: sphx-glr-script-out .. code-block:: none {np.int32(49745034): 'Main_Fault', np.int32(182903737): 'Sandstone_1', np.int32(247016868): 'Sandstone_2', np.int32(332313703): 'Shale', np.int32(426259438): 'Siltstone', 39541672: 'basement'} .. GENERATED FROM PYTHON SOURCE LINES 92-106 Structural groups and series ------------------------------- Geological units need to appear in the correct chronological order -- a sequence of deposition, unconformities, intrusions, and so on. In gempy this is expressed by assigning each unit (and each fault) to a **structural group**, using ``map_stack_to_surfaces``. Units in the same group share one continuous scalar field, so the order *within* a group only affects the default color; the order *between* groups is what encodes geological age, oldest at the bottom. Faults are always their own group and must be younger than whatever they affect. Where multiple faults are involved, their relative order encodes their tectonic relationship (the first entry is the youngest). This model has one fault and four stratigraphic layers, assigned to two groups: .. GENERATED FROM PYTHON SOURCE LINES 108-116 .. code-block:: Python gp.map_stack_to_surfaces( gempy_model=geo_model, mapping_object={ "Fault_Series": 'Main_Fault', "Strat_Series": ('Sandstone_2', 'Siltstone', 'Shale', 'Sandstone_1') } ) .. raw:: html
Structural Groups: StructuralGroup:
Name:Fault_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Main_Fault

StructuralGroup:
Name:Strat_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Sandstone_2

StructuralElement:
Name:Siltstone

StructuralElement:
Name:Shale

StructuralElement:
Name:Sandstone_1
Fault Relations:
Fault_Seri...Strat_Seri...
Fault_Series
Strat_Series
True
False


.. GENERATED FROM PYTHON SOURCE LINES 117-120 ``map_stack_to_surfaces`` doesn't yet mark ``Fault_Series`` as a fault -- every group defaults to an ``ERODE`` relation (the next section explains what that means). ``set_is_fault`` does that: .. GENERATED FROM PYTHON SOURCE LINES 122-124 .. code-block:: Python gp.set_is_fault(geo_model, ["Fault_Series"]) .. raw:: html
Structural Groups: StructuralGroup:
Name:Fault_Series
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Main_Fault

StructuralGroup:
Name:Strat_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Sandstone_2

StructuralElement:
Name:Siltstone

StructuralElement:
Name:Shale

StructuralElement:
Name:Sandstone_1
Fault Relations:
Fault_Seri...Strat_Seri...
Fault_Series
Strat_Series
True
False


.. GENERATED FROM PYTHON SOURCE LINES 125-128 Setting a group as a fault also populates ``fault_relations``: a boolean matrix of which groups each fault offsets. Here, ``Fault_Series`` (row 0) affects ``Strat_Series`` (column 1), and nothing affects the fault itself: .. GENERATED FROM PYTHON SOURCE LINES 130-132 .. code-block:: Python geo_model.structural_frame.fault_relations .. rst-class:: sphx-glr-script-out .. code-block:: none array([[False, True], [False, False]]) .. GENERATED FROM PYTHON SOURCE LINES 133-141 Building structural elements and groups directly --------------------------------------------------- Importing from a CSV is only one way to get data into a model. Since a ``StructuralElement`` is just a plain data class, it can be constructed directly from arrays -- useful when adding a unit that doesn't come from a file, or when building a model up incrementally (the next tutorial does exactly this, one borehole reading at a time). A new element needs at least two surface points and one orientation somewhere in its group before a model can be computed: .. GENERATED FROM PYTHON SOURCE LINES 143-156 .. code-block:: Python new_element = gp.data.StructuralElement( name='Example_Surface', color=next(geo_model.structural_frame.color_generator), surface_points=gp.data.SurfacePointsTable.from_arrays( x=np.array([500, 1500]), y=np.array([1000, 1000]), z=np.array([600, 600]), names='Example_Surface' ), orientations=gp.data.OrientationsTable.initialize_empty() ) new_element .. raw:: html
StructuralElement:
Name:Example_Surface


.. GENERATED FROM PYTHON SOURCE LINES 157-159 A ``StructuralGroup`` is likewise just a name, a list of elements, and a relation type: .. GENERATED FROM PYTHON SOURCE LINES 161-168 .. code-block:: Python new_group = gp.data.StructuralGroup( name='Example_Series', elements=[new_element], structural_relation=gp.data.StackRelationType.ERODE ) new_group .. raw:: html
StructuralGroup:
Name:Example_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Example_Surface


.. GENERATED FROM PYTHON SOURCE LINES 169-181 Adding either of these to a live model is a matter of inserting them into the structural frame -- ``existing_group.append_element(...)`` for an element joining an existing group, or ``structural_frame.insert_group(index, group)`` for a whole new group -- both covered as part of an actual worked example in the next tutorial. This example isn't inserted here, to keep the model above unchanged for the rest of this tutorial. Visualizing input data ------------------------- With the data imported and organized into groups, it can be checked visually before computing anything. ``plot_2d`` projects the input data onto a plane along a chosen ``direction`` (``'x'``, ``'y'``, or ``'z'``, default ``'y'``): .. GENERATED FROM PYTHON SOURCE LINES 183-185 .. code-block:: Python gpv.plot_2d(geo_model, show_lith=False, show_boundaries=False) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_001.png :alt: Cell Number: mid Direction: y :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 186-187 and ``plot_3d`` shows the same data in an interactive 3D view: .. GENERATED FROM PYTHON SOURCE LINES 189-191 .. code-block:: Python gpv.plot_3d(geo_model, show_lith=False) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_002.png :alt: a01 basics :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 192-197 Computing the model ---------------------- The interpolation parameters live in ``interpolation_options``, with sensible defaults (see the Grids tutorial for what ``number_octree_levels`` specifically controls) -- change them only if you understand the implications: .. GENERATED FROM PYTHON SOURCE LINES 199-201 .. code-block:: Python geo_model.interpolation_options .. rst-class:: sphx-glr-script-out .. code-block:: none 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=6, _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=, mesh_extraction_fancy=True, evaluation_chunk_size=500000, compute_scalar=True, compute_scalar_gradient=False, verbose=False), debug=False, cache_mode=, cache_model_name='Tutorial_Basics', block_solutions_type=, sigmoid_slope=5000000) .. GENERATED FROM PYTHON SOURCE LINES 202-204 ``compute_model`` runs the interpolation and returns a ``Solutions`` object, which is also stored on the model itself as ``geo_model.solutions`` for later reference: .. GENERATED FROM PYTHON SOURCE LINES 206-209 .. code-block:: Python gp.compute_model(geo_model) geo_model.solutions .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.PYTORCH GPU enabled. Using device: cuda GPU device count: 1 Current GPU device: 0 Chunking done: 18 chunks Chunking done: 16 chunks Chunking done: 89 chunks Chunking done: 7 chunks Chunking done: 41 chunks Chunking done: 7 chunks .. raw:: html
Solutions: 6 Octree Levels, 5 DualContouringMeshes


.. GENERATED FROM PYTHON SOURCE LINES 210-214 Visualizing the result -------------------------- The computed lithology block plots the same way as the input data, by default showing a section through the middle of the model: .. GENERATED FROM PYTHON SOURCE LINES 216-218 .. code-block:: Python gpv.plot_2d(geo_model, show_data=True, cell_number="mid", direction='y') .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_003.png :alt: Cell Number: mid Direction: y :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 219-221 Each structural group has its own scalar field, selectable via ``series_n`` (its position in ``map_stack_to_surfaces``, 0-indexed) -- series 0 is the fault: .. GENERATED FROM PYTHON SOURCE LINES 223-225 .. code-block:: Python gpv.plot_2d(geo_model, series_n=0, show_data=False, show_scalar=True, show_lith=False) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_004.png :alt: Cell Number: mid Direction: y :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 226-227 and series 1 is the stratigraphy, visibly offset by the fault: .. GENERATED FROM PYTHON SOURCE LINES 229-231 .. code-block:: Python gpv.plot_2d(geo_model, series_n=1, show_data=False, show_scalar=True, show_lith=False) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_005.png :alt: Cell Number: mid Direction: y :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 232-233 The same result in 3D, with the surfaces extracted via dual contouring: .. GENERATED FROM PYTHON SOURCE LINES 235-237 .. code-block:: Python gpv.plot_3d(geo_model, show_data=False) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_006.png :alt: a01 basics :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 238-243 Adding topography -------------------- gempy supports several other grid types for different purposes -- the Grids tutorial covers all of them in depth. A quick, practical one to see here is topography, which lets a model's surfaces be intersected with real (or, as below, synthetic) terrain: .. GENERATED FROM PYTHON SOURCE LINES 245-255 .. code-block:: Python gp.set_topography_from_random( grid=geo_model.grid, fractal_dimension=1.2, d_z=np.array([350, 750]), topography_resolution=np.array([50, 50]), ) gp.compute_model(geo_model) gpv.plot_2d(geo_model, show_topography=True) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_007.png :alt: Cell Number: mid Direction: y :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Active grids: GridTypes.OCTREE|TOPOGRAPHY|NONE Setting Backend To: AvailableBackends.PYTORCH GPU enabled. Using device: cuda GPU device count: 1 Current GPU device: 0 Chunking done: 18 chunks Chunking done: 16 chunks Chunking done: 89 chunks Chunking done: 7 chunks Chunking done: 41 chunks Chunking done: 7 chunks .. GENERATED FROM PYTHON SOURCE LINES 256-258 .. code-block:: Python gpv.plot_3d(geo_model, show_lith=True, show_topography=True) .. image-sg:: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_008.png :alt: a01 basics :srcset: /tutorials/b_fundamentals/images/sphx_glr_a01_basics_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 259-265 Extracting solutions ------------------------ Beyond plotting, ``geo_model.solutions`` holds the raw building blocks of the model for further analysis or export. ``dc_meshes`` is a list of the extracted surface meshes, in the same order as the structural frame -- index ``0`` is the youngest element, the fault: .. GENERATED FROM PYTHON SOURCE LINES 267-271 .. code-block:: Python vertices = geo_model.solutions.dc_meshes[0].vertices edges = geo_model.solutions.dc_meshes[0].edges vertices.shape, edges.shape .. rst-class:: sphx-glr-script-out .. code-block:: none ((1656, 3), (3179, 3)) .. GENERATED FROM PYTHON SOURCE LINES 272-275 These vertex coordinates are in gempy's internal, rescaled coordinate system rather than the model's real-world extent. ``input_transform`` (the same transform used to normalize input data before interpolation) maps them back: .. GENERATED FROM PYTHON SOURCE LINES 277-279 .. code-block:: Python geo_model.input_transform.apply_inverse(vertices) .. rst-class:: sphx-glr-script-out .. code-block:: none array([[4.14962628e+02, 2.09349265e+01, 1.89621196e+00], [4.14810780e+02, 6.25974383e+01, 2.05674592e+00], [4.14659464e+02, 1.04260007e+02, 2.21678183e+00], ..., [1.11106592e+03, 1.97907262e+03, 7.17182139e+02], [1.13428542e+03, 1.93741013e+03, 7.40676218e+02], [1.13433491e+03, 1.97907265e+03, 7.40615851e+02]], shape=(1656, 3)) .. GENERATED FROM PYTHON SOURCE LINES 280-283 ``raw_arrays`` holds the underlying arrays directly -- the lithology block (``lith_block``), for instance, comes back as a flat array that needs reshaping to the grid's actual resolution to index into as a volume: .. GENERATED FROM PYTHON SOURCE LINES 285-288 .. code-block:: Python lith_block = geo_model.solutions.raw_arrays.lith_block lith_block.shape .. rst-class:: sphx-glr-script-out .. code-block:: none (2359296,) .. GENERATED FROM PYTHON SOURCE LINES 289-291 .. code-block:: Python lith_block.reshape(geo_model.grid.regular_grid.resolution).shape .. rst-class:: sphx-glr-script-out .. code-block:: none (192, 192, 64) .. GENERATED FROM PYTHON SOURCE LINES 292-296 Saving and loading a model ------------------------------ A ``GeoModel`` can be saved to a single file and reloaded later, without needing to redo the setup above: .. GENERATED FROM PYTHON SOURCE LINES 298-300 .. code-block:: Python gp.save_model(geo_model, path='tutorial_basics_model.gempy') .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/buildAgent/work/3a8738c25f60c3c9/gempy/modules/serialization/save_load.py:33: UserWarning: This function is still in development. It may not work as expected. warnings.warn("This function is still in development. It may not work as expected.") 'tutorial_basics_model.gempy' .. GENERATED FROM PYTHON SOURCE LINES 301-304 .. code-block:: Python reloaded_model = gp.load_model('tutorial_basics_model.gempy') reloaded_model.structural_frame .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/buildAgent/work/3a8738c25f60c3c9/gempy/modules/serialization/save_load.py:118: UserWarning: This function is still in development. It may not work as expected. warnings.warn("This function is still in development. It may not work as expected.") .. raw:: html
Structural Groups: StructuralGroup:
Name:Fault_Series
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Main_Fault

StructuralGroup:
Name:Strat_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Sandstone_2

StructuralElement:
Name:Siltstone

StructuralElement:
Name:Shale

StructuralElement:
Name:Sandstone_1
Fault Relations:
Fault_Seri...Strat_Seri...
Fault_Series
Strat_Series
True
False


.. GENERATED FROM PYTHON SOURCE LINES 305-309 .. note:: Model serialization is still marked as under active development in gempy (you'll see a ``UserWarning`` when saving/loading) -- it works, but the format may still change in a future release. .. GENERATED FROM PYTHON SOURCE LINES 309-311 .. code-block:: Python # sphinx_gallery_thumbnail_number = -3 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 26.598 seconds) .. _sphx_glr_download_tutorials_b_fundamentals_a01_basics.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: a01_basics.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: a01_basics.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: a01_basics.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_