.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/tutorials/12-seismic/plot_fwd_1_tomography_2D.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_content_tutorials_12-seismic_plot_fwd_1_tomography_2D.py: Forward Simulation for Straight Ray Tomography in 2D ==================================================== Here we module *SimPEG.seismic.straight_ray_tomography* to predict arrival time data for a synthetic velocity/slowness model. In this tutorial, we focus on the following: - How to define the survey - How to define the forward simulation - How to predict arrival time data .. GENERATED FROM PYTHON SOURCE LINES 16-19 Import Modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 19-33 .. code-block:: Python import os import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from discretize import TensorMesh from SimPEG import maps from SimPEG.seismic import straight_ray_tomography as tomo from SimPEG.utils import model_builder save_file = False .. GENERATED FROM PYTHON SOURCE LINES 34-42 Defining the Survey ------------------- Here, we define survey that will be used for the forward simulation. The survey consists of a horizontal line of point receivers at Y = 100 m and a horizontal line of point sources at Y = -100 m. The shot by each source is measured by all receivers. .. GENERATED FROM PYTHON SOURCE LINES 42-65 .. code-block:: Python # Define the locations for the sources and receivers. x = np.linspace(-100, 100, 11) y_receivers = 100 * np.ones(len(x)) y_sources = -100 * np.ones(len(x)) receiver_locations = np.c_[x, y_receivers] source_locations = np.c_[x, y_sources] # Define the list of receivers used by each source receiver_list = [tomo.Rx(receiver_locations)] # Define an empty list to store sources objects. Define each source and # provide its corresponding receivers list source_list = [] for ii in range(0, len(y_sources)): source_list.append( tomo.Src(location=source_locations[ii, :], receiver_list=receiver_list) ) # Define they tomography survey survey = tomo.Survey(source_list) .. GENERATED FROM PYTHON SOURCE LINES 66-72 Defining a Tensor Mesh ---------------------- Here, we create the tensor mesh that will be used to predict arrival time data. .. GENERATED FROM PYTHON SOURCE LINES 72-80 .. code-block:: Python dh = 10.0 # cell width N = 21 # number of cells in X and Y direction hx = [(dh, N)] hy = [(dh, N)] mesh = TensorMesh([hx, hy], "CC") .. GENERATED FROM PYTHON SOURCE LINES 81-89 Model and Mapping on Tensor Mesh -------------------------------- Here, we create the velocity model that will be used to predict the data. Since the physical parameter for straight ray tomography is slowness, we must define a mapping which converts velocity values to slowness values. The model consists of a lower velocity block within a higher velocity background. .. GENERATED FROM PYTHON SOURCE LINES 89-122 .. code-block:: Python # Define velocity of each unit in m/s background_velocity = 3000.0 block_velocity = 1500.0 # Define the model. Models in SimPEG are vector arrays. model = background_velocity * np.ones(mesh.nC) ind_block = model_builder.get_indices_block(np.r_[-50, 20], np.r_[50, -20], mesh.gridCC) model[ind_block] = block_velocity # Define a mapping from the model (velocity) to the slowness. If your model # consists of slowness values, you can use *maps.IdentityMap*. model_mapping = maps.ReciprocalMap() # Plot Velocity Model fig = plt.figure(figsize=(6, 5.5)) ax1 = fig.add_axes([0.15, 0.15, 0.65, 0.75]) mesh.plot_image(model, ax=ax1, grid=True, pcolor_opts={"cmap": "viridis"}) ax1.set_xlabel("x (m)") ax1.set_ylabel("y (m)") ax1.plot(x, y_sources, "ro") # source locations ax1.plot(x, y_receivers, "ko") # receiver locations ax2 = fig.add_axes([0.82, 0.15, 0.05, 0.75]) norm = mpl.colors.Normalize(vmin=np.min(model), vmax=np.max(model)) cbar = mpl.colorbar.ColorbarBase( ax2, norm=norm, orientation="vertical", cmap=mpl.cm.viridis ) cbar.set_label("$Velocity (m/s)$", rotation=270, labelpad=15, size=12) .. image-sg:: /content/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_001.png :alt: plot fwd 1 tomography 2D :srcset: /content/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 123-129 Simulation: Arrival Time ------------------------ Here we demonstrate how to predict arrival time data for the 2D straight ray tomography problem using the 2D Integral formulation. .. GENERATED FROM PYTHON SOURCE LINES 129-137 .. code-block:: Python # Define the forward simulation. To do this we need the mesh, the survey and # the mapping from the model to the slowness values on the mesh. simulation = tomo.Simulation(mesh, survey=survey, slownessMap=model_mapping) # Compute predicted data for some model dpred = simulation.dpred(model) .. GENERATED FROM PYTHON SOURCE LINES 138-141 Plotting ----------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 141-161 .. code-block:: Python n_source = len(source_list) n_receiver = len(x) dpred_plotting = dpred.reshape(n_receiver, n_source) fig = plt.figure(figsize=(8, 5)) ax = fig.add_subplot(111) obs_string = [] for ii in range(0, n_source): ax.plot(x, dpred_plotting[:, ii]) obs_string.append("source {}".format(ii + 1)) ax.set_xlim(np.min(x), np.max(x)) ax.set_xlabel("x (m)") ax.set_ylabel("arrival time (s)") ax.set_title("Positions vs. Arrival Time") ax.legend(obs_string, loc="upper right") .. image-sg:: /content/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_002.png :alt: Positions vs. Arrival Time :srcset: /content/tutorials/12-seismic/images/sphx_glr_plot_fwd_1_tomography_2D_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 162-167 Optional: Exporting Results --------------------------- Write the data and true model .. GENERATED FROM PYTHON SOURCE LINES 167-189 .. code-block:: Python if save_file: dir_path = os.path.dirname(tomo.__file__).split(os.path.sep)[:-3] dir_path.extend(["tutorials", "seismic", "assets"]) dir_path = os.path.sep.join(dir_path) + os.path.sep noise = 0.05 * dpred * np.random.randn(len(dpred)) data_array = np.c_[ np.kron(x, np.ones(n_receiver)), np.kron(y_sources, np.ones(n_receiver)), np.kron(np.ones(n_source), x), np.kron(np.ones(n_source), y_receivers), dpred + noise, ] fname = dir_path + "tomography2D_data.obs" np.savetxt(fname, data_array, fmt="%.4e") output_model = model fname = dir_path + "true_model_2D.txt" np.savetxt(fname, output_model, fmt="%.4e") .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.669 seconds) **Estimated memory usage:** 8 MB .. _sphx_glr_download_content_tutorials_12-seismic_plot_fwd_1_tomography_2D.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_tomography_2D.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_fwd_1_tomography_2D.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_