From 16c6cc8d74ab7ce4e61980cfd009cfedab44dc8f Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 6 Apr 2016 22:17:34 -0700 Subject: [PATCH] Update speudo plot and allow app_res, app_con, volt --- SimPEG/DCIP/DCIPUtils.py | 80 +++++++++++++-------- SimPEG/Examples/DC_Forward_PseudoSection.py | 54 +++++++------- 2 files changed, 80 insertions(+), 54 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index af94ad4e..9b637d82 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName): return model -def plot_pseudoSection(DCsurvey, axs, stype, dtype="appr",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="appr",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,25 +221,43 @@ def plot_pseudoSection(DCsurvey, axs, stype, dtype="appr",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: - leg = np.log10(abs(1/leg)) + # 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 =="appr": - rho = np.hstack([rho,leg]) - elif dtype =="voltage": - rho = np.hstack([rho,data]) - ax = axs @@ -260,14 +278,20 @@ def plot_pseudoSection(DCsurvey, axs, stype, dtype="appr",clim=None): ticks = np.linspace(cmin,cmax,3) cbar.set_ticks(ticks) cbar.ax.tick_params(labelsize=10) - cbar.set_label("App. Conductivity",size=12) + + 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 ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax)) #ax.set_xticklabels([]) #ax.set_yticklabels([]) - + plt.gca().set_aspect('equal', adjustable='box') @@ -516,7 +540,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): """ - Read DC survey and projects the coordinate system + 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] @@ -575,19 +599,19 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): # 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) @@ -597,14 +621,14 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): 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]]) ) ) @@ -796,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]) @@ -817,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 diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index 904dcbd3..e1763da6 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -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,15 +59,11 @@ 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 = 100 - # Specify the survey type: "pdp" | "dpdp" - - # 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)] ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]] @@ -82,7 +85,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): #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,10 +149,9 @@ 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) , '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(figsize=(7,7)) @@ -158,29 +161,29 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): 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]) - - - pos = ax.get_position() + + + 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 + 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) - + cb.ax.tick_params(labelsize=12) + # Second plot for the predicted apparent resistivity data ax2 = plt.subplot(2,1,2, aspect='equal') @@ -189,16 +192,15 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): 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 - dat = DC.plot_pseudoSection(survey2D,ax2,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') - ax2.set_title('Apparent Conductivity data') - + ax2.set_title('Apparent Conductivity data') + plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) plt.show()