wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-09-09 11:34:26 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'dev' into ref/dev

Browse Source 
Conflicts:
	SimPEG/DCIP/DCIPUtils.py
This commit is contained in:
D Fournier
2016-05-27 10:01:28 -07:00
parent 0379df2bf2 e5ec512517
commit 7b72d3a92d
10 changed files with 386 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/DCIPUtils.py
+85 -38
View File
@@ -169,7 +169,8 @@ def readUBC_DC2DModel(fileName):
return model
def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None):
def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
"""
Read list of 2D tx-rx location and plot a speudo-section of apparent
resistivity.
@@ -187,9 +188,6 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
from scipy.interpolate import griddata
import pylab as plt
# Set depth to 0 for now
z0 = 0.
# Pre-allocate
midx = []
midz = []
@@ -254,38 +252,54 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
ax = axs
# Grid points
grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
# Scale the color scheme
if clim == None:
vmin, vmax = rho.min(), rho.max()
else:
vmin, vmax = clim[0], clim[1]
# Plot data
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 unitType == 'appConductivity':
cbar.set_label("App.Cond",size=12)
elif unitType == 'appResistivity':
cbar.set_label("App.Res.",size=12)
elif unitType == 'volt':
cbar.set_label("Potential (V)",size=12)
# Plot apparent resistivity
ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax)
plt.gca().tick_params(axis='both', which='major', labelsize=8)
#ax.set_xticklabels([])
#ax.set_yticklabels([])
if contour is not None:
plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax)
# Add scatter points
axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax)
if colorbar:
if unitType == 'volt':
cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
else:
cbar = plt.colorbar(ph, ax = axs, 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 unitType == 'appConductivity':
cbar.set_label("App.Cond",size=12)
elif unitType == 'appResistivity':
cbar.set_label("App.Res.",size=12)
elif unitType == 'volt':
cbar.set_label("Potential (V)",size=12)
if not axlabel:
axs.set_xticklabels([])
axs.set_yticklabels([])
plt.gca().set_aspect('equal', adjustable='box')
@@ -440,7 +454,8 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
survey = DC.SurveyDC(SrcList)
return survey, Tx, Rx
def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType):
def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
"""
Write UBC GIF DCIP 2D or 3D observation file
@@ -460,6 +475,12 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType):
fid = open(fileName,'w')
fid.write('! ' + surveyType + ' FORMAT\n')
if iptype!=0:
fid.write('IPTYPE=%i\n'%iptype)
else:
fid.write('! ' + stype + ' FORMAT\n')
count = 0
for ii in range(DCsurvey.nSrc):
@@ -491,18 +512,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType):
if surveyType == 'SURFACE':
fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
M = M[:,0]
N = N[:,0]
if surveyType == 'GENERAL':
# Flip sign for z-elevation to depth
tx[2::2,:] = -tx[2::2,:]
fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
M = M[:,0::2]
N = N[:,0::2]
# Flip sign for z-elevation to depth
M[:,1::2] = -M[:,1::2]
N[:,1::2] = -N[:,1::2]
fid.write('%i\n'% nD)
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
if dim=='3D':
@@ -514,10 +542,11 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType):
if surveyType == 'GENERAL':
fid.writelines("%e " % ii for ii in mkvc(tx))
fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
fid.write('%i\n'% nD)
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
fid.write('\n')
count += nD
@@ -620,43 +649,53 @@ def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'):
return DCsurvey2D
def readUBC_DC3Dobs(fileName):
def readUBC_DC3Dobs(fileName, rtype = 'DC'):
"""
Read UBC GIF DCIP 3D observation file and generate survey
Read UBC GIF IP 3D observation file and generate survey
:param string fileName:, path to the UBC GIF 3D obs file
:rtype: Survey
:return: DCIPsurvey
"""
zflag = True # Flag for z value provided
# Load file
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
if rtype == 'IP':
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
elif rtype == 'DC':
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
else:
print "rtype must be 'DC'(default) | 'IP'"
# Pre-allocate
srcLists = []
Rx = []
d = []
wd = []
zflag = True # Flag for z value provided
# Countdown for number of obs/tx
count = 0
for ii in range(obsfile.shape[0]):
# Skip if blank line
if not obsfile[ii]:
continue
# First line is transmitter with number of receivers
# First line or end of a transmitter block, read transmitter info
if count==0:
temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
# Read the line
temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T)
count = int(temp[-1])
# Check if z value is provided, if False -> nan
if len(temp)==5:
tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
zflag = False
tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
zflag = False # Pass on the flag to the receiver loc
else:
tx = temp[:-1]
@@ -664,8 +703,16 @@ def readUBC_DC3Dobs(fileName):
rx = []
continue
temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
# Filter out negative IP
# if temp[-2] < 0:
# count = count -1
# print "Negative!"
#
# else:
# If the Z-location is provided, otherwise put nan
if zflag:
rx.append(temp[:-2])
@@ -675,7 +722,7 @@ def readUBC_DC3Dobs(fileName):
wd.append(temp[-1])
else:
rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
# Check if there is data with the location
if len(temp)==6:
d.append(temp[-2])
@@ -683,7 +730,7 @@ def readUBC_DC3Dobs(fileName):
count = count -1
# Reach the end of transmitter block
# Reach the end of transmitter block, append the src, rx and continue
if count == 0:
rx = np.asarray(rx)
Rx = DC.RxDipole(rx[:,:3],rx[:,3:])