# 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

Out:

Downloading https://storage.googleapis.com/simpeg/Chile_GRAV_4_Miller/Chile_GRAV_4_Miller.tar.gz
saved to: /Users/lindseyjh/git/simpeg/simpeg/examples/04-grav/Chile_GRAV_4_Miller.tar.gz
Begin calculation of distance weighting for R= 3.0
Done 0.0 %
Done 10.0 %
Done 20.0 %
Done 30.0 %
Done 40.0 %
Done 50.0 %
Done 60.0 %
Done 70.0 %
Done 80.0 %
Done 90.0 %
Done 100% ...distance weighting completed!!

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
Begin linear forward calculation: z
Done 0.0 %
Done 10.0 %
Done 20.0 %
Done 30.0 %
Done 40.0 %
Done 50.0 %
Done 60.0 %
Done 70.0 %
Done 80.0 %
Done 90.0 %
Linear forward calculation ended in: 14.069777250289917 sec

SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***
Approximated diag(JtJ) with linear operator
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  2.32e-03  1.30e+06  1.87e+01  1.30e+06    1.95e+02      0
1  1.16e-03  3.71e+03  7.00e+06  1.18e+04    1.04e+02      0
2  5.81e-04  1.85e+03  8.10e+06  6.56e+03    8.71e+01      0   Skip BFGS
3  2.91e-04  8.71e+02  9.27e+06  3.56e+03    7.20e+01      0   Skip BFGS
4  1.45e-04  3.74e+02  1.04e+07  1.89e+03    5.69e+01      0   Skip BFGS
5  7.26e-05  1.48e+02  1.15e+07  9.84e+02    5.21e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.5000837959458159 eps_q: 0.5000837959458159
delta phim: 1.314e+06
6  3.63e-05  5.83e+01  2.21e+07  8.61e+02    3.72e+01      0   Skip BFGS
delta phim: 3.064e+00
7  8.53e-05  3.54e+01  2.66e+07  2.31e+03    6.81e+01      0   Skip BFGS
delta phim: 3.104e-01
8  5.52e-05  1.75e+02  2.88e+07  1.77e+03    5.10e+01      0
delta phim: 2.697e-01
9  5.52e-05  9.25e+01  3.34e+07  1.94e+03    1.39e+02      0
delta phim: 2.665e-01
10  3.43e-05  1.94e+02  3.42e+07  1.37e+03    2.18e+02      0
delta phim: 2.283e-01
11  3.43e-05  7.61e+01  3.88e+07  1.41e+03    1.03e+02      0
delta phim: 2.376e-01
12  3.43e-05  1.37e+02  3.98e+07  1.50e+03    1.94e+02      0
delta phim: 1.844e-01
13  3.43e-05  1.03e+02  4.22e+07  1.55e+03    1.79e+02      0
delta phim: 1.698e-01
14  2.38e-05  1.49e+02  4.22e+07  1.16e+03    8.75e+01      0
delta phim: 1.448e-01
15  2.38e-05  7.86e+01  4.57e+07  1.17e+03    7.36e+01      0
delta phim: 1.433e-01
16  2.38e-05  1.26e+02  4.49e+07  1.20e+03    6.94e+01      0
delta phim: 9.322e-02
17  2.38e-05  1.24e+02  4.56e+07  1.21e+03    1.71e+02      3
delta phim: 1.200e-01
18  2.38e-05  1.30e+02  4.45e+07  1.19e+03    8.30e+01      0
delta phim: 2.610e-02
19  1.63e-05  1.54e+02  4.38e+07  8.70e+02    1.16e+02      1   Skip BFGS
delta phim: 6.995e-02
20  1.63e-05  7.77e+01  4.58e+07  8.26e+02    1.63e+02      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 4.3789e+01 <= tolF*(1+|f0|) = 1.3018e+05
0 : |xc-x_last| = 1.4757e+00 <= tolX*(1+|x0|) = 1.0400e-01
0 : |proj(x-g)-x|    = 1.6269e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.6269e+02 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      20    <= iter          =     20
------------------------- DONE! -------------------------
/Users/lindseyjh/git/simpeg/simpeg/examples/04-grav/plot_laguna_del_maule_inversion.py:260: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
plt.show()


import os
import shutil
import tarfile
import SimPEG.PF as PF
from SimPEG import (
Maps, Regularization, Optimization, DataMisfit,
InvProblem, Directives, Inversion, Utils
)
import matplotlib.pyplot as plt
import numpy as np

def run(plotIt=True, cleanAfterRun=True):

# unzip the tarfile
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.
# %%
# Read in the input file which included all parameters at once
# (mesh, topo, model, survey, inv param, etc.)
driver = PF.GravityDriver.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

# define gravity survey locations
rxLoc = survey.srcField.rxList[0].locs

# define gravity data and errors
d = survey.dobs
wd = survey.std

# Get the active cells
active = driver.activeCells
nC = len(active)  # 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.nCx/2)
# %%
# Now that we have a model and a survey we can build the linear system ...
# Create the forward model operator
prob = PF.Gravity.GravityIntegral(mesh, rhoMap=staticCells,
actInd=active)
prob.solverOpts['accuracyTol'] = 1e-4

# Pair the survey and problem
survey.pair(prob)

# Apply depth weighting
wr = PF.Magnetics.get_dist_wgt(mesh, rxLoc, active, wgtexp,
np.min(mesh.hx)/4.)
wr = wr**2.

# %% Create inversion objects
reg = Regularization.Sparse(mesh, indActive=active,
reg.mref = driver.mref[dynamic]
reg.cell_weights = wr * mesh.vol[active]
reg.norms = np.c_[0., 1., 1., 1.]
# 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-3)

# Define misfit function (obs-calc)
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.W = 1./wd

# create the default L2 inverse problem from the above objects
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)

# Specify how the initial beta is found
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e-2)

# 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.Update_IRLS(
f_min_change=1e-4, maxIRLSiter=40, beta_tol=5e-1,
betaSearch=False)

# Preconditioning refreshing for each IRLS iteration
update_Jacobi = Directives.UpdatePreconditioner()

# Create combined the L2 and Lp problem
inv = Inversion.BaseInversion(invProb,
directiveList=[IRLS, update_Jacobi, betaest])

# %%
# Run L2 and Lp inversion
mrec = inv.run(mstart)

if cleanAfterRun:
shutil.rmtree(basePath)

# %%
if plotIt:
# Plot observed data
Utils.PlotUtils.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.plotSlice(L2out, ax=ax, normal='Z', ind=-16,
clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
c='gray', linestyle='--')
plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
plt.title('Z: ' + str(mesh.vectorCCz[-16]) + ' m')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)')
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.plotSlice(L2out, ax=ax, normal='Z', ind=-27,
clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
c='gray', linestyle='--')
plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
plt.title('Z: ' + str(mesh.vectorCCz[-27]) + ' m')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)')
cb = plt.colorbar(dat1[0], orientation="vertical",
ticks=np.linspace(vmin, vmax, 4))
cb.set_label('Density (g/cc$^3$)')

ax = plt.subplot(212)
mesh.plotSlice(L2out, ax=ax, normal='Y', ind=yslice,
clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
plt.title('Cross Section')
plt.xlabel('Easting(m)')
plt.ylabel('Elevation')
cb = plt.colorbar(dat1[0], 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) +
': $\epsilon_p$ = ' + str(round(reg.eps_p, 1)) +
': $\epsilon_q$ = ' + str(round(reg.eps_q, 2)))
ax = plt.subplot(221)
dat = mesh.plotSlice(Lpout, ax=ax, normal='Z', ind=-16,
clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
c='gray', linestyle='--')
plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
plt.title('Z: ' + str(mesh.vectorCCz[-16]) + ' m')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)')
cb = plt.colorbar(dat[0], orientation="vertical",
ticks=np.linspace(vmin, vmax, 4))
cb.set_label('Density (g/cc$^3$)')

ax = plt.subplot(222)
dat = mesh.plotSlice(Lpout, ax=ax, normal='Z', ind=-27,
clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
plt.plot(np.array([mesh.vectorCCx[0], mesh.vectorCCx[-1]]),
np.array([mesh.vectorCCy[yslice], mesh.vectorCCy[yslice]]),
c='gray', linestyle='--')
plt.scatter(rxLoc[0:, 0], rxLoc[0:, 1], color='k', s=1)
plt.title('Z: ' + str(mesh.vectorCCz[-27]) + ' m')
plt.xlabel('Easting (m)')
plt.ylabel('Northing (m)')
cb = plt.colorbar(dat[0], orientation="vertical",
ticks=np.linspace(vmin, vmax, 4))
cb.set_label('Density (g/cc$^3$)')

ax = plt.subplot(212)
dat = mesh.plotSlice(Lpout, ax=ax, normal='Y', ind=yslice,
clim=(vmin, vmax), pcolorOpts={'cmap': 'bwr'})
plt.title('Cross Section')
plt.xlabel('Easting (m)')
plt.ylabel('Elevation (m)')
cb.set_label('Density (g/cc$^3$)')