.. 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 Click :ref:`here ` 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-16 .. code-block:: python3 import gempy as gp from gempy.assets import topology as tp import numpy as np import matplotlib.pyplot as plt import warnings warnings.filterwarnings("ignore") .. GENERATED FROM PYTHON SOURCE LINES 17-25 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 27-46 .. code-block:: python3 geo_model = gp.create_model("Model_Tutorial6") data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/' gp.init_data( geo_model, [0, 3000, 0, 20, 0, 2000], [50, 10, 67], path_i=data_path+"data/input_data/tut_chapter6/ch6_data_interf.csv", path_o=data_path+"data/input_data/tut_chapter6/ch6_data_fol.csv" ) gp.map_stack_to_surfaces( geo_model, { "fault": "Fault", "Rest": ('Layer 2', 'Layer 3', 'Layer 4', 'Layer 5') } ) geo_model.set_is_fault(["fault"]); gp.set_interpolator(geo_model) sol = gp.compute_model(geo_model, compute_mesh=True) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Active grids: ['regular'] Fault colors changed. If you do not like this behavior, set change_color to False. Setting kriging parameters to their default values. Compiling theano function... Level of Optimization: fast_compile Device: cpu Precision: float64 Number of faults: 1 Compilation Done! Kriging values: values range 3605.61 $C_o$ 309533.33 drift equations [3, 3, 3] .. GENERATED FROM PYTHON SOURCE LINES 47-50 .. code-block:: python3 gp.plot_2d(geo_model, cell_number=[5]) .. image:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_001.png :alt: Cell Number: 5 Direction: y :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 51-67 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 69-72 .. code-block:: python3 edges, centroids = tp.compute_topology(geo_model) .. GENERATED FROM PYTHON SOURCE LINES 73-78 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 80-83 .. code-block:: python3 edges .. rst-class:: sphx-glr-script-out 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 84-88 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 90-93 .. code-block:: python3 centroids .. rst-class:: sphx-glr-script-out 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 94-97 After computing the model topology, we can overlay the topology graph over a model section: .. GENERATED FROM PYTHON SOURCE LINES 100-106 Visualizing topology ~~~~~~~~~~~~~~~~~~~~ 2-D Visualization of the Topology Graph ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 108-111 .. code-block:: python3 gp.plot.plot_topology(geo_model, edges, centroids) plt.show() .. image:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_002.png :alt: ch4 1 Topology :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 112-116 .. code-block:: python3 gp.plot_2d(geo_model, cell_number=[5], show=False) gp.plot.plot_topology(geo_model, edges, centroids, scale=True) plt.show() .. image:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_003.png :alt: Cell Number: 5 Direction: y :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 117-123 Adjacency Matrix ~~~~~~~~~~~~~~~~ Another way to encode and visualize the geomodel topology is using an adjacency graph: .. GENERATED FROM PYTHON SOURCE LINES 125-128 .. code-block:: python3 M = tp.get_adjacency_matrix(geo_model, edges, centroids) print(M) .. rst-class:: sphx-glr-script-out 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 129-132 .. code-block:: python3 tp.plot_adjacency_matrix(geo_model, M) .. image:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_004.png :alt: Topology Adjacency Matrix :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 133-136 3-D Visualization of the Topology Graph ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 139-142 You can also plot the topology in 3-D using GemPy's 3-D visualization toolkit powered by ``pyvista``: .. GENERATED FROM PYTHON SOURCE LINES 144-150 .. code-block:: python3 from gempy.plot._vista import Vista gpv = Vista(geo_model) gpv.plot_topology(edges, centroids) gpv.show() .. image:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_005.png :alt: ch4 1 Topology :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none surface Fault #527682 Layer 2 #9f0052 Layer 3 #ffbe00 Layer 4 #728f02 Layer 5 #443988 basement #ff3f20 Name: color, dtype: object .. GENERATED FROM PYTHON SOURCE LINES 151-154 Look-up tables ~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 157-162 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 164-168 .. code-block:: python3 lith_lot = tp.get_lot_node_to_lith_id(geo_model, centroids) lith_lot .. rst-class:: sphx-glr-script-out 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 169-171 Figuring out which nodes are in which fault block: .. GENERATED FROM PYTHON SOURCE LINES 173-177 .. code-block:: python3 fault_lot = tp.get_lot_node_to_fault_block(geo_model, centroids) fault_lot .. rst-class:: sphx-glr-script-out 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 178-181 We can also easily map the lithology id to the corresponding topology id's: .. GENERATED FROM PYTHON SOURCE LINES 183-186 .. code-block:: python3 tp.get_lot_lith_to_node_id(lith_lot) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none {2: [1, 6], 3: [2, 7], 4: [3, 8], 5: [4, 9], 6: [5, 10]} .. GENERATED FROM PYTHON SOURCE LINES 187-190 Detailed node labeling ~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 192-193 sphinx_gallery_thumbnail_number = 4 .. GENERATED FROM PYTHON SOURCE LINES 193-198 .. code-block:: python3 dedges, dcentroids = tp.get_detailed_labels(geo_model, edges, centroids) gp.plot_2d(geo_model, cell_number=[5], show=False) gp.plot.plot_topology(geo_model, dedges, dcentroids, scale=True) plt.show() .. image:: /tutorials/ch4-Topology/images/sphx_glr_ch4-1-Topology_006.png :alt: Cell Number: 5 Direction: y :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 199-201 .. code-block:: python3 dedges .. rst-class:: sphx-glr-script-out Out: .. code-block:: none {('3_0', '4_0'), ('3_1', '4_1'), ('2_0', '2_1'), ('4_1', '5_1'), ('6_0', '6_1'), ('4_0', '6_1'), ('4_0', '5_1'), ('2_0', '4_1'), ('2_0', '3_1'), ('2_0', '3_0'), ('5_0', '6_1'), ('3_0', '5_1'), ('5_0', '6_0'), ('2_1', '3_1'), ('3_0', '4_1'), ('5_1', '6_1'), ('4_0', '5_0')} .. GENERATED FROM PYTHON SOURCE LINES 202-205 .. code-block:: python3 dcentroids .. rst-class:: sphx-glr-script-out 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 206-209 Checking adjacency ~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 212-216 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 218-221 .. code-block:: python3 tp.check_adjacency(dedges, "5_1", "3_0") .. rst-class:: sphx-glr-script-out Out: .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 222-225 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 227-229 .. code-block:: python3 tp.get_adjacencies(dedges, "5_1") .. rst-class:: sphx-glr-script-out Out: .. code-block:: none {'6_1', '4_1', '4_0', '3_0'} .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 5.379 seconds) .. _sphx_glr_download_tutorials_ch4-Topology_ch4-1-Topology.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ch4-1-Topology.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ch4-1-Topology.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_