.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "Tutorials\01_POC_generate-model.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_01_POC_generate-model.py: 01 - Geological model creation and gravity simulation ================================================ The following tutorial will step-by-step lead you through an example workflow on creating a GemPy model from interface and orientation data, assigning densities to geological units, and model their gravity response. .. GENERATED FROM PYTHON SOURCE LINES 9-15 Create the base Proof-of-Concept Model -------------------------------------- Based on a seismic section from the NAGRA report `NAGRA NAB 14-17 `_ [1], we extracted interface and orientation points for lithological units and faults. The lithological units comprise the permo-carboniferous filling (divided in three stages based on the report results), Mesozoic, Tertiary, and Quaternary strata, as well as the Palaeozoic crystalline basement rocks. .. GENERATED FROM PYTHON SOURCE LINES 15-33 .. code-block:: default import warnings warnings.filterwarnings("ignore") # Importing GemPy import gempy as gp from gempy.plot import visualization_2d as vv # Importing auxilary libraries import numpy as np import pandas as pd import matplotlib.pyplot as plt plt.style.use('seaborn-talk') # What GemPy version was used print(f"Code run with GemPy version: {gp.__version__}") .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Code run with GemPy version: 2.2.9 .. GENERATED FROM PYTHON SOURCE LINES 34-37 Initialize the model -------------------- We start with modelling the trough by generating a gempy model object. This will use interface points and orientations, which we previously stored in a `.csv` file. .. GENERATED FROM PYTHON SOURCE LINES 37-52 .. code-block:: default # Fix random number seed to get the same topography np.random.seed(333) # Import data # Create a model instance geo_model = gp.create_model('POC_model') # Initialize the model, set dimension and load interface and orientation data gp.init_data(geo_model, [0, 28000., 0, 14000., -6500, 1000.], [100, 50, 60], path_i = '../../data/GemPy/line82_interfaces_wo_middle_MC.csv', path_o = '../../data/GemPy/line82_foliations_wo_middle_MC.csv') geo_model.set_topography(source='random', d_z=np.array([300,1000])) gp.plot_2d(geo_model, show_data=True, show_topography=True) .. image:: /Tutorials/images/sphx_glr_01_POC_generate-model_001.png :alt: Cell Number: mid Direction: y :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Active grids: ['regular'] Active grids: ['regular' 'topography'] .. GENERATED FROM PYTHON SOURCE LINES 53-56 Adding information to the model ------------------------------- Only loading interface and orientation points is not enough. First, let's assign colors to the different model units, e.g. for coloring faults similarly. .. GENERATED FROM PYTHON SOURCE LINES 56-72 .. code-block:: default col_dict = {'basement': '#c7848f', 'Lower-filling': '#a5d490', 'Upper-filling': '#cfc199', 'Unconformity': '#725c9a', 'Orange': '#ff792b', 'Pink': '#e588f3', 'Tertiary': '#dbdbac', 'Fault2': '#015482', 'Fault5': '#015482', 'Fault6': '#015482', 'Thrust1_south': '#5DA629', 'Thrust2_south': '#5DA629'} geo_model.surfaces.colors.change_colors(col_dict) geo_model.surfaces .. raw:: html
surface series order_surfaces color id
0 Fault2 Default series 1 #015482 1
1 Fault5 Default series 2 #015482 2
2 Fault6 Default series 3 #015482 3
3 Lower-filling Default series 4 #a5d490 4
4 Upper-filling Default series 5 #cfc199 5
5 Unconformity Default series 6 #725c9a 6
6 Tertiary Default series 7 #dbdbac 7
7 Orange Default series 8 #ff792b 8
8 Pink Default series 9 #e588f3 9
9 Thrust1_south Default series 10 #5DA629 10
10 Thrust2_south Default series 11 #5DA629 11
11 basement Basement 1 #c7848f 12


.. GENERATED FROM PYTHON SOURCE LINES 73-77 Model Characteristics --------------------- Main features of the model is the asymetric graben system, with the major fault (denoted with **A**), and the graben fill, which is not present beyond the graben shoulders. This, as well as the stop of major faults beneath the mesozoic units (blue units) are important considerations for the modelling process. These could be caught, for instance, in likelihood functions if we model the PCT as a Bayesian inference problem. .. GENERATED FROM PYTHON SOURCE LINES 77-93 .. code-block:: default # Assign formations to series gp.map_series_to_surfaces(geo_model, {"Thrust1_series": 'Thrust1_south', "Thrust2_series": 'Thrust2_south', "Fault2_series": 'Fault2', "Fault5_series": 'Fault5', "Fault6_series": 'Fault6', "Post_tectonic_series": ('Tertiary', 'Pink', 'Orange'), "Detachement": 'Unconformity', "Syn_tectonic_series2": 'Upper-filling', #"Syn_tectonic_series1": 'Middle-filling', "Pre_tectonic_series": 'Lower-filling'}, remove_unused_series=True) geo_model.surfaces .. raw:: html
surface series order_surfaces color id
9 Thrust1_south Thrust1_series 1 #5DA629 1
10 Thrust2_south Thrust2_series 1 #5DA629 2
0 Fault2 Fault2_series 1 #015482 3
1 Fault5 Fault5_series 1 #015482 4
2 Fault6 Fault6_series 1 #015482 5
6 Tertiary Post_tectonic_series 1 #dbdbac 6
7 Orange Post_tectonic_series 2 #ff792b 7
8 Pink Post_tectonic_series 3 #e588f3 8
5 Unconformity Detachement 1 #725c9a 9
4 Upper-filling Syn_tectonic_series2 1 #cfc199 10
3 Lower-filling Pre_tectonic_series 1 #a5d490 11
11 basement Basement 1 #c7848f 12


.. GENERATED FROM PYTHON SOURCE LINES 94-96 After assigning units to stacks or series, we have so define which of those series is a fault. Here, we see that it is usually important to assign each fault its own series, as faults may have very different scalar fields (in which the fault surfaces are interpolated). .. GENERATED FROM PYTHON SOURCE LINES 96-101 .. code-block:: default geo_model.set_is_fault(['Thrust1_series', 'Thrust2_series', 'Fault2_series', 'Fault5_series', 'Fault6_series'], change_color=False) .. raw:: html
order_series BottomRelation isActive isFault isFinite
Thrust1_series 1 Fault True True False
Thrust2_series 2 Fault True True False
Fault2_series 3 Fault True True False
Fault5_series 4 Fault True True False
Fault6_series 5 Fault True True False
Post_tectonic_series 6 Erosion True False False
Detachement 7 Erosion True False False
Syn_tectonic_series2 8 Erosion True False False
Pre_tectonic_series 9 Erosion True False False
Basement 10 Erosion False False False


.. GENERATED FROM PYTHON SOURCE LINES 102-103 Further we have to set bottom relations, if a series is **not** erosive. For instance, the Units in the Graben are most likely onlapping units. .. GENERATED FROM PYTHON SOURCE LINES 103-107 .. code-block:: default geo_model.set_bottom_relation(series=['Post_tectonic_series', 'Pre_tectonic_series', 'Syn_tectonic_series2'], bottom_relation='Onlap') #, .. raw:: html
order_series BottomRelation isActive isFault isFinite
Thrust1_series 1 Fault True True False
Thrust2_series 2 Fault True True False
Fault2_series 3 Fault True True False
Fault5_series 4 Fault True True False
Fault6_series 5 Fault True True False
Post_tectonic_series 6 Onlap True False False
Detachement 7 Erosion True False False
Syn_tectonic_series2 8 Onlap True False False
Pre_tectonic_series 9 Onlap True False False
Basement 10 Erosion False False False


.. GENERATED FROM PYTHON SOURCE LINES 108-110 The following table shows the fault relations, i.e. which unit (or fault) is affected by a fault. If the respective entry in the table is set to `True`, the fault on the left displaces the unit (or fault) in a respective column. .. GENERATED FROM PYTHON SOURCE LINES 110-113 .. code-block:: default geo_model.faults.faults_relations_df .. raw:: html
Thrust1_series Thrust2_series Fault2_series Fault5_series Fault6_series Post_tectonic_series Detachement Syn_tectonic_series2 Pre_tectonic_series Basement
Thrust1_series False False False False False True True True True True
Thrust2_series False False False False False True True True True True
Fault2_series False False False False False True True True True True
Fault5_series False False False False False True True True True True
Fault6_series False False False False False True True True True True
Post_tectonic_series False False False False False False False False False False
Detachement False False False False False False False False False False
Syn_tectonic_series2 False False False False False False False False False False
Pre_tectonic_series False False False False False False False False False False
Basement False False False False False False False False False False


.. GENERATED FROM PYTHON SOURCE LINES 114-115 Per default, faults displace all lithological units. However, the normal faults of the graben do not affect the younger units, so we define a boolean matrix, which sets the fault relations correctly. .. GENERATED FROM PYTHON SOURCE LINES 115-129 .. code-block:: default fr = np.array([[False, True, False, False, False, True, False, False, False, False], [False, False, False, False, False, True, False, False, False, False], [False, False, False, False, False, False, True, True, True, True], [False, False, False, False, False, False, True, True, True, True], [False, False, False, False, False, False, True, True, True, True], [False, False, False, False, False, False, False, False, False, False], [False, False, False, False, False, False, False, False, False, False], [False, False, False, False, False, False, False, False, False, False], [False, False, False, False, False, False, False, False, False, False], [False, False, False, False, False, False, False, False, False, False]]) geo_model.set_fault_relation(fr) .. raw:: html
Thrust1_series Thrust2_series Fault2_series Fault5_series Fault6_series Post_tectonic_series Detachement Syn_tectonic_series2 Pre_tectonic_series Basement
Thrust1_series False True False False False True False False False False
Thrust2_series False False False False False True False False False False
Fault2_series False False False False False False True True True True
Fault5_series False False False False False False True True True True
Fault6_series False False False False False False True True True True
Post_tectonic_series False False False False False False False False False False
Detachement False False False False False False False False False False
Syn_tectonic_series2 False False False False False False False False False False
Pre_tectonic_series False False False False False False False False False False
Basement False False False False False False False False False False


.. GENERATED FROM PYTHON SOURCE LINES 130-133 Creating the model ------------------ Now that we set the parameters and fault relations, it is time to start the modeling process: .. GENERATED FROM PYTHON SOURCE LINES 133-148 .. code-block:: default # decrease the kriging range geo_model.modify_kriging_parameters('range', 20000.) geo_model.modify_kriging_parameters('$C_o$', 2e5) # Set the interpolator function gp.set_interpolator(geo_model, compile_theano=True, theano_optimizer='fast_compile', verbose=[], update_kriging=False) # Compute the model sol = gp.compute_model(geo_model) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Compiling theano function... Level of Optimization: fast_compile Device: cpu Precision: float64 Number of faults: 5 Compilation Done! Kriging values: values range 20000 $C_o$ 200000 drift equations [3, 3, 3, 3, 3, 3, 3, 3, 3, 3] .. GENERATED FROM PYTHON SOURCE LINES 149-150 Saving the model is straight forward. It can optionally also be compressed in a zip archive, or be _pickled_. An example on how to save a model is shown next. There, we give the saving path and the model name. .. GENERATED FROM PYTHON SOURCE LINES 150-155 .. code-block:: default geo_model.save_model(name='POC_PCT_model', path='../../models/2021-06-04_POC_base_model') .. rst-class:: sphx-glr-script-out Out: .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 156-157 Let's have a look how the created model looks like: .. GENERATED FROM PYTHON SOURCE LINES 157-160 .. code-block:: default gp.plot_2d(geo_model, cell_number=25, direction='y', show_data=False, show_topography=False, show_lith=True, show_results=True, show_boundaries=True) .. image:: /Tutorials/images/sphx_glr_01_POC_generate-model_002.png :alt: Cell Number: 25 Direction: y :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 161-163 The displayed lithologies are stored in a regular grid, which we defined at the beginning. We can store this grid, containing the lithology IDs in space for further postprocessing .. GENERATED FROM PYTHON SOURCE LINES 163-166 .. code-block:: default np.save('../../models/POC_base_model_lith_blocks.npy', np.round(geo_model.solutions.lith_block,0).astype('int')) .. GENERATED FROM PYTHON SOURCE LINES 167-171 Simulate Gravity ================ Using the now generated POC-model, we simulate its gravity at different locations. These locations will be treated as observations later on in the workflow. In a first step, we distribute 15 points randomly across the topography of our model. Those will be the station locations, where we pick up the gravity signal of the POC-model. .. GENERATED FROM PYTHON SOURCE LINES 171-182 .. code-block:: default # distribute stations import random np.random.seed(58) station_indices = np.random.randint(0, high=4999, size=15) station_coordinates = geo_model._grid.topography.values[station_indices, :] cs = plt.scatter(station_coordinates[:,0], station_coordinates[:,1], c=station_coordinates[:,2], cmap='viridis') plt.colorbar(cs) .. image:: /Tutorials/images/sphx_glr_01_POC_generate-model_003.png :alt: 01 POC generate model :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 183-184 Next, we create centered grids around each station. The centered grid here has 10 cells in x- and y-direction, and extends 15 cells down in the z-direction. .. GENERATED FROM PYTHON SOURCE LINES 184-190 .. code-block:: default from gempy.assets.geophysics import GravityPreprocessing geo_model.set_centered_grid(station_coordinates, resolution = [10, 10, 15], radius=6000) g = GravityPreprocessing(geo_model.grid.centered_grid) tz = g.set_tz_kernel() .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Active grids: ['regular' 'topography' 'centered'] .. GENERATED FROM PYTHON SOURCE LINES 191-193 The gravity response cannot be modeled without assigning a density to the model units. Theoretically, one could also assign different petrophyiscal properties here. They will be added as separate columns to the surfaces dataframe. .. GENERATED FROM PYTHON SOURCE LINES 193-198 .. code-block:: default densities = [0, 0, 0, 0, 0, 2.466, 2.61, 2.53, 2.61, 2.47, 2.55, 2.67] geo_model.add_surface_values(densities, ['density']) .. raw:: html
surface series order_surfaces color id density
9 Thrust1_south Thrust1_series 1 #5DA629 1 0.000000
10 Thrust2_south Thrust2_series 1 #5DA629 2 0.000000
0 Fault2 Fault2_series 1 #015482 3 0.000000
1 Fault5 Fault5_series 1 #015482 4 0.000000
2 Fault6 Fault6_series 1 #015482 5 0.000000
6 Tertiary Post_tectonic_series 1 #dbdbac 6 2.466000
8 Pink Post_tectonic_series 2 #e588f3 7 2.610000
7 Orange Post_tectonic_series 3 #ff792b 8 2.530000
5 Unconformity Detachement 1 #725c9a 9 2.610000
4 Upper-filling Syn_tectonic_series2 1 #cfc199 10 2.470000
3 Lower-filling Pre_tectonic_series 1 #a5d490 11 2.550000
11 basement Basement 1 #c7848f 12 2.670000


.. GENERATED FROM PYTHON SOURCE LINES 199-200 Modeling the lithology on all grids (regular, topography, centered) can get time consuming. So here, we only activate the centered grid to catch the gravity response. .. GENERATED FROM PYTHON SOURCE LINES 200-213 .. code-block:: default geo_model.set_active_grid('centered', reset=True) gp.set_interpolator(geo_model, output=['gravity'], theano_optimizer='fast_run', update_kriging=False) sol = gp.compute_model(geo_model) # reshape solved gravity and add coordinates grav = sol.fw_gravity grav1 = grav.reshape(len(grav),1) station_forw_grav = np.round(np.append(station_coordinates, grav1, axis=1),4) # make everything into a dataframe df_stations = pd.DataFrame(station_forw_grav, columns=["X", "Y", "Z", "grav"]) # round X Y and Z to 2 decimals df_stations[['X','Y','Z']] = np.around(df_stations[['X','Y','Z']], 2) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Active grids: ['centered'] Compiling theano function... Level of Optimization: fast_run Device: cpu Precision: float64 Number of faults: 5 Compilation Done! Kriging values: values range 20000 $C_o$ 200000 drift equations [3, 3, 3, 3, 3, 3, 3, 3, 3, 3] .. GENERATED FROM PYTHON SOURCE LINES 214-215 and finally, we save the modeled gravity to be used as observations later on: .. GENERATED FROM PYTHON SOURCE LINES 215-218 .. code-block:: default df_stations.to_csv('../../data/Data_for_MC/20210322_forw_grav_seed58.csv', index=False) .. GENERATED FROM PYTHON SOURCE LINES 219-221 References ---------- [1] Naef, H., & Madritsch, H. (2014). Tektonische Karte des Nordschweizer Permokarbontrogs: Aktualisierung basierend auf 2D-Seismik und Schweredaten. Nagra Arbeitsbericht (NAB 14-17). Wettingen: Nagra. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 59.426 seconds) .. _sphx_glr_download_Tutorials_01_POC_generate-model.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: 01_POC_generate-model.py <01_POC_generate-model.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01_POC_generate-model.ipynb <01_POC_generate-model.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_