.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/tutorials/05-dcr/plot_fwd_1_dcr_sounding.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_content_tutorials_05-dcr_plot_fwd_1_dcr_sounding.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_content_tutorials_05-dcr_plot_fwd_1_dcr_sounding.py:


Simulate a 1D Sounding over a Layered Earth
===========================================

Here we use the module *SimPEG.electromangetics.static.resistivity* to predict
sounding data over a 1D layered Earth. In this tutorial, we focus on the following:

    - General definition of sources and receivers
    - How to define the survey
    - How to predict voltage or apparent resistivity data
    - The units of the model and resulting data

For this tutorial, we will simulate sounding data over a layered Earth using
a Wenner array. The end product is a sounding curve which tells us how the
electrical resistivity changes with depth.

.. GENERATED FROM PYTHON SOURCE LINES 22-25

Import modules
--------------


.. GENERATED FROM PYTHON SOURCE LINES 25-42

.. code-block:: Python


    import os
    import numpy as np
    import matplotlib as mpl
    import matplotlib.pyplot as plt

    from SimPEG import maps
    from SimPEG.electromagnetics.static import resistivity as dc
    from SimPEG.utils import plot_1d_layer_model

    mpl.rcParams.update({"font.size": 16})

    write_output = False

    # sphinx_gallery_thumbnail_number = 2









.. GENERATED FROM PYTHON SOURCE LINES 43-51

Create Survey
-------------

Here we demonstrate a general way to define sources and receivers.
For pole and dipole sources, we must define the A or AB electrode locations,
respectively. For the pole and dipole receivers, we must define the M or
MN electrode locations, respectively.


.. GENERATED FROM PYTHON SOURCE LINES 51-86

.. code-block:: Python


    a_min = 20.0
    a_max = 500.0
    n_stations = 25

    # Define the 'a' spacing for Wenner array measurements for each reading
    electrode_separations = np.linspace(a_min, a_max, n_stations)

    source_list = []  # create empty array for sources to live

    for ii in range(0, len(electrode_separations)):
        # Extract separation parameter for sources and receivers
        a = electrode_separations[ii]

        # AB electrode locations for source. Each is a (1, 3) numpy array
        A_location = np.r_[-1.5 * a, 0.0, 0.0]
        B_location = np.r_[1.5 * a, 0.0, 0.0]

        # MN electrode locations for receivers. Each is an (N, 3) numpy array
        M_location = np.r_[-0.5 * a, 0.0, 0.0]
        N_location = np.r_[0.5 * a, 0.0, 0.0]

        # Create receivers list. Define as pole or dipole.
        receiver_list = dc.receivers.Dipole(
            M_location, N_location, data_type="apparent_resistivity"
        )
        receiver_list = [receiver_list]

        # Define the source properties and associated receivers
        source_list.append(dc.sources.Dipole(receiver_list, A_location, B_location))

    # Define survey
    survey = dc.Survey(source_list)









.. GENERATED FROM PYTHON SOURCE LINES 87-96

Defining a 1D Layered Earth Model
---------------------------------

Here, we define the layer thicknesses and electrical resistivities for our
1D simulation. If we have N layers, we define N electrical resistivity
values and N-1 layer thicknesses. The lowest layer is assumed to extend to
infinity. In the case of a halfspace, the layer thicknesses would be
an empty array.


.. GENERATED FROM PYTHON SOURCE LINES 96-106

.. code-block:: Python


    # Define layer thicknesses.
    layer_thicknesses = np.r_[100.0, 100.0]

    # Define layer resistivities.
    model = np.r_[1e3, 4e3, 2e2]

    # Define mapping from model to 1D layers.
    model_map = maps.IdentityMap(nP=len(model))








.. GENERATED FROM PYTHON SOURCE LINES 107-112

Plot Resistivity Model
----------------------

Here we plot the 1D resistivity model.


.. GENERATED FROM PYTHON SOURCE LINES 112-117

.. code-block:: Python


    # Plot the 1D model
    ax = plot_1d_layer_model(layer_thicknesses, model_map * model)
    ax.set_xlabel(r"Resistivity ($\Omega m$)")




.. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_fwd_1_dcr_sounding_001.png
   :alt: plot fwd 1 dcr sounding
   :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_fwd_1_dcr_sounding_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 49.52222222222221, 'Resistivity ($\\Omega m$)')



.. GENERATED FROM PYTHON SOURCE LINES 118-125

Define the Forward Simulation and Predict DC Resistivity Data
-------------------------------------------------------------

Here we predict DC resistivity data. If the keyword argument *rhoMap* is
defined, the simulation will expect a resistivity model. If the keyword
argument *sigmaMap* is defined, the simulation will expect a conductivity model.


.. GENERATED FROM PYTHON SOURCE LINES 125-144

.. code-block:: Python


    simulation = dc.simulation_1d.Simulation1DLayers(
        survey=survey,
        rhoMap=model_map,
        thicknesses=layer_thicknesses,
    )

    # Predict data for a given model
    dpred = simulation.dpred(model)

    # Plot apparent resistivities on sounding curve
    fig = plt.figure(figsize=(11, 5))
    ax1 = fig.add_axes([0.1, 0.1, 0.75, 0.85])
    ax1.semilogy(1.5 * electrode_separations, dpred, "b")
    ax1.set_xlabel("AB/2 (m)")
    ax1.set_ylabel(r"Apparent Resistivity ($\Omega m$)")
    plt.show()





.. image-sg:: /content/tutorials/05-dcr/images/sphx_glr_plot_fwd_1_dcr_sounding_002.png
   :alt: plot fwd 1 dcr sounding
   :srcset: /content/tutorials/05-dcr/images/sphx_glr_plot_fwd_1_dcr_sounding_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 145-150

Optional: Export Data
---------------------

Export data and true model


.. GENERATED FROM PYTHON SOURCE LINES 150-172

.. 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)

        np.random.seed(145)
        noise = 0.025 * dpred * np.random.randn(len(dpred))

        data_array = np.c_[
            survey.locations_a,
            survey.locations_b,
            survey.locations_m,
            survey.locations_n,
            dpred + noise,
        ]

        fname = dir_path + "app_res_1d_data.dobs"
        np.savetxt(fname, data_array, fmt="%.4e")








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 4.295 seconds)

**Estimated memory usage:**  8 MB


.. _sphx_glr_download_content_tutorials_05-dcr_plot_fwd_1_dcr_sounding.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_fwd_1_dcr_sounding.ipynb <plot_fwd_1_dcr_sounding.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_fwd_1_dcr_sounding.py <plot_fwd_1_dcr_sounding.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_