From f799733a9d9c59961b34e5d6dfcb558daed3bc60 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 6 Apr 2016 14:43:33 -0700 Subject: [PATCH] Fix TwoSphere example and Utils.pseudoPlot --- SimPEG/DCIP/DCIPUtils.py | 50 ++++++++++--------- SimPEG/Examples/DC_Forward_PseudoSection.py | 55 ++++++++++++++------- 2 files changed, 64 insertions(+), 41 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index d8cf04bf..af94ad4e 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="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 diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index 2f198238..904dcbd3 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -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