.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/user-guide/tutorials/04-magnetics/plot_2a_magnetics_induced.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
:ref:`Go to the end <sphx_glr_download_content_user-guide_tutorials_04-magnetics_plot_2a_magnetics_induced.py>`
to download the full example code.
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_content_user-guide_tutorials_04-magnetics_plot_2a_magnetics_induced.py:
Forward Simulation of Total Magnetic Intensity Data
===================================================
Here we use the module *simpeg.potential_fields.magnetics* to predict magnetic
data for a magnetic susceptibility model. We simulate the data on a tensor mesh.
For this tutorial, we focus on the following:
- How to define the survey
- How to predict magnetic data for a susceptibility model
- How to include surface topography
- The units of the physical property model and resulting data
.. GENERATED FROM PYTHON SOURCE LINES 18-21
Import Modules
--------------
.. GENERATED FROM PYTHON SOURCE LINES 21-39
.. code-block:: Python
import numpy as np
from scipy.interpolate import LinearNDInterpolator
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
from discretize import TensorMesh
from discretize.utils import mkvc, active_from_xyz
from simpeg.utils import plot2Ddata, model_builder
from simpeg import maps
from simpeg.potential_fields import magnetics
write_output = False
# sphinx_gallery_thumbnail_number = 2
.. GENERATED FROM PYTHON SOURCE LINES 40-46
Topography
----------
Surface topography is defined as an (N, 3) numpy array. We create it here but
topography could also be loaded from a file.
.. GENERATED FROM PYTHON SOURCE LINES 46-52
.. code-block:: Python
[x_topo, y_topo] = np.meshgrid(np.linspace(-200, 200, 41), np.linspace(-200, 200, 41))
z_topo = -15 * np.exp(-(x_topo**2 + y_topo**2) / 80**2)
x_topo, y_topo, z_topo = mkvc(x_topo), mkvc(y_topo), mkvc(z_topo)
xyz_topo = np.c_[x_topo, y_topo, z_topo]
.. GENERATED FROM PYTHON SOURCE LINES 53-61
Defining the Survey
-------------------
Here, we define survey that will be used for the simulation. Magnetic
surveys are simple to create. The user only needs an (N, 3) array to define
the xyz locations of the observation locations, the list of field components
which are to be modeled and the properties of the Earth's field.
.. GENERATED FROM PYTHON SOURCE LINES 61-97
.. code-block:: Python
# Define the observation locations as an (N, 3) numpy array or load them.
x = np.linspace(-80.0, 80.0, 17)
y = np.linspace(-80.0, 80.0, 17)
x, y = np.meshgrid(x, y)
x, y = mkvc(x.T), mkvc(y.T)
fun_interp = LinearNDInterpolator(np.c_[x_topo, y_topo], z_topo)
z = fun_interp(np.c_[x, y]) + 10 # Flight height 10 m above surface.
receiver_locations = np.c_[x, y, z]
# Define the component(s) of the field we want to simulate as a list of strings.
# Here we simulation total magnetic intensity data.
components = ["tmi"]
# Use the observation locations and components to define the receivers. To
# simulate data, the receivers must be defined as a list.
receiver_list = magnetics.receivers.Point(receiver_locations, components=components)
receiver_list = [receiver_list]
# Define the inducing field H0 = (intensity [nT], inclination [deg], declination [deg])
inclination = 90
declination = 0
strength = 50000
source_field = magnetics.sources.UniformBackgroundField(
receiver_list=receiver_list,
amplitude=strength,
inclination=inclination,
declination=declination,
)
# Define the survey
survey = magnetics.survey.Survey(source_field)
.. GENERATED FROM PYTHON SOURCE LINES 98-103
Defining a Tensor Mesh
----------------------
Here, we create the tensor mesh that will be used for the forward simulation.
.. GENERATED FROM PYTHON SOURCE LINES 103-111
.. code-block:: Python
dh = 5.0
hx = [(dh, 5, -1.3), (dh, 40), (dh, 5, 1.3)]
hy = [(dh, 5, -1.3), (dh, 40), (dh, 5, 1.3)]
hz = [(dh, 5, -1.3), (dh, 15)]
mesh = TensorMesh([hx, hy, hz], "CCN")
.. GENERATED FROM PYTHON SOURCE LINES 112-119
Defining a Susceptibility Model
-------------------------------
Here, we create the model that will be used to predict magnetic data
and the mapping from the model to the mesh. The model
consists of a susceptible sphere in a less susceptible host.
.. GENERATED FROM PYTHON SOURCE LINES 119-164
.. code-block:: Python
# Define susceptibility values for each unit in SI
background_susceptibility = 0.0001
sphere_susceptibility = 0.01
# Find cells that are active in the forward modeling (cells below surface)
ind_active = active_from_xyz(mesh, xyz_topo)
# Define mapping from model to active cells
nC = int(ind_active.sum())
model_map = maps.IdentityMap(nP=nC) # model is a vlue for each active cell
# Define model. Models in SimPEG are vector arrays
model = background_susceptibility * np.ones(ind_active.sum())
ind_sphere = model_builder.get_indices_sphere(
np.r_[0.0, 0.0, -45.0], 15.0, mesh.cell_centers
)
ind_sphere = ind_sphere[ind_active]
model[ind_sphere] = sphere_susceptibility
# Plot Model
fig = plt.figure(figsize=(9, 4))
plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan)
ax1 = fig.add_axes([0.1, 0.12, 0.73, 0.78])
mesh.plot_slice(
plotting_map * model,
normal="Y",
ax=ax1,
ind=int(mesh.shape_cells[1] / 2),
grid=True,
clim=(np.min(model), np.max(model)),
)
ax1.set_title("Model slice at y = 0 m")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
ax2 = fig.add_axes([0.85, 0.12, 0.05, 0.78])
norm = mpl.colors.Normalize(vmin=np.min(model), vmax=np.max(model))
cbar = mpl.colorbar.ColorbarBase(ax2, norm=norm, orientation="vertical")
cbar.set_label("Magnetic Susceptibility (SI)", rotation=270, labelpad=15, size=12)
plt.show()
.. image-sg:: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_2a_magnetics_induced_001.png
:alt: Model slice at y = 0 m
:srcset: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_2a_magnetics_induced_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 165-171
Simulation: TMI Data for a Susceptibility Model
-----------------------------------------------
Here we demonstrate how to predict magnetic data for a magnetic
susceptibility model using the integral formulation.
.. GENERATED FROM PYTHON SOURCE LINES 173-175
Define the forward simulation. By setting the 'store_sensitivities' keyword
argument to "forward_only", we simulate the data without storing the sensitivities
.. GENERATED FROM PYTHON SOURCE LINES 175-186
.. code-block:: Python
simulation = magnetics.simulation.Simulation3DIntegral(
survey=survey,
mesh=mesh,
model_type="scalar",
chiMap=model_map,
active_cells=ind_active,
store_sensitivities="forward_only",
engine="choclo",
)
.. GENERATED FROM PYTHON SOURCE LINES 187-194
.. tip::
Since SimPEG v0.22.0 we can use `Choclo
<https://www.fatiando.org/choclo>`_ as the engine for running the magnetic
simulations, which results in faster and more memory efficient runs. Just
pass ``engine="choclo"`` when constructing the simulation.
.. GENERATED FROM PYTHON SOURCE LINES 196-197
Compute predicted data for a susceptibility model
.. GENERATED FROM PYTHON SOURCE LINES 197-227
.. code-block:: Python
dpred = simulation.dpred(model)
# Plot
fig = plt.figure(figsize=(6, 5))
v_max = np.max(np.abs(dpred))
ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.85])
plot2Ddata(
receiver_list[0].locations,
dpred,
ax=ax1,
ncontour=30,
clim=(-v_max, v_max),
contourOpts={"cmap": "bwr"},
)
ax1.set_title("TMI Anomaly")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("y (m)")
ax2 = fig.add_axes([0.87, 0.1, 0.03, 0.85])
norm = mpl.colors.Normalize(vmin=-np.max(np.abs(dpred)), vmax=np.max(np.abs(dpred)))
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.bwr
)
cbar.set_label("$nT$", rotation=270, labelpad=15, size=12)
plt.show()
.. image-sg:: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_2a_magnetics_induced_002.png
:alt: TMI Anomaly
:srcset: /content/user-guide/tutorials/04-magnetics/images/sphx_glr_plot_2a_magnetics_induced_002.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 228-233
Optional: Export Data
---------------------
Write the data and topography
.. GENERATED FROM PYTHON SOURCE LINES 233-250
.. code-block:: Python
if write_output:
dir_path = os.path.dirname(__file__).split(os.path.sep)
dir_path.extend(["outputs"])
dir_path = os.path.sep.join(dir_path) + os.path.sep
if not os.path.exists(dir_path):
os.mkdir(dir_path)
fname = dir_path + "magnetics_topo.txt"
np.savetxt(fname, np.c_[xyz_topo], fmt="%.4e")
np.random.seed(211)
maximum_anomaly = np.max(np.abs(dpred))
noise = 0.02 * maximum_anomaly * np.random.randn(len(dpred))
fname = dir_path + "magnetics_data.obs"
np.savetxt(fname, np.c_[receiver_locations, dpred + noise], fmt="%.4e")
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 8.071 seconds)
**Estimated memory usage:** 301 MB
.. _sphx_glr_download_content_user-guide_tutorials_04-magnetics_plot_2a_magnetics_induced.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_2a_magnetics_induced.ipynb <plot_2a_magnetics_induced.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_2a_magnetics_induced.py <plot_2a_magnetics_induced.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: plot_2a_magnetics_induced.zip <plot_2a_magnetics_induced.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_