.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/models/GemPy-I.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_gallery_models_GemPy-I.py: GemPy-I: *Simple Fault Model* ============================= This example uses `GemPy `_ to create a geological model as input to emg3d, utilizing `discretize `_. Having it in discretize allows us also to plot it with `PyVista `_. The starting point is the *simple_fault_model* as used in `Chapter 1.1 `_ of the GemPy documentation. It is a nice, made-up model of a folded structure with a fault. Here we slightly modify it (convert it into a shallow marine setting), and create a resisistivity model out of the lithological model. The result is what is referred to in other examples as model `GemPy-I`, a synthetic, shallow-marine resistivity model consisting of a folded structure with a fault. It is one of a few models created to be used in other examples. .. note:: The original model (*simple_fault_model*) hosted on https://github.com/cgre-aachen/gempy_data is released under the `LGPL-3.0 License `_. .. GENERATED FROM PYTHON SOURCE LINES 28-40 .. code-block:: default import os import pooch import emg3d import matplotlib.pyplot as plt from matplotlib.colors import LogNorm plt.style.use('bmh') # Adjust this path to a folder of your choice. data_path = os.path.join('..', 'download', '') .. GENERATED FROM PYTHON SOURCE LINES 41-45 Fetch the model --------------- Retrieve and load the pre-computed resistivity model. .. GENERATED FROM PYTHON SOURCE LINES 45-57 .. code-block:: default fname = "GemPy-I.h5" pooch.retrieve( 'https://raw.github.com/emsig/data/2021-05-21/emg3d/models/'+fname, '06f522a69c94dc02ca3da0ea4ca7b60f7a9c764cdcbf6699ef4155621d70b3bb', fname=fname, path=data_path, ) fmodel = emg3d.load(data_path + fname)['model'] fgrid = fmodel.grid .. rst-class:: sphx-glr-script-out .. code-block:: none Data loaded from «/home/dtr/3-GitHub/emsig/emg3d-gallery/examples/download/GemPy-I.h5» [emg3d v1.0.0rc3.dev5+g0cd9e09 (format 1.0) on 2021-05-21T14:06:32.551618]. .. GENERATED FROM PYTHON SOURCE LINES 58-60 QC resistivity model -------------------- .. GENERATED FROM PYTHON SOURCE LINES 60-67 .. code-block:: default fgrid.plot_3d_slicer( fmodel.property_x, zslice=-1000, pcolor_opts={'norm': LogNorm(vmin=0.3, vmax=500)} ) .. image-sg:: /gallery/models/images/sphx_glr_GemPy-I_001.png :alt: GemPy I :srcset: /gallery/models/images/sphx_glr_GemPy-I_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 68-70 Compute some example CSEM data with it -------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 70-107 .. code-block:: default # Source: x-directed electric-source at (1000, 1000, -500) src_coo = [1000, 1000, -500, 0, 0] frequency = 1.0 # Hz # Computational grid grid = emg3d.construct_mesh( frequency=frequency, center=src_coo[:3], properties=[0.3, 200, 1000], domain=([0, 2000], [0, 2000], [-2000, 0]), seasurface=0, center_on_edge=False, ) grid # Get the computational model model = fmodel.interpolate_to_grid(grid) # Compute the response efield = emg3d.solve_source( model=model, source=emg3d.TxElectricDipole(src_coo), frequency=frequency, verb=1, ) # Plot the response grid.plot_3d_slicer( efield.fx.ravel('F'), view='abs', v_type='Ex', zslice=-1000, xlim=(-500, 2500), ylim=(-500, 2500), zlim=(-2000, 50), pcolor_opts={'norm': LogNorm(vmin=5e-12, vmax=5e-9)}, ) .. image-sg:: /gallery/models/images/sphx_glr_GemPy-I_002.png :alt: GemPy I :srcset: /gallery/models/images/sphx_glr_GemPy-I_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d :: 7.9e-07; 2(9); 0:00:27; CONVERGED .. GENERATED FROM PYTHON SOURCE LINES 108-149 Reproduce the model ------------------- .. note:: The following sections are about how to reproduce the model. For this you have to install ``gempy``. The code example and the ``GemPy-I.h5``-file used in the gallery were created on 2021-05-21 with ``gempy=2.2.9`` and ``pandas=1.2.4``. Get and initiate the *simple_fault_model* ----------------------------------------- **Changes made to the original model** (differences between the files `simple_fault_model_*.csv` and `simple_fault_model_*_geophy.csv`): Changed the stratigraphic unit names, and moved the model 2 km down. Instead of reading a csv-file we could also initiate an empty instance and then add points and orientations after that by, e.g., providing numpy arrays. .. code-block:: python import gempy as gempy import numpy as np # Initiate a model geo_model = gempy.create_model('GemPy-I') # Location of data files. data_url = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/' data_url += 'master/data/input_data/tut_chapter1/' # Importing the data from CSV-files and setting extent and resolution # This is a regular grid, mainly for plotting purposes gempy.init_data( geo_model, [0, 2000., 0, 2000., -2000, 40.], [50, 50, 51], path_o=data_url+"simple_fault_model_orientations_geophy.csv", path_i=data_url+"simple_fault_model_points_geophy.csv", ) .. GENERATED FROM PYTHON SOURCE LINES 152-185 Initiate the stratigraphies and faults, and add an air layer ------------------------------------------------------------ .. code-block:: python # Add an air-layer: Horizontal layer at z=0m geo_model.add_surfaces('air') geo_model.add_surface_points(0, 0, 0, 'air') geo_model.add_surface_points(0, 0, 0, 'air') geo_model.add_orientations(0, 0, 0, 'air', [0, 0, 1]) # Add Series for the air layer; this series will not be cut by the fault geo_model.add_series('Air_Series') geo_model.modify_order_series(2, 'Air_Series') gempy.map_series_to_surfaces(geo_model, {'Air_Series': 'air'}) # Map the different series gempy.map_series_to_surfaces( geo_model, { "Fault_Series": 'fault', "Air_Series": ('air'), "Strat_Series": ('seawater', 'overburden', 'target', 'underburden', 'basement') }, remove_unused_series=True ) geo_model.rename_series({'Main_Fault': 'Fault_Series'}) # Set which series the fault series is cutting geo_model.set_is_fault('Fault_Series') geo_model.faults.faults_relations_df .. GENERATED FROM PYTHON SOURCE LINES 188-207 Compute the model with GemPy ---------------------------- .. code-block:: python # Set the interpolator. gempy.set_interpolator( geo_model, compile_theano=True, theano_optimizer='fast_compile', verbose=[] ) # Compute it. sol = gempy.compute_model(geo_model, compute_mesh=True) # Plot lithologies (colour-code corresponds to lithologies) _ = gempy.plot_2d(geo_model, cell_number=25, direction='y', show_data=True) .. GENERATED FROM PYTHON SOURCE LINES 210-232 Get id's for a discretize mesh ------------------------------ We could define the resistivities before, but currently it is difficult for GemPy to interpolate for something like resistivities with a very wide range of values (several orders of magnitudes). So we can simply map it here to the ``id`` (Note: GemPy does not do interpolation for cells which lie in different stratigraphies, so the id is always in integer). .. code-block:: python # First we create a detailed discretize-mesh to store the resistivity # model and use it in other examples as well. hxy = np.ones(100)*100 hz = np.ones(100)*25 fgrid = emg3d.TensorMesh([hxy, hxy, hz], origin=(-4000, -4000, -2400)) # Get the solution at cell centers of our grid. sol = gempy.compute_model(geo_model, at=fgrid.gridCC) # Show the surfaces. geo_model.surfaces .. GENERATED FROM PYTHON SOURCE LINES 235-256 Replace id's by resistivities ----------------------------- .. code-block:: python # Now, we convert the id's to resistivities res = sol.custom[0][0, :fgrid.n_cells] res[res == 1] = 1e8 # air # id=2 is the fault res[np.round(res) == 3] = 0.3 # sea water res[np.round(res) == 4] = 1.0 # overburden res[np.round(res) == 5] = 50 # resistive layer res[np.round(res) == 6] = 1.5 # underburden res[np.round(res) == 7] = 200 # resistive basement # Create an emg3d-model. fmodel = emg3d.Model(fgrid, property_x=res, mapping='Resistivity') # Store model. emg3d.save('GemPy-I.h5', model=fmodel) .. GENERATED FROM PYTHON SOURCE LINES 259-301 PyVista plot ------------ .. note:: The final cell is about how to plot the model in 3D using PyVista, for which you have to install ``pyvista``. The code example was created on 2021-05-21 with ``pyvista=0.30.1``. .. code-block:: python import pyvista import numpy as np dataset = fgrid.toVTK({'res': np.log10(fmodel.property_x.ravel('F'))}) # Create the rendering scene and add a grid axes p = pyvista.Plotter(notebook=True) p.show_grid(location='outer') # Add spatially referenced data to the scene dparams = {'rng': np.log10([0.3, 500]), 'show_edges': False} xyz = (1500, 500, -1500) p.add_mesh(dataset.slice('x', xyz), name='x-slice', **dparams) p.add_mesh(dataset.slice('y', xyz), name='y-slice', **dparams) # Add a layer as 3D p.add_mesh(dataset.threshold([1.69, 1.7]), name='vol', **dparams) # Show the scene! p.camera_position = [ (-10000, 25000, 4000), (1000, 1000, -1000), (0, 0, 1) ] p.show() .. figure:: ../../_static/images/GemPy-I.png :scale: 66 % :align: center :alt: GemPy-I model with PyVista :name: gempy-i .. GENERATED FROM PYTHON SOURCE LINES 304-306 .. code-block:: default emg3d.Report() .. raw:: html
Wed Aug 31 21:49:32 2022 CEST
OS Linux CPU(s) 4 Machine x86_64
Architecture 64bit RAM 15.5 GiB Environment Python
File system ext4
Python 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:22:55) [GCC 10.3.0]
numpy 1.22.4 scipy 1.9.0 numba 0.55.2
emg3d 1.8.0 empymod 2.2.0 xarray 2022.6.0
discretize 0.8.2 h5py 3.7.0 matplotlib 3.4.3
tqdm 4.64.0 IPython 8.4.0
Intel(R) oneAPI Math Kernel Library Version 2022.0-Product Build 20211112 for Intel(R) 64 architecture applications


.. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 29.164 seconds) **Estimated memory usage:** 124 MB .. _sphx_glr_download_gallery_models_GemPy-I.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: GemPy-I.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: GemPy-I.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_