wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-12 12:30:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'dcip/dev' of https://github.com/simpeg/simpeg into dcip/dev

Browse Source
This commit is contained in:
seogi_macbook
2016-04-13 16:52:44 -07:00
parent fb1ff4e867 16c6cc8d74
commit 05e3b02b3a
2 changed files with 121 additions and 82 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/DCIPUtils.py
+73 -57
View File
@@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName):
return model
def plot_pseudoSection(DCsurvey, axs, stype, dtype="voltage",clim=None):
def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
"""
Read list of 2D tx-rx location and plot a speudo-section of apparent
resistivity.
@@ -179,7 +179,7 @@ def plot_pseudoSection(DCsurvey, axs, stype, dtype="voltage",clim=None):
Input:
:param d2D, z0
:switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
:switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
Output:
:figure scatter plot overlayed on image
@@ -221,23 +221,43 @@ def plot_pseudoSection(DCsurvey, axs, stype, dtype="voltage",clim=None):
Cmid = (Tx[0][0] + Tx[1][0])/2
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
# Compute pant leg of apparent rho
if stype == 'pdp':
leg = data * 2*np.pi * MA * ( MA + MN ) / MN
# Change output for dtype
if dtype == 'volt':
leg = np.log10(abs(1/leg))
rho = np.hstack([rho,data])
elif stype == 'dpdp':
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
else:
# Compute pant leg of apparent rho
if stype == 'pdp':
leg = data * 2*np.pi * MA * ( MA + MN ) / MN
elif stype == 'dpdp':
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
else:
print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
break
if dtype == 'appc':
leg = np.log10(abs(1./leg))
rho = np.hstack([rho,leg])
elif dtype == 'appr':
leg = np.log10(abs(leg))
rho = np.hstack([rho,leg])
else:
print """dtype must be 'appr' | 'appc' | 'volt' """
break
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
#TODO ... let stick to list then finally convert to array.
if dtype =="voltage":
rho = np.hstack([rho,leg])
elif dtype =="appr":
rho = np.hstack([rho,data])
ax = axs
@@ -250,28 +270,33 @@ def plot_pseudoSection(DCsurvey, axs, stype, dtype="voltage",clim=None):
else:
vmin, vmax = clim[0], clim[1]
plt.imshow(grid_rho.T, extent = (np.min(midx),np.max(midx),np.min(midz),np.max(midz)), origin='lower', alpha=0.8, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
cbar = plt.colorbar(format = '%.2f',fraction=0.04,orientation="horizontal")
grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax))
cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
cmin,cmax = cbar.get_clim()
ticks = np.linspace(cmin,cmax,3)
cbar.set_ticks(ticks)
cbar.ax.tick_params(labelsize=10)
if dtype == 'appc':
cbar.set_label("App.Cond",size=12)
elif dtype == 'appr':
cbar.set_label("App.Res.",size=12)
elif dtype == 'volt':
cbar.set_label("Potential (V)",size=12)
# Plot apparent resistivity
plt.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
#ax.set_xticklabels([])
#ax.set_yticklabels([])
ax.set_xticklabels([])
ax.set_yticklabels([])
#ax.set_ylabel('Z')
#ax.yaxis.tick_right()
#ax.yaxis.set_label_position('right')
plt.gca().set_aspect('equal', adjustable='box')
return ax
return ph
def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
"""
@@ -371,16 +396,6 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
SrcList.append(srcClass)
#==============================================================================
# elif re.match(stype,'dpdp'):
#
# for ii in range(0, int(nstn)-2):
#
# indx = np.min([ii+n+1,nstn])
# Tx.append(np.c_[M[ii,:],N[ii,:]])
# Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]])
#==============================================================================
elif stype == 'gradient':
# Gradient survey only requires Tx at end of line and creates a square
@@ -523,22 +538,22 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
fid.close()
def convertObs_DC3D_to_2D(DCsurvey,lineID,flag):
def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
"""
Read DC survey and data and change
coordinate system to distance along line assuming
all data is acquired along line.
First transmitter pole is assumed to be at the origin
Read DC survey and projects the coordinate system
according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
In the 'local' system, station coordinates are referenced
to distance from the first srcLoc[0].loc[0]
Assumes flat topo for now...
The Z value is preserved, but Y coordinates zeroed.
Input:
:param Tx, Rx
:param survey3D
Output:
:figure Tx2d, Rx2d
:figure survey2D
Edited Feb 17th, 2016
Edited April 6th, 2016
@author: dominiquef
@@ -584,19 +599,19 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID,flag):
# Find A electrode along line
vec, r = r_unit(x0,Tx[ii][0,0:2])
A = stn_id(vecTx,vec,r)
# Find B electrode along line
vec, r = r_unit(x0,Tx[ii][1,0:2])
B = stn_id(vecTx,vec,r)
M = np.zeros(nrx)
N = np.zeros(nrx)
for kk in range(nrx):
# Find all M electrodes along line
vec, r = r_unit(x0,Rx[0][kk,0:2])
M[kk] = stn_id(vecTx,vec,r)
# Find all N electrodes along line
vec, r = r_unit(x0,Rx[1][kk,0:2])
N[kk] = stn_id(vecTx,vec,r)
@@ -606,14 +621,14 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID,flag):
B = Tx[ii][1,1]
M = Rx[0][:,1]
N = Rx[1][:,1]
elif flag == 'Xloc':
""" Copy the rx-tx locs"""
A = Tx[ii][0,0]
B = Tx[ii][1,0]
M = Rx[0][:,0]
N = Rx[1][:,0]
Rx = DC.RxDipole(np.c_[M,np.zeros(nrx),Rx[0][:,2]],np.c_[N,np.zeros(nrx),Rx[1][:,2]])
srcLists.append( DC.SrcDipole( [Rx], np.asarray([A,0,Tx[ii][0,2]]),np.asarray([B,0,Tx[ii][1,2]]) ) )
@@ -628,16 +643,16 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID,flag):
def readUBC_DC3Dobs(fileName):
"""
Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
Read UBC GIF DCIP 3D observation file and generate survey
Input:
:param fileName, path to the UBC GIF 3D obs file
Output:
:param rx, tx, d, wd
:param DCIPsurvey
:return
Created on Mon December 7th, 2015
Created on Mon April 6th, 2015
@author: dominiquef
@@ -712,6 +727,7 @@ def readUBC_DC3Dobs(fileName):
def readUBC_DC2Dobs(fileName):
"""
------- NEEDS TO BE UPDATED ------
Read UBC GIF 2D observation file and generate arrays for tx-rx location
Input:
@@ -761,7 +777,7 @@ def readUBC_DC2Dobs(fileName):
def readUBC_DC2Dpre(fileName):
"""
Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location
Input:
:param fileName, path to the UBC GIF 3D obs file
@@ -804,15 +820,15 @@ def readUBC_DC2Dpre(fileName):
else:
tx = np.r_[temp[0],np.nan,temp[1],temp[2],np.nan,temp[3]]
if zflag:
rx = np.c_[temp[4],np.nan,temp[5],temp[6],np.nan,temp[7]]
else:
rx = np.c_[temp[2],np.nan,np.nan,temp[3],np.nan,np.nan]
# Check if there is data with the location
d.append(temp[-1])
@@ -825,7 +841,7 @@ def readUBC_DC2Dpre(fileName):
survey.dobs = np.asarray(d)
return {'DCsurvey':survey}
def readUBC_DC2DMesh(fileName):
"""
Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
SimPEG/Examples/DC_Forward_PseudoSection.py
+48 -25
View File
@@ -2,19 +2,27 @@ from SimPEG import Mesh, Utils, np, sp
import SimPEG.DCIP as DC
import time
def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
"""
DC Forward Simulation
=====================
Forward model conductive spheres in a half-space and plot a pseudo-section
Forward model two conductive spheres in a half-space and plot a
pseudo-section. Assumes an infinite line source and measures along the
center of the spheres.
INPUT:
loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
radi = Radius of spheres [r1,r2]
param = Conductivity of background and two spheres [m0,m1,m2]
stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
Created by @fourndo on Mon Feb 01 19:28:06 2016
"""
assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)"
assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)"
if loc is None:
loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
@@ -27,7 +35,6 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
# First we need to create a mesh and a model.
# This is our mesh
dx = 5.
@@ -52,14 +59,10 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
# Get index of the center
indy = int(mesh.nCy/2)
# Plot the model for reference
# Define core mesh extent
xlim = 200
zlim = 125
# Specify the survey type: "pdp" | "dpdp"
zlim = 100
# Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
ends = [(-175,0),(175,0)]
@@ -77,12 +80,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
dl_y = ( Tx[-1][1,1] - Tx[0][1,0] ) / dl_len
azm = np.arctan(dl_y/dl_x)
#azm = np.arctan(dl_y/dl_x)
#Set boundary conditions
mesh.setCellGradBC('neumann')
# Define the differential operators needed for the DC problem
# Define the linear system needed for the DC problem. We assume an infitite
# line source for simplicity.
Div = mesh.faceDiv
Grad = mesh.cellGrad
Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
@@ -145,16 +149,23 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
print 'Forward completed'
# Let's just convert the 3D format into 2D (distance along line) and plot
# [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
survey2D.dobs =np.hstack(data)
# Here is an example for the first tx-rx array
if plotIt:
import matplotlib.pyplot as plt
fig = plt.figure()
fig = plt.figure(figsize=(7,7))
ax = plt.subplot(2,1,1, aspect='equal')
mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
# Plot the location of the spheres for reference
circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
ax.add_artist(circle1)
ax.add_artist(circle2)
dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y',
ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()]))
ax.set_title('3-D model')
plt.gca().set_aspect('equal', adjustable='box')
plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
@@ -163,22 +174,34 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
ax = plt.subplot(2,1,2, aspect='equal')
pos = ax.get_position()
ax.set_position([pos.x0 , pos.y0 + 0.025 , pos.width, pos.height])
pos = ax.get_position()
cbarax = fig.add_axes([pos.x0 , pos.y0 + 0.025 , pos.width, pos.height * 0.04]) ## the parameters are the specified position you set
cb = fig.colorbar(dat[0],cax=cbarax, orientation="horizontal",
ax = ax, ticks=np.linspace(np.log10(sig.min()),
np.log10(sig.max()), 3), format="$10^{%.1f}$")
cb.set_label("Conductivity (S/m)",size=12)
cb.ax.tick_params(labelsize=12)
# Second plot for the predicted apparent resistivity data
ax2 = plt.subplot(2,1,2, aspect='equal')
# Plot the location of the spheres for reference
circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
ax.add_artist(circle1)
ax.add_artist(circle2)
circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
ax2.add_artist(circle1)
ax2.add_artist(circle2)
# Add the speudo section
DC.plot_pseudoSection(survey2D,ax,stype)
dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype)
# plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
# plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
# plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
ax2.set_title('Apparent Conductivity data')
plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
plt.show()
return fig, ax