In [1]:
from SimPEG import Utils, Maps, Mesh
%pylab inline
In [2]:
# temp = np.loadtxt("../data/pohang/dem.dat")
# inds = np.logical_and(temp[:,0]>230000, temp[:,0]<233000) & np.logical_and(temp[:,1]>286500, temp[:,1]<289500)
# topo = temp[inds,:]
# np.save("../data/pohang/topo", topo)
In [3]:
topo = np.load("../data/pohang/topo.npy")
In [4]:
# temp = np.loadtxt("../data/pohang/SP_cor_allf.dat")
fid = open("../data/pohang/finalSPfield.dat")
lines = fid.readlines()
temp = []
temp_skip = []
for line in lines:
if len(line.split()) == 4:
temp.append(map(float, line.split()[1:]))
else:
print "skip data"
# temp_skip.append(map(float, line.split()[1:]))
print line
data = np.vstack(temp)
# "refe1" 231829.19 287921.942
# "refe2" 231785.361 287857.776
In [5]:
plt.plot(topo[:,2]-topo[:,2].min())
Out[5]:
In [6]:
zmax = (topo[:,2]-topo[:,2].min()).max()
In [7]:
actdatinds = (data[:,0] < 233000) & (data[:,0] > 230281) & (data[:,1] > 286500) & (data[:,1] < 289500) & (data[:,2] < 40.)
xmax, xmin = data[actdatinds,0].max(), data[actdatinds,0].min()
ymax, ymin = data[actdatinds,1].max(), data[actdatinds,1].min()
xmed, ymed = 0.5*(xmax-xmin) + xmin, 0.5*(ymax-ymin) + ymin
xref, yref = 231829.19, 287921.942
xref_Inv, yref_Inv = 231829.19-xmed, 287921.942-ymed
topoInv = np.c_[topo[:,0]-xmed, topo[:,1]-ymed, topo[:,2]-topo[:,2].min()]
dataInv = np.c_[data[actdatinds,0]-xmed, data[actdatinds,1]-ymed, data[actdatinds,2]]
dataInv = np.vstack((dataInv, np.r_[xref_Inv, yref_Inv, 0.]))
In [8]:
Utils.plot2Ddata(topo[:,:2], topo[:,2], ncontour=100, contourOpts={"cmap":"viridis"})
plt.plot(data[actdatinds,0], data[actdatinds,1], '.')
plt.plot(xref, yref, 'ro')
Out[8]:
In [9]:
fig = plt.figure(figsize = (5,5))
ax = plt.subplot(111)
Utils.plot2Ddata(topoInv[:,:2], topoInv[:,2], ncontour=50, ax=ax, contourOpts={"cmap":"binary_r", "alpha":0.2})
Utils.plot2Ddata(dataInv[:,:2], dataInv[:,2], ncontour=100, contourOpts={"cmap":"jet", "alpha":0.6}, dataloc=True, ax=ax)
Out[9]:
In [10]:
plt.plot(dataInv[:,2])
Out[10]:
In [11]:
dx = 50.
dz = 20.
npad = 5
hxind = [(dx, npad, -1.3), (dx, 60), (dx, npad, 1.3)]
hyind = [(dx, npad, -1.3), (dx, 60), (dx, npad, 1.3)]
hzind = [(dz, npad, -1.3), (dz, 10), (dz*0.8, 16)]
mesh = Mesh.TensorMesh([hxind, hyind, hzind], "CCN")
mesh._x0 = np.r_[mesh.x0[0], mesh.x0[1], mesh.x0[2]+zmax]
In [12]:
mesh.hz[:3].sum()
Out[12]:
In [13]:
from SimPEG.EM.Static.Utils import StaticUtils
In [14]:
print mesh
In [15]:
from scipy.sparse import block_diag
In [16]:
actind = Utils.surface2ind_topo(mesh, topoInv, gridLoc="N")
actMap = Maps.InjectActiveCells(mesh, actind, 0.)
mesh2D, topoCC = StaticUtils.gettopoCC(mesh, ~actind)
In [17]:
mesh2D.plotImage(topoCC)
Out[17]:
In [19]:
xyzlocInv = StaticUtils.drapeTopotoLoc(mesh, topoInv, dataInv[:,:2], actind=actind)
xyzlocInv = np.c_[xyzlocInv[:,:2], xyzlocInv[:,2]]
In [20]:
sigma = np.ones(mesh.nC)*1e-3
sigma[~actind] = 1e-8
In [21]:
indy = 40
indstemp = abs(xyzlocInv[:,1]-mesh.vectorCCy[indy]) < 10.
mesh.plotSlice(actind, normal="Z", ind=27, grid=True)
plt.plot(dataInv[:,0], dataInv[:,1], 'w.')
plt.plot(dataInv[indstemp,0], dataInv[indstemp,1], 'ro')
mesh.plotSlice(sigma, normal="Y", grid=True, ind = indy)
plt.plot(xyzlocInv[indstemp,0], xyzlocInv[indstemp,2], 'w.')
Out[21]:
In [22]:
mesh.plotSlice(actind, normal="Z", ind=20, grid=True)
plt.plot(dataInv[:,0], dataInv[:,1], 'w.')
mesh.plotSlice(actind, normal="Y", grid=True)
plt.plot(dataInv[:,0], dataInv[:,2], 'w.')
Out[22]:
In [23]:
zlocCC = Utils.mkvc(topoCC.reshape([-1,1]).repeat(mesh.nCz, axis=1))
depthweight = 1./ ((abs(mesh.gridCC[:,2] - zlocCC)+2.5)**1.)
depthweight /= depthweight.max()
depthweight[~actind] = np.nan
zlocCC[~actind] = np.nan
In [24]:
mesh.plotSlice(np.log10(depthweight), normal="Y", grid=True, clim=(-3, 0), ind=26)
Out[24]:
In [25]:
import simpegSP as SP
from pymatsolver import PardisoSolver
from SimPEG.EM.Static.SIP.Regularization import MultiRegularization
xyzM = xyzlocInv[:-1,:]
xyzN = np.atleast_2d(xyzlocInv[-1,:]).repeat(xyzM.shape[0], axis=0)
wires = Maps.Wires(('jsx', actMap.nP), ('jsy', actMap.nP), ('jsz', actMap.nP))
prb = SP.Problem_CC(mesh, sigma=sigma, jsxMap=actMap*wires.jsx, jsyMap=actMap*wires.jsy, jszMap=actMap*wires.jsz, Solver=PardisoSolver)
rx = SP.Rx.Dipole(xyzN, xyzM)
src = SP.Src.StreamingCurrents([rx], L=np.ones(mesh.nC)*1e-6, mesh=mesh, modelType="CurrentDensity")
survey = SP.Survey([src])
survey.pair(prb)
In [26]:
survey.dobs = dataInv[:-1,2]
In [29]:
from SimPEG import (Mesh, Maps, Utils, DataMisfit, Regularization,
Optimization, Inversion, InvProblem, Directives)
survey.std = 0.
survey.eps = abs(survey.dobs).max() * 0.05
dmisfit = DataMisfit.l2_DataMisfit(survey)
regmap = Maps.IdentityMap(nP = actMap.nP*3)
reg = MultiRegularization(mesh, mapping=regmap ,nModels=3, indActive=actind)
reg.cell_weights = depthweight[actind]
reg.alpha_s = 1.
reg.alpha_x = 1.
reg.alpha_y = 1.
reg.alpha_z = 1.
reg.ratios = [1., 1., 0.1]
opt = Optimization.InexactGaussNewton(maxIter=15, tolX=1e-20, tolF=1e-20)
opt.maxIterCG = 3
opt.maxIterLS = 30
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
target = Directives.TargetMisfit()
# Create an inversion object
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
betaest = Directives.BetaEstimate_ByEig()
betaest.beta0_ratio = 2.
invProb.beta = 1e0
inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
prb.counter = opt.counter = Utils.Counter()
opt.LSshorten = 0.5
opt.remember('xc')
m0 = np.ones(actMap.nP*3)*0.
reg.mref = m0
opt.maxIterCG = 10
mopt = inv.run(m0)
In [30]:
target.target
Out[30]:
In [31]:
fig = plt.figure(figsize = (5,5))
ax = plt.subplot(111)
vmin, vmax = invProb.dpred.min(), invProb.dpred.max()*1.3
Utils.plot2Ddata(xyzM[:,:2], invProb.dpred, ncontour=100, contourOpts={"cmap":"jet", "alpha":0.6, "vmin":vmin, "vmax":vmax}, dataloc=True, ax=ax)
Out[31]:
In [32]:
xc = opt.recall("xc")
In [44]:
dpred10 = survey.dpred(xc[10])
In [45]:
fig = plt.figure(figsize = (5,5))
ax = plt.subplot(111)
vmin, vmax = dpred10.min(), dpred10.max()*1.3
Utils.plot2Ddata(xyzM[:,:2], dpred10, ncontour=100, contourOpts={"cmap":"jet", "alpha":0.6, "vmin":vmin, "vmax":vmax}, dataloc=True, ax=ax)
Out[45]:
In [35]:
fig = plt.figure(figsize = (5,5))
ax = plt.subplot(111)
Utils.plot2Ddata(xyzM[:,:2], survey.dobs, ncontour=100, contourOpts={"cmap":"jet", "alpha":0.6, "vmin":vmin, "vmax":vmax}, dataloc=True, ax=ax)
Out[35]:
In [36]:
plt.plot(dpred10)
plt.plot(survey.dobs)
Out[36]:
In [46]:
iteration = 10
xc = opt.recall("xc")
temp = xc[iteration].reshape((actMap.nP, 3), order="F")
tempx = actMap * temp[:,0]
tempy = actMap * temp[:,1]
tempz = actMap * temp[:,2]
temp_amp = np.sqrt(temp[:,0]**2 + temp[:,1]**2 + temp[:,2]**2)
tempx[~actind] = np.nan
tempy[~actind] = np.nan
tempz[~actind] = np.nan
temp = np.c_[tempx, tempy, tempz]
In [47]:
temp_amp.min()
temp_amp.max()
Out[47]:
In [52]:
actMap_plot = Maps.InjectActiveCells(mesh, actind, np.nan)
In [53]:
indz, indy = 17, 31
mesh.plotSlice((actMap_plot*temp_amp), ind = indz, clim=(0, temp_amp.max()), normal="Z")
plt.ylim(-1500, 1500)
plt.xlim(-1500, 1500)
plt.gca().set_aspect("equal")
In [54]:
mesh.plotSlice((actMap_plot*temp_amp), clim=(0, temp_amp.max()), normal="Y", ind=indy)
plt.xlim(-1500, 1500)
plt.ylim(-300, 90)
plt.gca().set_aspect(2)
print mesh.vectorCCz[indz]
print mesh.vectorCCy[indy]
In [59]:
mesh.plotSlice(temp, clim=(0, temp_amp.max()), normal="Y", vType='CCv', view='vec', ind=32, streamOpts={"color":'w'})
print mesh.vectorCCz[indz]
plt.xlim(-1500, 1500)
plt.ylim(-300, 90)
plt.gca().set_aspect(2)
In [60]:
mesh.plotSlice(tempx, normal="X", grid=True, ind=19)
Out[60]:
In [61]:
mesh.writeVTK("pohang1108", models={"sigma":sigma, "amp":actMap_plot*temp_amp, "j":np.c_[tempx, tempy, tempz]})
# np.savetxt("topoxyz.csv", topoInv, delimiter=",", header="xcoord, ycoord, zcoord")
# np.savetxt("dataxyz.csv", out, delimiter=",", header="xcoord, ycoord, zcoord")
In [62]:
actMap_plot*temp_amp
Out[62]:
In [63]:
temp_amp.min()
Out[63]:
In [64]:
temp_amp.max()
Out[64]: