.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tutorials/layered.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_tutorials_layered.py: 10. Layered =========== Since ``v1.8.0`` there exists the possibility to use locally approximated, layered models together with ``empymod`` to compute responses and the gradient instead of the 3D model with ``emg3d``. Here we look at this possibility by reproducing the example :ref:`sphx_glr_gallery_tutorials_gradient.py`. (Please note, however, that the layered option does not only work for the gradient, but also for pure forward modelling.) For this example we use the resistivity model created in the example :ref:`sphx_glr_gallery_models_gempy-ii.py`, together with the survey created in :ref:`sphx_glr_gallery_tutorials_simulation.py`. Limitations ----------- The limitations of the layered modelling are: - Only point and dipole sources. - Only isotropic and VTI models. Also, using the layered modelling has a few other implications, e.g.: - There are no {e;h}fields, only the fields at receiver locations. - ``gridding``, most of ``gridding_opts``, ``solver_opts``, ``receiver_interpolation``, and ``file_dir`` have no effect. - The attribute :attr:`emg3d.simulations.Simulation.jvec` is not implemented. Some considerations about speed ------------------------------- If you are concerned about speed take into account that using layered computations changes the things one has to consider. - "Regular" 3D: - Linear with number of source-frequency pairs; parallelized over it. - Depends on gridding options; linear in number of cells. - Extra receivers at almost no extra cost. - "Layered" 1D: - Linear with number of source-receiver pairs; parallelized over sources; serial over receivers. - Only for the gradient: depends linearly on the number of layers in the provided model. - Extra frequencies at little extra cost. .. GENERATED FROM PYTHON SOURCE LINES 56-67 .. code-block:: default import os import pooch import emg3d import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LogNorm, SymLogNorm # Adjust this path to a folder of your choice. data_path = os.path.join('..', 'download', '') .. GENERATED FROM PYTHON SOURCE LINES 69-76 Load survey, data, and model ---------------------------- First we load the model, survey, and accompanying data as obtained in the examples :ref:`sphx_glr_gallery_models_gempy-ii.py` and :ref:`sphx_glr_gallery_tutorials_simulation.py` (check these examples for more information). .. GENERATED FROM PYTHON SOURCE LINES 76-99 .. code-block:: default isurvey = ( 'GemPy-II-survey-A.h5', '5f2ed0b959a4f80f5378a071e6f729c6b7446898be7689ddc9bbd100a8f5bce7', ) imodel = ( 'GemPy-II.h5', 'ea8c23be80522d3ca8f36742c93758370df89188816f50cb4e1b2a6a3012d659', ) # Download model and survey. for data in [isurvey, imodel]: pooch.retrieve( 'https://raw.github.com/emsig/data/2021-05-21/emg3d/surveys/'+data[0], data[1], fname=data[0], path=data_path, ) # Load them. survey = emg3d.load(os.path.join(data_path, isurvey[0]))['survey'] true_model = emg3d.load(os.path.join(data_path, imodel[0]))['model'] grid = true_model.grid .. rst-class:: sphx-glr-script-out .. code-block:: none Data loaded from «/home/dtr/3-GitHub/emsig/emg3d-gallery/examples/download/GemPy-II-survey-A.h5» [emg3d v0.17.1.dev22+ge23c468 (format 1.0) on 2021-04-14T22:24:03.229285]. Data loaded from «/home/dtr/3-GitHub/emsig/emg3d-gallery/examples/download/GemPy-II.h5» [emg3d v1.0.0rc3.dev5+g0cd9e09 (format 1.0) on 2021-05-21T18:40:16.721968]. .. GENERATED FROM PYTHON SOURCE LINES 100-106 Create an initial model ----------------------- To create an initial model we take the true model, but set all subsurface resistivities to 1 Ohm.m. So we are left with a homogeneous three-layer model air-seawater-subsurface, which includes the topography of the seafloor. .. GENERATED FROM PYTHON SOURCE LINES 106-135 .. code-block:: default # Overwrite all subsurface resistivity values with 1.0 model = true_model.copy() res = model.property_x subsurface = (res > 0.5) & (res < 1000) res[subsurface] = 1.0 model.property_x = res # QC the initial model and the survey. popts = {'norm': LogNorm(vmin=0.3, vmax=200)} grid.plot_3d_slicer( model.property_x, xslice=12000, yslice=7500, pcolor_opts=popts ) # Plot survey in figure above fig = plt.gcf() fig.suptitle('Initial resistivity model (Ohm.m)') axs = fig.get_children() rec_coords = survey.receiver_coordinates() src_coords = survey.source_coordinates() axs[1].plot(rec_coords[0], rec_coords[1], 'bv') axs[2].plot(rec_coords[0], rec_coords[2], 'bv') axs[3].plot(rec_coords[2], rec_coords[1], 'bv') axs[1].plot(src_coords[0], src_coords[1], 'r*') axs[2].plot(src_coords[0], src_coords[2], 'r*') axs[3].plot(src_coords[2], src_coords[1], 'r*') plt.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_layered_001.png :alt: Initial resistivity model (Ohm.m) :srcset: /gallery/tutorials/images/sphx_glr_layered_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 136-155 Options for layered computation ------------------------------- If ``layered=True``, one can provide a ``layered_opts``-dict to a simulation. For more information about this see towards the end of this example. The simulation uses some defaults values, if they are not provided. By default a local 1D model is generated by creating a right elliptic cylinder around source and receiver position, and average the values for each layer within this cylinder. You can read more about the local extraction and the ellipse options in the functions :attr:`emg3d.models.Model.extract_1d` and :func:`emg3d.maps.ellipse_indices`. Note that using layered computations means that we do not have to provide any gridding options, as there is no gridding happening. Create the Simulation --------------------- .. GENERATED FROM PYTHON SOURCE LINES 155-169 .. code-block:: default simulation = emg3d.simulations.Simulation( name="Initial Model", survey=survey, model=model, max_workers=4, layered=True, # Set layered to True! # layered_opts=..., # Parameters for the local 1D extractions. ) # Let's QC our Simulation instance simulation .. raw:: html

Simulation «Initial Model»



.. GENERATED FROM PYTHON SOURCE LINES 170-177 The simulation indicates that it uses layered computations with the method ``'cylinder'`` and certain ellipse-parameters. Compute Gradient ---------------- Let's compute the gradient (with ``layered=True``!). .. GENERATED FROM PYTHON SOURCE LINES 177-181 .. code-block:: default grad = simulation.gradient .. rst-class:: sphx-glr-script-out .. code-block:: none Compute empymod 0/3 [00:00] Compute empymod ###3 1/3 [00:01] Compute empymod ########## 3/3 [00:01] Compute empymod 0/3 [00:00] Compute empymod ###3 1/3 [00:32] Compute empymod ######6 2/3 [00:36] Compute empymod ########## 3/3 [00:36] .. GENERATED FROM PYTHON SOURCE LINES 182-184 QC Gradient ''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 184-212 .. code-block:: default # Set the gradient of air and water to NaN. # This will eventually move directly into emgd3 (active and inactive cells). grad[~subsurface] = np.nan # Plot the gradient grid.plot_3d_slicer( grad.ravel('F'), xslice=12000, yslice=7500, zslice=-4000, pcolor_opts={ 'cmap': 'RdBu_r', 'norm': SymLogNorm(linthresh=1e-2, base=10, vmin=-1e1, vmax=1e1) } ) # Add survey fig = plt.gcf() fig.suptitle('Gradient of the misfit function') axs = fig.get_children() axs[1].plot(rec_coords[0], rec_coords[1], 'bv') axs[2].plot(rec_coords[0], rec_coords[2], 'bv') axs[3].plot(rec_coords[2], rec_coords[1], 'bv') axs[1].plot(src_coords[0], src_coords[1], 'r*') axs[2].plot(src_coords[0], src_coords[2], 'r*') axs[3].plot(src_coords[2], src_coords[1], 'r*') plt.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_layered_002.png :alt: Gradient of the misfit function :srcset: /gallery/tutorials/images/sphx_glr_layered_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 213-252 Even though we use layered (1D) models, the gradient looks three dimensional. This is because each source-receiver pair has its own cylinder, from which it generates a 1D model. The resulting gradient is then projected back to the same cylinder. This yields, summing up all source-receiver pairs, a three-dimensional gradient. Note that this gradient is obtained using the finite-difference method, by perturbing each layer by a tiny bit. This is in contrast to the gradient obtained with ``emg3d``, where the adjoint-state method is used. Compare this gradient to the adjoint-state gradient, obtained by true 3D computations, in example :ref:`sphx_glr_gallery_tutorials_gradient.py`. A few words about the 1D model extraction ----------------------------------------- There are five possible methods: cylinder, prism, midpoint, source, and receiver. However, in reality they reduce to really two different possibility: - ``'midpoint'``: If midpoint is chosen, the selected 1D model is simply the column midway between the two provided points, usually source and receiver. ``'source'`` and ``'receiver'`` only differ in that both points are set to source or receiver, respectively, meaning that the column at source or receiver location is selected. - ``'cylinder'``: If cylinder is chosen, a right elliptic cylinder around the two points is selected. How the ellipse is selected depends on a few other parameters, see functions :attr:`emg3d.models.Model.extract_1d` and :func:`emg3d.maps.ellipse_indices`. The only difference for ``'prism'`` is that it uses the rectangular, coordinate-system aligned envelope of the ellipse. The default values are: ``method='cylinder'``. The corresponding default values of the ellipse-options are ``factor=1.2``, ``minor=0.8``, and ``radius`` is set to one skin depth using the lowest frequency and the lowest conductivity value in the lowest layer in vertical direction. Let's look at the default value the simulation has chosen in this case: .. GENERATED FROM PYTHON SOURCE LINES 252-257 .. code-block:: default lopts = simulation.layered_opts lopts .. rst-class:: sphx-glr-script-out .. code-block:: none {'method': 'cylinder', 'ellipse': {'radius': 711.762543223444, 'factor': 1.2, 'minor': 0.8}} .. GENERATED FROM PYTHON SOURCE LINES 258-263 So we see that it chose the default method ``'cylinder'``, the default values ``factor=1.2`` and ``minor=0.8``, and the ``radius`` is set to roughly 712 m. We select now one source and one receiver, to QC how the selected cylinder looks like. .. GENERATED FROM PYTHON SOURCE LINES 263-294 .. code-block:: default src = survey.sources['TxED-1'] rec = survey.receivers['RxEP-30'] layered, imat = true_model.extract_1d( lopts['method'], src.center[:2], rec.center[:2], lopts['ellipse'], return_imat=True ) values = true_model.property_x.copy() values[~np.array(imat, dtype=bool), :] = np.nan xy = (src.center[:2] + rec.center[:2])/2 grid.plot_3d_slicer( values, xslice=xy[0], yslice=xy[1], zslice=-4000, pcolor_opts=popts ) fig = plt.gcf() fig.suptitle('Right elliptic cylinder for this source and receiver location') axs = fig.get_children() axs[1].plot(rec.center[0], rec.center[1], 'bv') axs[1].plot(src.center[0], src.center[1], 'r*') axs[2].plot(rec.center[0], rec.center[2], 'bv') axs[3].plot(rec.center[2], rec.center[1], 'bv') axs[2].plot(src.center[0], src.center[2], 'r*') axs[3].plot(src.center[2], src.center[1], 'r*') plt.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_layered_003.png :alt: Right elliptic cylinder for this source and receiver location :srcset: /gallery/tutorials/images/sphx_glr_layered_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 295-297 We can compare the five different methods, and see what 1D models they extract from the 3D model. .. GENERATED FROM PYTHON SOURCE LINES 297-332 .. code-block:: default fig, axs = plt.subplots( 1, 5, figsize=(9, 4), sharex=True, sharey=True, constrained_layout=True ) axs = axs.ravel() pdepth = np.repeat(true_model.grid.nodes_z[1:-1], 2) ellipse = lopts['ellipse'] methods = ['source', 'receiver', 'midpoint', 'cylinder', 'prism'] for i, method in enumerate(methods): axs[i].set_title(method) p0 = src.center[:2] p1 = rec.center[:2] if method == 'source': p1 = p0 method = 'midpoint' elif method == 'receiver': p0 = p1 method = 'midpoint' oned = true_model.extract_1d(method, p0, p1, ellipse) res = np.repeat(oned.property_x, 2)[1:-1] axs[i].plot(res, pdepth) axs[0].set_xscale('log') axs[0].set_xlim([0.2, 60]) plt.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_layered_004.png :alt: source, receiver, midpoint, cylinder, prism :srcset: /gallery/tutorials/images/sphx_glr_layered_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 333-340 Ellipse ------- As mentioned before, the ellipse is selected using the function :func:`emg3d.maps.ellipse_indices`. The documentation of that function contains a lot of useful information about it. However, here is a function that illustrates the ellipse with all its parameters. .. GENERATED FROM PYTHON SOURCE LINES 340-422 .. code-block:: default def plot_ellipse(ax, p0, p1, radius, factor=1.0, minor=1.0, check_foci=True): """Plot annotated explanation of emg3d.maps.ellipse_ind().""" # Cast input points to numpy arrays. p0, p1 = np.array(p0), np.array(p1) # Center coordinates cx, cy = (p0 + p1) / 2 # Adjacent and opposite sides dx, dy = (p1 - p0) / 2 # c: linear eccentricity dxy = np.linalg.norm([dx, dy]) # a: semi-major axis major = max(dxy * factor, dxy + radius) # b: semi-minor axis minor = max(minor * major, radius) if check_foci: minor = max(minor, np.sqrt(abs(major**2 - dxy**2))) # Distance to foci df = np.sqrt(major**2 - minor**2) # Take care of division by zero if dy == 0.0: cos, sin = 1.0, 0.0 else: cos, sin = dx/dxy, dy/dxy # Parametrized ellipse t = np.linspace(0, 2*np.pi, 101) cost, sint = np.cos(t), np.sin(t) x = cx + major*cos*cost - minor*sin*sint y = cy + major*sin*cost + minor*cos*sint # Verteces vr = [cx + major*cos, cy + major*sin] vl = [cx - major*cos, cy - major*sin] # Upper co-vertex vu = [cx - minor*sin, cy + minor*cos] # Foci fr = [cx + df*cos, cy + df*sin] fl = [cx - df*cos, cy - df*sin] # PLOTTING # Title with radius and _actual_ minor factor. ax.set_title(f"r={radius:.1f}; " f"f={factor:.2f}; " f"m={1.0 if major == 0.0 else minor/major:.2f}", fontsize=10) # Draw a, b, c, and r ax.plot([cx, vr[0]], [cy, vr[1]], 'C4', label=r'a = c + r') ax.plot([cx, vu[0]], [cy, vu[1]], 'C3', label=r'b = minor · a') ax.plot([p0[0], cx], [p0[1], cy], 'C0', label=r'c = off / 2') ax.plot([vl[0], p0[0]], [vl[1], p0[1]], 'C5', label=r'r = radius') # Draw the two points and their circles with radius for p in [p0, p1]: ax.plot(*p, 'o', c='k', ms=8) ax.plot(p[0] + radius*cost, p[1] + radius*sint, '.7') # Draw the foci ax.plot([fr[0], fl[0]], [fr[1], fl[1]], 'C2x', label='Foci') # Draw the envelopping square ax.plot([x.min(), x.max(), x.max(), x.min(), x.min()], [y.min(), y.min(), y.max(), y.max(), y.min()], '.7') # Draw the ellipse ax.plot(x, y, 'C1') # Ensure square axes ax.axis('equal') .. GENERATED FROM PYTHON SOURCE LINES 423-425 Here the annotation for the above selection. However, you can use this function to test the influence of the different parameters. .. GENERATED FROM PYTHON SOURCE LINES 425-438 .. code-block:: default fig, ax = plt.subplots(constrained_layout=True) plot_ellipse( ax, p0=src.center[:2], p1=rec.center[:2], radius=lopts['ellipse']['radius'], factor=lopts['ellipse']['factor'], minor=lopts['ellipse']['minor'], ) ax.legend(bbox_to_anchor=(1.05, 0.7)) plt.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_layered_005.png :alt: r=711.8; f=1.20; m=0.80 :srcset: /gallery/tutorials/images/sphx_glr_layered_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 439-441 .. code-block:: default emg3d.Report() .. raw:: html
Wed Aug 31 22:03:18 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 43.191 seconds) **Estimated memory usage:** 165 MB .. _sphx_glr_download_gallery_tutorials_layered.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: layered.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: layered.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_