notebook.community

Edit and run



In [1]:

    
from SimPEG import *
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [17]:

    
import matplotlib
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('png')
matplotlib.rcParams['savefig.dpi'] = 100



In [18]:

    
cs, ncx, ncy, ncz, npad = 20., 30, 20, 30, 12
hx = [(cs,npad,-1.4), (cs,ncx), (cs,npad,1.4)]
hy = [(cs,npad,-1.4), (cs,ncy), (cs,npad,1.4)]
hz = [(cs,npad,-1.4), (cs,ncz), (cs,npad,1.4)]
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
print ("Padding distance x: %10.5f m") % (np.sum(mesh.hx[:npad]))
print ("Padding distance z: %10.5f m") % (np.sum(mesh.hz[:npad]))
print ("Min dx: %10.5f m") % (mesh.hx.min())
print ("Min dz: %10.5f m") % (mesh.hz.min())









    



Padding distance x: 3898.57387 m
Padding distance z: 3898.57387 m
Min dx:   20.00000 m
Min dz:   20.00000 m



In [19]:

    
print mesh









    



  ---- 3-D TensorMesh ----  
   x0: -4198.57
   y0: -4098.57
   z0: -4198.57
  nCx: 54
  nCy: 44
  nCz: 54
   hx: 1133.88, 809.91, 578.51, 413.22, 295.16, 210.83, 150.59, 107.56, 76.83, 54.88, 39.20, 28.00, 30*20.00, 28.00, 39.20, 54.88, 76.83, 107.56, 150.59, 210.83, 295.16, 413.22, 578.51, 809.91, 1133.88
   hy: 1133.88, 809.91, 578.51, 413.22, 295.16, 210.83, 150.59, 107.56, 76.83, 54.88, 39.20, 28.00, 20*20.00, 28.00, 39.20, 54.88, 76.83, 107.56, 150.59, 210.83, 295.16, 413.22, 578.51, 809.91, 1133.88
   hz: 1133.88, 809.91, 578.51, 413.22, 295.16, 210.83, 150.59, 107.56, 76.83, 54.88, 39.20, 28.00, 30*20.00, 28.00, 39.20, 54.88, 76.83, 107.56, 150.59, 210.83, 295.16, 413.22, 578.51, 809.91, 1133.88



In [20]:

    
def circfun(xc, yc, r, npoint):
    theta = np.linspace(np.pi, -np.pi, npoint)
    x = r*np.cos(theta)
    y = r*np.sin(theta)
    return x+xc, y+yc



In [21]:

    
xcirc1, ycirc1 = circfun(-150., 0., 250., 60)
xcirc2, ycirc2 = circfun(150., 0., 250., 60)



In [22]:

    
Utils.meshutils.writeUBCTensorMesh(mesh, 'mesh.msh')



In [23]:

    
sigma = Utils.meshutils.readUBCTensorModel('sigma_realistic.con', mesh)



In [24]:

    
mopt = np.load('./inv3D_FD_realistic_1/model_13.npy')
active = mesh.gridCC[:,2] < 0.
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nC)
mapping = Maps.ExpMap(mesh) * actMap



In [25]:

    
Utils.meshutils.writeUBCTensorModel(mesh, mapping*mopt, 'sigest3D_realisticFD.con')



In [26]:

    
x1 = np.arange(30)*10 - 300.
y1 = np.arange(30)*10 - 150.
xyz1 = Utils.ndgrid(x1, y1, np.r_[0.])
x2 = np.arange(30)*10 + 10.
y2 = np.arange(30)*10 - 150.
xyz2 = Utils.ndgrid(x2, y2, np.r_[0.])



In [27]:

    
fig, ax = plt.subplots(1,2, figsize=(18,6))
indz1 = 21
indz2 = indz1
print mesh.vectorCCz[indz1]
mesh.plotSlice(np.log10(sigma), ind = indz1, ax = ax[0], clim=(-3, -0.5))
mesh.plotSlice(np.log10(mapping*mopt), ind = indz2, ax = ax[1])
for i in range(2):
    ax[i].plot(xyz1[:,0], xyz1[:,1], 'r.')
    ax[i].plot(xyz2[:,0], xyz2[:,1], 'b.')
#     ax[i].plot(xcirc1, ycirc1, 'r-')
#     ax[i].plot(xcirc2, ycirc2, 'b-')
    ax[i].set_xlim(-320, 320)
    ax[i].set_ylim(-150, 150)



In [28]:

    
indy=23
fig, ax = plt.subplots(1,2, figsize=(16,4))
print mesh.vectorCCy[indy]
dat0 = mesh.plotSlice(np.log10(sigma), normal='Y', ind = indy, ax = ax[0], clim=(-3, -0.5))
dat1 = mesh.plotSlice(np.log10(mapping*mopt), normal='Y', ind = indy, ax = ax[1], clim=(-3, -0.5))
for i in range(2):
    plt.colorbar(dat0[0], ax = ax[i])
    ax[i].set_xlim(-300, 300)
    ax[i].set_ylim(-700, 0.)



In [29]:

    
import simpegEM as EM
# mapping = Maps.IdentityMap(mesh)



In [30]:

    
ntx = 2
nrx1 = xyz1.shape[0]
frequency = np.r_[1., 10., 100.]



In [31]:

    
txList = []
for itrx in range(2):
    if itrx == 0:
        for freq in frequency:
            rxr = EM.FDEM.RxFDEM(xyz1,'bzr')
            rxi = EM.FDEM.RxFDEM(xyz1,'bzi')
            tx = EM.FDEM.TxFDEM(np.array([0., -150., 0.]), 'CircularLoop', freq, [rxr, rxi])
            tx.radius = 250.
            txList.append(tx)
    elif itrx == 1:
        for freq in frequency:
            rxr = EM.FDEM.RxFDEM(xyz2,'bzr')
            rxi = EM.FDEM.RxFDEM(xyz2,'bzi')
            tx = EM.FDEM.TxFDEM(np.array([0., 150., 0.]), 'CircularLoop', freq, [rxr, rxi])
            tx.radius = 250.
            txList.append(tx)

survey = EM.FDEM.SurveyFDEM(txList)

prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping, verbose=True)



In [32]:

    
survey = EM.FDEM.SurveyFDEM(txList)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping, verbose=True)
# prb.Solver = MumpsSolver
if prb.ispaired:
    prb.unpair()
if survey.ispaired:
    survey.unpair()
prb.pair(survey)



In [33]:

    
dpred = np.load('bz_FD_realistic.npy')



In [34]:

    
noise = abs(dpred)*np.random.randn(dpred.size)*0.05 
dobs = dpred+noise



In [35]:

    
dest = np.load('inv3D_FD_realistic_1/dpred_14.npy')



In [36]:

    
Dpred = dobs.reshape((900, 2, frequency.size, 2), order='F')
Dpredr1 = Dpred[:,0,:,0]
Dpredi1 = Dpred[:,1,:,0]
Dpredr2 = Dpred[:,0,:,1]
Dpredi2 = Dpred[:,1,:,1]



In [46]:

    
Dpredamp1 = np.sqrt(Dpred[:,0,:,0]**2+Dpred[:,1,:,0]**2)
Dpredamp2 = np.sqrt(Dpred[:,0,:,1]**2+Dpred[:,1,:,1]**2)



In [47]:

    
Dest = dest.reshape((900, 2, frequency.size, 2), order='F')
Destr1 = Dest[:,0,:,0]
Desti1 = Dest[:,1,:,0]
Destr2 = Dest[:,0,:,1]
Desti2 = Dest[:,1,:,1]



In [48]:

    
Destamp1 = np.sqrt(Dest[:,0,:,0]**2+Dest[:,1,:,0]**2)
Destamp2 = np.sqrt(Dest[:,0,:,1]**2+Dest[:,1,:,1]**2)



In [52]:

    
ifreq = 1
fig, ax = plt.subplots(1,2, figsize = (12, 7))

for i in range(2):
    ax[i].set_xlabel('Easting (m)', fontsize = 16)
    ax[i].set_ylabel('Northing (m)', fontsize = 16)
    cb = plt.colorbar(dat1, ax=ax[i], orientation='horizontal')
    cb.set_label('|Magnetic flux density| (T)', fontsize = 14)
ax[0].set_title('Observed data from Tx#1', fontsize = 16)
ax[1].set_title('Predicted data from Tx#1', fontsize = 16)

vmin = Utils.mkvc(Dpredamp1[:,ifreq]).min()
vmax = Utils.mkvc(Dpredamp1[:,ifreq]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpredamp1[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
# plt.colorbar(dat1, ax=ax[0], orientation = 'horizontal')
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Destamp1[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
# plt.colorbar(dat1, ax=ax[1], orientation = 'horizontal')
fig.savefig('figures/obspredFDamp10Hztx1.png')



In [45]:

    
ifreq = 1
fig, ax = plt.subplots(1,2, figsize = (12, 7))

for i in range(2):
    ax[i].set_xlabel('Easting (m)', fontsize = 16)
    ax[i].set_ylabel('Northing (m)', fontsize = 16)
    cb = plt.colorbar(dat1, ax=ax[i], orientation='horizontal')
    cb.set_label('Magnetic flux density (T)', fontsize = 14)
ax[0].set_title('Observed data from Tx#1', fontsize = 16)
ax[1].set_title('Predicted data from Tx#1', fontsize = 16)

vmin = Utils.mkvc(Dpredr1[:,ifreq]).min()
vmax = Utils.mkvc(Dpredr1[:,ifreq]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpredr1[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
# plt.colorbar(dat1, ax=ax[0], orientation = 'horizontal')
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Destr1[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
# plt.colorbar(dat1, ax=ax[1], orientation = 'horizontal')
fig.savefig('figures/obspredFDamp10Hztx1.png')



In [39]:

    
ifreq = 1
fig, ax = plt.subplots(1,2, figsize = (12, 7))
for i in range(2):
    ax[i].set_xlabel('x (m)', fontsize = 16)
    ax[i].set_ylabel('y (m)', fontsize = 16)
vmin = Utils.mkvc(Dpredi1[:,ifreq]).min()
vmax = Utils.mkvc(Dpredi1[:,ifreq]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpredi1[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
plt.colorbar(dat1, ax=ax[0], orientation = 'horizontal')
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Desti1[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
plt.colorbar(dat1, ax=ax[1], orientation = 'horizontal')









    Out[39]:





<matplotlib.colorbar.Colorbar instance at 0x000000000BEA2A48>



In [40]:

    
ifreq = 1
fig, ax = plt.subplots(1,2, figsize = (12, 7))
for i in range(2):
    ax[i].set_xlabel('x (m)', fontsize = 16)
    ax[i].set_ylabel('y (m)', fontsize = 16)
vmin = Utils.mkvc(Dpredr2[:,ifreq]).min()
vmax = Utils.mkvc(Dpredr2[:,ifreq]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpredr2[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
plt.colorbar(dat1, ax=ax[0], orientation = 'horizontal')
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Destr2[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
plt.colorbar(dat1, ax=ax[1], orientation = 'horizontal')









    Out[40]:





<matplotlib.colorbar.Colorbar instance at 0x000000000EBFD688>



In [41]:

    
ifreq = 1
fig, ax = plt.subplots(1,2, figsize = (12, 7))
for i in range(2):
    ax[i].set_xlabel('x (m)', fontsize = 16)
    ax[i].set_ylabel('y (m)', fontsize = 16)
vmin = Utils.mkvc(Dpredi2[:,ifreq]).min()
vmax = Utils.mkvc(Dpredi2[:,ifreq]).max()
dat1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpredi2[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
plt.colorbar(dat1, ax=ax[0], orientation = 'horizontal')
dat1 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Desti2[:,ifreq].reshape((30, 30), order='F'), 30, vmin=vmin, vmax=vmax)
plt.colorbar(dat1, ax=ax[1], orientation = 'horizontal')









    Out[41]:





<matplotlib.colorbar.Colorbar instance at 0x000000000F15EE48>



In [98]:

    
ifreq = 2
fig, ax = plt.subplots(1,2, figsize = (10, 4))
ax[0].hist(Dpredr1[:,ifreq])
ax[1].hist(Dpredi1[:,ifreq])









    Out[98]:





(array([   5.,    9.,   18.,   31.,   32.,   39.,   57.,   91.,  158.,  460.]),
 array([ -4.79040616e-10,  -4.31161655e-10,  -3.83282693e-10,
         -3.35403731e-10,  -2.87524769e-10,  -2.39645808e-10,
         -1.91766846e-10,  -1.43887884e-10,  -9.60089227e-11,
         -4.81299610e-11,  -2.50999259e-13]),
 <a list of 10 Patch objects>)



In [99]:

    
# from JSAnimation import IPython_display
# from matplotlib import animation
# fig, ax = plt.subplots(1,2, figsize = (12, 5))
# for i in range(2):
#     ax[i].set_xlabel('x (m)', fontsize = 16)
#     ax[i].set_ylabel('y (m)', fontsize = 16)

# def animate(itime):
#     frame1 = ax[0].contourf(xyz1[:,0].reshape((30, 30), order='F'), xyz1[:,1].reshape((30, 30), order='F'),  Dpred[:,itime,0].reshape((30, 30), order='F'), 30)
#     frame2 = ax[1].contourf(xyz2[:,0].reshape((30, 30), order='F'), xyz2[:,1].reshape((30, 30), order='F'),  Dpred[:,itime,1].reshape((30, 30), order='F'), 30)
# #     cb1 = plt.colorbar(frame1, ax = ax[0])
# #     cb2 = plt.colorbar(frame2, ax = ax[1])    
#     return frame1, frame2
# animation.FuncAnimation(fig, animate, frames=31, interval=40, blit=True)



In [28]:

    
# np.save('bzobs_FD_realistic', dobs)
# dmis = DataMisfit.l2_DataMisfit(survey)



In [28]: