from SimPEG import np from SimPEG.EM.Static import DC, IP 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. Assumes flat topo for now... 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 Edited Feb 17th, 2016 @author: dominiquef """ from SimPEG import np from scipy.interpolate import griddata import pylab as plt # Set depth to 0 for now z0 = 0. # Pre-allocate midx = [] midz = [] rho = [] LEG = [] count = 0 # Counter for data for ii in range(DCsurvey.nSrc): Tx = DCsurvey.srcList[ii].loc Rx = DCsurvey.srcList[ii].rxList[0].locs nD = DCsurvey.srcList[ii].rxList[0].nD data = DCsurvey.dobs[count:count+nD] count += nD # Get distances between each poles A-B-M-N if stype == 'pdp': MA = np.abs(Tx[0] - Rx[0][:,0]) NA = np.abs(Tx[0] - Rx[1][:,0]) MN = np.abs(Rx[1][:,0] - Rx[0][:,0]) # Create mid-point location Cmid = Tx[0] Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 if DCsurvey.mesh.dim == 2: zsrc = Tx[1] elif DCsurvey.mesh.dim ==3: zsrc = Tx[2] elif stype == 'dpdp': MA = np.abs(Tx[0][0] - Rx[0][:,0]) MB = np.abs(Tx[1][0] - Rx[0][:,0]) NA = np.abs(Tx[0][0] - Rx[1][:,0]) NB = np.abs(Tx[1][0] - Rx[1][:,0]) # Create mid-point location Cmid = (Tx[0][0] + Tx[1][0])/2 Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 if DCsurvey.mesh.dim == 2: zsrc = (Tx[0][1] + Tx[1][1])/2 elif DCsurvey.mesh.dim ==3: zsrc = (Tx[0][2] + Tx[1][2])/2 # Change output for dtype if dtype == 'volt': rho = np.hstack([rho,data]) 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 ) LEG.append(1./(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 ]) if DCsurvey.mesh.dim==3: midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ]) elif DCsurvey.mesh.dim==2: midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ]) 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') if clim == None: vmin, vmax = rho.min(), rho.max() else: vmin, vmax = clim[0], clim[1] 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), vmin=vmin, vmax=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 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') return ph, LEG def gen_DCIPsurvey(endl, mesh, stype, a, b, n): """ Load in endpoints and survey specifications to generate Tx, Rx location stations. Assumes flat topo for now... Input: :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] :object mesh -> SimPEG mesh object :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient' : param a, n -> pole seperation, number of rx dipoles per tx Output: :param Tx, Rx -> List objects for each tx location Lines: P1x, P1y, P1z, P2x, P2y, P2z Created on Wed December 9th, 2015 @author: dominiquef !! Require clean up to deal with DCsurvey """ from SimPEG import np def xy_2_r(x1,x2,y1,y2): r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) ) return r ## Evenly distribute electrodes and put on surface # Mesure survey length and direction dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1]) dl_x = ( endl[1,0] - endl[0,0] ) / dl_len dl_y = ( endl[1,1] - endl[0,1] ) / dl_len nstn = np.floor( dl_len / a ) # Compute discrete pole location along line stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a if mesh.dim==2: ztop = mesh.vectorNy[-1] # Create line of P1 locations M = np.c_[stn_x, np.ones(nstn).T*ztop] # Create line of P2 locations N = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop] elif mesh.dim==3: ztop = mesh.vectorNz[-1] # Create line of P1 locations M = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop] # Create line of P2 locations N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop] ## Build list of Tx-Rx locations depending on survey type # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn] # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B] SrcList = [] if stype != 'gradient': for ii in range(0, int(nstn)-1): if stype == 'dpdp': tx = np.c_[M[ii,:],N[ii,:]] elif stype == 'pdp': tx = np.c_[M[ii,:],M[ii,:]] # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]]) # Current elctrode seperation AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1]) # Number of receivers to fit nstn = np.min([np.floor( (AB - b) / a ) , n]) # Check if there is enough space, else break the loop if nstn <= 0: continue # Compute discrete pole location along line stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a # Create receiver poles if mesh.dim==3: # Create line of P1 locations P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop] # Create line of P2 locations P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop] rxClass = DC.Rx.Dipole(P1, P2) elif mesh.dim==2: # Create line of P1 locations P1 = np.c_[stn_x, np.ones(nstn).T*ztop] # Create line of P2 locations P2 = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop] rxClass = DC.Rx.Dipole_ky(P1, P2) if stype == 'dpdp': srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:]) elif stype == 'pdp': srcClass = DC.Src.Pole([rxClass], M[ii,:]) SrcList.append(srcClass) elif stype == 'gradient': # Gradient survey only requires Tx at end of line and creates a square # grid of receivers at in the middle at a pre-set minimum distance # Get the edge limit of survey area min_x = endl[0,0] + dl_x * b min_y = endl[0,1] + dl_y * b max_x = endl[1,0] - dl_x * b max_y = endl[1,1] - dl_y * b box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 ) box_w = box_l/2. nstn = np.floor( box_l / a ) # Compute discrete pole location along line stn_x = min_x + np.array(range(int(nstn)))*dl_x*a stn_y = min_y + np.array(range(int(nstn)))*dl_y*a # Define number of cross lines nlin = int(np.floor( box_w / a )) lind = range(-nlin,nlin+1) ngrad = nstn * len(lind) rx = np.zeros([ngrad,6]) for ii in range( len(lind) ): # Move line in perpendicular direction by dipole spacing lxx = stn_x - lind[ii]*a*dl_y lyy = stn_y + lind[ii]*a*dl_x M = np.c_[ lxx, lyy , np.ones(nstn).T*ztop] N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*ztop] rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N] if mesh.dim==3: rxClass = DC.Rx.Dipole(rx[:,:3], rx[:,3:]) elif mesh.dim==2: M = M[:,[0,2]] N = N[:,[0,2]] rxClass = DC.Rx.Dipole_ky(rx[:,[0,2]], rx[:,[3,5]]) srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:]) SrcList.append(srcClass) else: print """stype must be either 'pdp', 'dpdp' or 'gradient'. """ return SrcList def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0): """ Write UBC GIF DCIP 2D or 3D observation file Input: :string fileName -> including path where the file is written out :DCsurvey DC survey class object :string dtype -> either '2D' | '3D' :string stype -> either 'SURFACE' | 'GENERAL' Output: :param UBC2D-Data file :return Last edit: February 16th, 2016 @author: dominiquef """ from SimPEG import mkvc assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'" assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" fid = open(fileName,'w') if iptype!=0: fid.write('IPTYPE=%i\n'%iptype) else: fid.write('! ' + stype + ' FORMAT\n') count = 0 for ii in range(DCsurvey.nSrc): tx = np.c_[DCsurvey.srcList[ii].loc] rx = DCsurvey.srcList[ii].rxList[0].locs nD = DCsurvey.srcList[ii].nD M = rx[0] N = rx[1] # Adapt source-receiver location for dtype and stype if dtype=='2D': if stype == 'SIMPLE': #fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) A = np.repeat(tx[0,0],M.shape[0],axis=0) B = np.repeat(tx[0,1],M.shape[0],axis=0) M = M[:,0] N = N[:,0] np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') else: if stype == 'SURFACE': fid.writelines("%f " % ii for ii in mkvc(tx[0,:])) M = M[:,0] N = N[:,0] if stype == '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='%f',delimiter=' ',newline='\n') if dtype=='3D': if stype == 'SURFACE': fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:])) M = M[:,0:2] N = N[:,0:2] if stype == 'GENERAL': 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 fid.close()