.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/ch4-Topology/ch4-1-Topology.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_ch4-Topology_ch4-1-Topology.py: Chapter 4: Analyzing Geomodel Topology ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 6-17 .. code-block:: Python import gempy as gp import gempy_viewer as gpv from gempy_viewer.modules.plot_2d.visualization_2d import Plot2D from gempy_plugins.topology_analysis import topology as tp import os import warnings warnings.filterwarnings("ignore") .. GENERATED FROM PYTHON SOURCE LINES 18-26 Load example Model ^^^^^^^^^^^^^^^^^^ First let's set up a very simple example model. For that we initialize the geo_data object with the correct model extent and the resolution we like. Then we load our data points from csv files and set the series and order the formations (stratigraphic pile). .. GENERATED FROM PYTHON SOURCE LINES 28-55 .. code-block:: Python data_path = os.path.abspath('../../') geo_model: gp.data.GeoModel = gp.create_geomodel( project_name='Model_Tutorial6', extent= [0, 3000, 0, 20, 0, 2000], resolution=[50, 10, 67], refinement=1, # * For this model is better not to use octrees because we want to see what is happening in the scalar fields importer_helper=gp.data.ImporterHelper( path_to_orientations=data_path + "/data/input_data/tut_chapter6/ch6_data_fol.csv", path_to_surface_points=data_path + "/data/input_data/tut_chapter6/ch6_data_interf.csv", ) ) gp.map_stack_to_surfaces( gempy_model=geo_model, mapping_object= { "fault": "Fault", "Rest": ('Layer 2', 'Layer 3', 'Layer 4', 'Layer 5') } ) gp.set_is_fault(geo_model, ['fault']) geo_model.interpolation_options.mesh_extraction = False sol = gp.compute_model(geo_model) .. rst-class:: sphx-glr-script-out .. code-block:: none Setting Backend To: AvailableBackends.numpy .. GENERATED FROM PYTHON SOURCE LINES 56-59 .. code-block:: Python gpv.plot_2d(geo_model, cell_number=[5]) .. image-sg:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_001.png :alt: Cell Number: 5 Direction: y :srcset: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 60-76 Analyzing Topology ^^^^^^^^^^^^^^^^^^ GemPy sports in-built functionality to analyze the topology of its models. All we need for this is our geo_data object, lithology block and the fault block. We input those into *gp.topology_compute* and get several useful outputs: - an adjacency graph **G**, representing the topological relationships of the model - the **centroids** of the all the unique topological regions in the model (x,y,z coordinates of their center) - a list of all the unique labels (labels_unique) - two look-up-tables from the lithology id's to the node labels, and vice versa .. GENERATED FROM PYTHON SOURCE LINES 78-81 .. code-block:: Python edges, centroids = tp.compute_topology(geo_model) .. GENERATED FROM PYTHON SOURCE LINES 82-87 The first output of the topology function is the ``set`` of edges representing topology relationships between unique geobodies of the block model. An edge is represented by a ``tuple`` of two ``int`` geobody (or node) labels: .. GENERATED FROM PYTHON SOURCE LINES 89-92 .. code-block:: Python edges .. rst-class:: sphx-glr-script-out .. code-block:: none {(9, 10), (4, 10), (1, 2), (3, 4), (1, 8), (3, 10), (2, 3), (2, 9), (1, 7), (4, 5), (3, 9), (5, 10), (6, 7), (8, 9), (1, 6), (7, 8), (2, 8)} .. GENERATED FROM PYTHON SOURCE LINES 93-97 The second output is the centroids ``dict``, mapping the unique geobody id's (graph node id's) to the geobody centroid position in grid coordinates: .. GENERATED FROM PYTHON SOURCE LINES 99-102 .. code-block:: Python centroids .. rst-class:: sphx-glr-script-out .. code-block:: none {1: array([35.27893175, 4.5 , 50.19485658]), 2: array([36.46666667, 4.5 , 29.14444444]), 3: array([37.59756098, 4.5 , 21.62195122]), 4: array([38.84563758, 4.5 , 14.00671141]), 5: array([39.09550562, 4.5 , 5.37640449]), 6: array([ 9.79081633, 4.5 , 60.10204082]), 7: array([10.17687075, 4.5 , 51.02721088]), 8: array([11.37804878, 4.5 , 43.47560976]), 9: array([12.51098901, 4.5 , 35.90659341]), 10: array([13.659857 , 4.5 , 15.34320735])} .. GENERATED FROM PYTHON SOURCE LINES 103-106 After computing the model topology, we can overlay the topology graph over a model section: .. GENERATED FROM PYTHON SOURCE LINES 109-115 Visualizing topology ~~~~~~~~~~~~~~~~~~~~ 2-D Visualization of the Topology Graph ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 118-124 .. code-block:: Python gpv.plot_topology( regular_grid=geo_model.grid.regular_grid, edges=edges, centroids=centroids ) .. image-sg:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_002.png :alt: ch4 1 Topology :srcset: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 125-133 .. code-block:: Python plot_2d: Plot2D = gpv.plot_2d(geo_model, cell_number=[5], show=False) gpv.plot_topology( regular_grid=geo_model.grid.regular_grid, edges=edges, centroids=centroids, ax=plot_2d.axes[0] ) .. image-sg:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_003.png :alt: Cell Number: 5 Direction: y :srcset: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 134-140 Adjacency Matrix ~~~~~~~~~~~~~~~~ Another way to encode and visualize the geomodel topology is using an adjacency graph: .. GENERATED FROM PYTHON SOURCE LINES 142-145 .. code-block:: Python M = tp.get_adjacency_matrix(geo_model, edges, centroids) print(M) .. rst-class:: sphx-glr-script-out .. code-block:: none [[False True False False False True True True False False] [ True False True False False False False True True False] [False True False True False False False False True True] [False False True False True False False False False True] [False False False True False False False False False True] [ True False False False False False True False False False] [ True False False False False True False True False False] [ True True False False False False True False True False] [False True True False False False False True False True] [False False True True True False False False True False]] .. GENERATED FROM PYTHON SOURCE LINES 146-149 .. code-block:: Python tp.plot_adjacency_matrix(geo_model, M) .. image-sg:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_004.png :alt: Topology Adjacency Matrix :srcset: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 150-153 Look-up tables ~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 156-161 The ``topology`` asset provides several look-up tables to work with the unique geobody topology id's. Mapping node id's back to lithology / surface id's: .. GENERATED FROM PYTHON SOURCE LINES 163-167 .. code-block:: Python lith_lot = tp.get_lot_node_to_lith_id(geo_model, centroids) lith_lot .. rst-class:: sphx-glr-script-out .. code-block:: none {1: 2, 2: 3, 3: 4, 4: 5, 5: 6, 6: 2, 7: 3, 8: 4, 9: 5, 10: 6} .. GENERATED FROM PYTHON SOURCE LINES 168-170 Figuring out which nodes are in which fault block: .. GENERATED FROM PYTHON SOURCE LINES 172-176 .. code-block:: Python fault_lot = tp.get_lot_node_to_fault_block(geo_model, centroids) fault_lot .. rst-class:: sphx-glr-script-out .. code-block:: none {1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 1, 7: 1, 8: 1, 9: 1, 10: 1} .. GENERATED FROM PYTHON SOURCE LINES 177-180 We can also easily map the lithology id to the corresponding topology id's: .. GENERATED FROM PYTHON SOURCE LINES 182-185 .. code-block:: Python tp.get_lot_lith_to_node_id(lith_lot) .. rst-class:: sphx-glr-script-out .. code-block:: none {2: [1, 6], 3: [2, 7], 4: [3, 8], 5: [4, 9], 6: [5, 10]} .. GENERATED FROM PYTHON SOURCE LINES 186-189 Detailed node labeling ~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 191-192 sphinx_gallery_thumbnail_number = 4 .. GENERATED FROM PYTHON SOURCE LINES 192-193 .. code-block:: Python dedges, dcentroids = tp.get_detailed_labels(geo_model, edges, centroids) .. GENERATED FROM PYTHON SOURCE LINES 194-202 .. code-block:: Python plot_2d: Plot2D = gpv.plot_2d(geo_model, cell_number=[5], show=False) gpv.plot_topology( regular_grid=geo_model.grid.regular_grid, edges=dedges, centroids=dcentroids, ax=plot_2d.axes[0] ) .. image-sg:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_005.png :alt: Cell Number: 5 Direction: y :srcset: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-205 .. code-block:: Python dedges .. rst-class:: sphx-glr-script-out .. code-block:: none {('4_0', '6_1'), ('2_0', '3_1'), ('5_0', '6_0'), ('3_0', '4_0'), ('3_1', '4_1'), ('3_0', '5_1'), ('2_0', '2_1'), ('2_0', '4_1'), ('3_0', '4_1'), ('4_0', '5_1'), ('5_0', '6_1'), ('4_0', '5_0'), ('6_0', '6_1'), ('2_1', '3_1'), ('4_1', '5_1'), ('5_1', '6_1'), ('2_0', '3_0')} .. GENERATED FROM PYTHON SOURCE LINES 206-209 .. code-block:: Python dcentroids .. rst-class:: sphx-glr-script-out .. code-block:: none {'2_0': array([35.27893175, 4.5 , 50.19485658]), '3_0': array([36.46666667, 4.5 , 29.14444444]), '4_0': array([37.59756098, 4.5 , 21.62195122]), '5_0': array([38.84563758, 4.5 , 14.00671141]), '6_0': array([39.09550562, 4.5 , 5.37640449]), '2_1': array([ 9.79081633, 4.5 , 60.10204082]), '3_1': array([10.17687075, 4.5 , 51.02721088]), '4_1': array([11.37804878, 4.5 , 43.47560976]), '5_1': array([12.51098901, 4.5 , 35.90659341]), '6_1': array([13.659857 , 4.5 , 15.34320735])} .. GENERATED FROM PYTHON SOURCE LINES 210-213 Checking adjacency ~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 216-220 So lets say we want to check if the purple layer (id 5) is connected across the fault to the yellow layer (id 3). For this we can make easy use of the detailed labeling and the ``check_adjacency`` function: .. GENERATED FROM PYTHON SOURCE LINES 222-225 .. code-block:: Python tp.check_adjacency(dedges, "5_1", "3_0") .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 226-229 We can also check all geobodies that are adjacent to the purple layer (id 5) on the left side of the fault (fault id 1): .. GENERATED FROM PYTHON SOURCE LINES 231-233 .. code-block:: Python tp.get_adjacencies(dedges, "5_1") .. rst-class:: sphx-glr-script-out .. code-block:: none {'4_1', '3_0', '6_1', '4_0'} .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.739 seconds) .. _sphx_glr_download_tutorials_ch4-Topology_ch4-1-Topology.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ch4-1-Topology.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ch4-1-Topology.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_