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

Fix TwoSphere example and Utils.pseudoPlot

Browse Source
This commit is contained in:
D Fournier
2016-04-06 14:43:33 -07:00
parent 8b94cd4dfe
commit f799733a9d
2 changed files with 64 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/DCIPUtils.py
+26 -24
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, dtype="appr",clim=None):
"""
Read list of 2D tx-rx location and plot a speudo-section of apparent
resistivity.
@@ -230,12 +230,14 @@ def plot_pseudoSection(DCsurvey, axs, stype, dtype="voltage",clim=None):
elif stype == 'dpdp':
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
leg = np.log10(abs(1/leg))
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":
if dtype =="appr":
rho = np.hstack([rho,leg])
elif dtype =="appr":
elif dtype =="voltage":
rho = np.hstack([rho,data])
@@ -250,28 +252,27 @@ 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)
# Plot apparent resistivity
plt.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
cbar.set_label("App. Conductivity",size=12)
ax.set_xticklabels([])
ax.set_yticklabels([])
# Plot apparent resistivity
ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
#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):
"""
@@ -515,12 +516,12 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
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 survey3D
@@ -528,7 +529,7 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
Output:
:figure survey2D
Edited Feb 17th, 2016
Edited April 6th, 2016
@author: dominiquef
@@ -618,16 +619,16 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
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
@@ -702,6 +703,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:
@@ -751,7 +753,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
SimPEG/Examples/DC_Forward_PseudoSection.py
+38 -17
View File
@@ -56,7 +56,7 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
# Plot the model for reference
# Define core mesh extent
xlim = 200
zlim = 125
zlim = 100
# Specify the survey type: "pdp" | "dpdp"
@@ -77,7 +77,7 @@ 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')
@@ -146,39 +146,60 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
# 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')
plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y')
plt.xlim([-xlim,xlim])
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)
# 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')
ax2.set_title('Apparent Conductivity data')
plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
plt.show()
return fig, ax