.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/20-published/plot_laguna_del_maule_inversion.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_examples_20-published_plot_laguna_del_maule_inversion.py: PF: Gravity: Laguna del Maule Bouguer Gravity ============================================= This notebook illustrates the SimPEG code used to invert Bouguer gravity data collected at Laguna del Maule volcanic field, Chile. Refer to Miller et al 2016 EPSL for full details. We run the inversion in two steps. Firstly creating a L2 model and then applying an Lp norm to produce a compact model. Craig Miller .. GENERATED FROM PYTHON SOURCE LINES 14-339 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_laguna_del_maule_inversion_001.png :alt: plot laguna del maule inversion :srcset: /content/examples/20-published/images/sphx_glr_plot_laguna_del_maule_inversion_001.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_laguna_del_maule_inversion_002.png :alt: Smooth Inversion: Depth weight = 3.0, Z: 1087.5 m, Z: -1125.0 m, Cross Section :srcset: /content/examples/20-published/images/sphx_glr_plot_laguna_del_maule_inversion_002.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_laguna_del_maule_inversion_003.png :alt: Compact Inversion: Depth weight = 3.0: $\epsilon_p$ = 0.0: $\epsilon_q$ = 0.0, Z: 1087.5 m, Z: -1125.0 m, Cross Section :srcset: /content/examples/20-published/images/sphx_glr_plot_laguna_del_maule_inversion_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Downloading https://storage.googleapis.com/simpeg/Chile_GRAV_4_Miller/Chile_GRAV_4_Miller.tar.gz saved to: /home/vsts/work/1/s/examples/20-published/Chile_GRAV_4_Miller.tar.gz Download completed! Running inversion with SimPEG v0.23.0 simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using the default solver Pardiso and no solver_opts.*** model has any nan: 0 =============================== Projected GNCG =============================== # beta phi_d phi_m f |proj(x-g)-x| LS Comment ----------------------------------------------------------------------------- x0 has any nan: 0 0 9.10e-04 2.60e+06 1.89e+02 2.60e+06 1.95e+02 0 1 4.55e-04 1.09e+04 1.05e+08 5.87e+04 1.41e+02 0 2 2.27e-04 4.09e+03 1.15e+08 3.03e+04 1.33e+02 0 Skip BFGS 3 1.14e-04 1.45e+03 1.23e+08 1.55e+04 1.22e+02 0 Skip BFGS 4 5.69e-05 4.86e+02 1.29e+08 7.82e+03 1.06e+02 0 Skip BFGS Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 0.3497475132140182 5 5.69e-05 1.56e+02 2.69e+08 1.55e+04 1.73e+02 0 Skip BFGS 6 2.39e-05 5.29e+02 3.65e+08 9.24e+03 9.04e+01 0 7 1.00e-05 2.20e+02 4.80e+08 5.03e+03 9.94e+01 0 8 4.20e-06 1.73e+02 4.55e+08 2.08e+03 1.77e+02 1 Skip BFGS 9 9.53e-06 2.50e+01 4.07e+08 3.90e+03 1.65e+02 0 Skip BFGS 10 2.17e-05 1.11e+02 3.08e+08 6.79e+03 1.90e+02 0 11 1.20e-05 3.07e+02 2.30e+08 3.07e+03 1.16e+02 0 12 6.66e-06 1.30e+02 1.79e+08 1.32e+03 1.12e+02 0 Skip BFGS 13 1.41e-05 3.74e+01 1.44e+08 2.07e+03 1.87e+02 0 Skip BFGS 14 2.43e-05 7.80e+01 1.24e+08 3.08e+03 1.89e+02 0 15 4.20e-05 1.71e+02 9.35e+07 4.10e+03 1.88e+02 0 16 2.34e-05 3.03e+02 7.79e+07 2.13e+03 8.96e+01 0 17 1.30e-05 1.57e+02 6.80e+07 1.04e+03 6.24e+01 0 18 2.46e-05 5.95e+01 6.58e+07 1.68e+03 1.18e+02 0 Skip BFGS 19 4.06e-05 8.79e+01 6.03e+07 2.54e+03 1.54e+02 0 20 6.71e-05 1.51e+02 5.55e+07 3.88e+03 1.69e+02 0 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 1.3416e+03 <= tolF*(1+|f0|) = 2.6036e+05 0 : |xc-x_last| = 1.0259e+00 <= tolX*(1+|x0|) = 1.0400e-01 0 : |proj(x-g)-x| = 1.6866e+02 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 1.6866e+02 <= 1e3*eps = 1.0000e-02 1 : maxIter = 20 <= iter = 20 ------------------------- DONE! ------------------------- | .. code-block:: Python import os import shutil import tarfile from simpeg.potential_fields import gravity from simpeg import ( data_misfit, maps, regularization, optimization, inverse_problem, directives, inversion, ) from simpeg.utils import download, plot2Ddata import matplotlib.pyplot as plt import numpy as np from simpeg.utils.drivers.gravity_driver import GravityDriver_Inv def run(plotIt=True, cleanAfterRun=True): # Start by downloading files from the remote repository # directory where the downloaded files are url = "https://storage.googleapis.com/simpeg/Chile_GRAV_4_Miller/Chile_GRAV_4_Miller.tar.gz" downloads = download(url, overwrite=True) basePath = downloads.split(".")[0] # unzip the tarfile tar = tarfile.open(downloads, "r") tar.extractall() tar.close() input_file = basePath + os.path.sep + "LdM_input_file.inp" # %% User input # Plotting parameters, max and min densities in g/cc vmin = -0.6 vmax = 0.6 # weight exponent for default weighting wgtexp = 3.0 # %% # Read in the input file which included all parameters at once # (mesh, topo, model, survey, inv param, etc.) driver = GravityDriver_Inv(input_file) # %% # Now we need to create the survey and model information. # Access the mesh and survey information mesh = driver.mesh # survey = driver.survey data_object = driver.data # [survey, data_object] = driver.survey # define gravity survey locations rxLoc = survey.source_field.receiver_list[0].locations # define gravity data and errors d = data_object.dobs # Get the active cells active = driver.activeCells nC = int(active.sum()) # Number of active cells # Create active map to go from reduce set to full activeMap = maps.InjectActiveCells(mesh, active, -100) # Create static map static = driver.staticCells dynamic = driver.dynamicCells staticCells = maps.InjectActiveCells(None, dynamic, driver.m0[static], nC=nC) mstart = driver.m0[dynamic] # Get index of the center midx = int(mesh.shape_cells[0] / 2) # %% # Now that we have a model and a survey we can build the linear system ... # Create the forward model operator simulation = gravity.simulation.Simulation3DIntegral( survey=survey, mesh=mesh, rhoMap=staticCells, active_cells=active ) # %% Create inversion objects reg = regularization.Sparse( mesh, active_cells=active, mapping=staticCells, gradient_type="total" ) reg.reference_model = driver.mref[dynamic] reg.norms = [0.0, 1.0, 1.0, 1.0] # reg.norms = driver.lpnorms # Specify how the optimization will proceed opt = optimization.ProjectedGNCG( maxIter=20, lower=driver.bounds[0], upper=driver.bounds[1], maxIterLS=10, maxIterCG=20, tolCG=1e-4, ) # Define misfit function (obs-calc) dmis = data_misfit.L2DataMisfit(data=data_object, simulation=simulation) # create the default L2 inverse problem from the above objects invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) # Specify how the initial beta is found betaest = directives.BetaEstimate_ByEig(beta0_ratio=0.5, random_seed=518936) # IRLS sets up the Lp inversion problem # Set the eps parameter parameter in Line 11 of the # input file based on the distribution of model (DEFAULT = 95th %ile) IRLS = directives.UpdateIRLS( f_min_change=1e-4, max_irls_iterations=40, irls_cooling_factor=1.5, misfit_tolerance=5e-1, ) # Preconditioning refreshing for each IRLS iteration update_Jacobi = directives.UpdatePreconditioner() sensitivity_weights = directives.UpdateSensitivityWeights() # Create combined the L2 and Lp problem inv = inversion.BaseInversion( invProb, directiveList=[sensitivity_weights, IRLS, update_Jacobi, betaest] ) # %% # Run L2 and Lp inversion mrec = inv.run(mstart) if cleanAfterRun: os.remove(downloads) shutil.rmtree(basePath) # %% if plotIt: # Plot observed data # The sign of the data is flipped here for the change of convention # between Cartesian coordinate system (internal SimPEG format that # expects "positive up" gravity signal) and traditional gravity data # conventions (positive down). For example a traditional negative # gravity anomaly is described as "positive up" in Cartesian coordinates # and hence the sign needs to be flipped for use in SimPEG. plot2Ddata(rxLoc, -d) # %% # Write output model and data files and print misfit stats. # reconstructing l2 model mesh with air cells and active dynamic cells L2out = activeMap * invProb.l2model # reconstructing lp model mesh with air cells and active dynamic cells Lpout = activeMap * mrec # %% # Plot out sections and histograms of the smooth l2 model. # The ind= parameter is the slice of the model from top down. yslice = midx + 1 L2out[L2out == -100] = np.nan # set "air" to nan plt.figure(figsize=(10, 7)) plt.suptitle("Smooth Inversion: Depth weight = " + str(wgtexp)) ax = plt.subplot(221) dat1 = mesh.plot_slice( L2out, ax=ax, normal="Z", ind=-16, clim=(vmin, vmax), pcolor_opts={"cmap": "bwr"}, ) plt.plot( np.array([mesh.cell_centers_x[0], mesh.cell_centers_x[-1]]), np.array([mesh.cell_centers_y[yslice], mesh.cell_centers_y[yslice]]), c="gray", linestyle="--", ) plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color="k", s=1) plt.title("Z: " + str(mesh.cell_centers_z[-16]) + " m") plt.xlabel("Easting (m)") plt.ylabel("Northing (m)") plt.gca().set_aspect("equal", adjustable="box") cb = plt.colorbar( dat1[0], orientation="vertical", ticks=np.linspace(vmin, vmax, 4) ) cb.set_label("Density (g/cc$^3$)") ax = plt.subplot(222) dat = mesh.plot_slice( L2out, ax=ax, normal="Z", ind=-27, clim=(vmin, vmax), pcolor_opts={"cmap": "bwr"}, ) plt.plot( np.array([mesh.cell_centers_x[0], mesh.cell_centers_x[-1]]), np.array([mesh.cell_centers_y[yslice], mesh.cell_centers_y[yslice]]), c="gray", linestyle="--", ) plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color="k", s=1) plt.title("Z: " + str(mesh.cell_centers_z[-27]) + " m") plt.xlabel("Easting (m)") plt.ylabel("Northing (m)") plt.gca().set_aspect("equal", adjustable="box") cb = plt.colorbar( dat1[0], ax=ax, orientation="vertical", ticks=np.linspace(vmin, vmax, 4) ) cb.set_label("Density (g/cc$^3$)") ax = plt.subplot(212) mesh.plot_slice( L2out, ax=ax, normal="Y", ind=yslice, clim=(vmin, vmax), pcolor_opts={"cmap": "bwr"}, ) plt.title("Cross Section") plt.xlabel("Easting(m)") plt.ylabel("Elevation") plt.gca().set_aspect("equal", adjustable="box") cb = plt.colorbar( dat1[0], ax=ax, orientation="vertical", ticks=np.linspace(vmin, vmax, 4), cmap="bwr", ) cb.set_label("Density (g/cc$^3$)") # %% # Make plots of Lp model yslice = midx + 1 Lpout[Lpout == -100] = np.nan # set "air" to nan plt.figure(figsize=(10, 7)) plt.suptitle( "Compact Inversion: Depth weight = " + str(wgtexp) + r": $\epsilon_p$ = " + str(round(reg.objfcts[0].irls_threshold, 1)) + r": $\epsilon_q$ = " + str(round(reg.objfcts[1].irls_threshold, 2)) ) ax = plt.subplot(221) dat = mesh.plot_slice( Lpout, ax=ax, normal="Z", ind=-16, clim=(vmin, vmax), pcolor_opts={"cmap": "bwr"}, ) plt.plot( np.array([mesh.cell_centers_x[0], mesh.cell_centers_x[-1]]), np.array([mesh.cell_centers_y[yslice], mesh.cell_centers_y[yslice]]), c="gray", linestyle="--", ) plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color="k", s=1) plt.title("Z: " + str(mesh.cell_centers_z[-16]) + " m") plt.xlabel("Easting (m)") plt.ylabel("Northing (m)") plt.gca().set_aspect("equal", adjustable="box") cb = plt.colorbar( dat[0], ax=ax, orientation="vertical", ticks=np.linspace(vmin, vmax, 4) ) cb.set_label("Density (g/cc$^3$)") ax = plt.subplot(222) dat = mesh.plot_slice( Lpout, ax=ax, normal="Z", ind=-27, clim=(vmin, vmax), pcolor_opts={"cmap": "bwr"}, ) plt.plot( np.array([mesh.cell_centers_x[0], mesh.cell_centers_x[-1]]), np.array([mesh.cell_centers_y[yslice], mesh.cell_centers_y[yslice]]), c="gray", linestyle="--", ) plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color="k", s=1) plt.title("Z: " + str(mesh.cell_centers_z[-27]) + " m") plt.xlabel("Easting (m)") plt.ylabel("Northing (m)") plt.gca().set_aspect("equal", adjustable="box") cb = plt.colorbar( dat[0], ax=ax, orientation="vertical", ticks=np.linspace(vmin, vmax, 4) ) cb.set_label("Density (g/cc$^3$)") ax = plt.subplot(212) dat = mesh.plot_slice( Lpout, ax=ax, normal="Y", ind=yslice, clim=(vmin, vmax), pcolor_opts={"cmap": "bwr"}, ) plt.title("Cross Section") plt.xlabel("Easting (m)") plt.ylabel("Elevation (m)") plt.gca().set_aspect("equal", adjustable="box") cb = plt.colorbar( dat[0], ax=ax, orientation="vertical", ticks=np.linspace(vmin, vmax, 4) ) cb.set_label("Density (g/cc$^3$)") if __name__ == "__main__": run() plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 7.929 seconds) **Estimated memory usage:** 423 MB .. _sphx_glr_download_content_examples_20-published_plot_laguna_del_maule_inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_laguna_del_maule_inversion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_laguna_del_maule_inversion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_laguna_del_maule_inversion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_