Update speudo plot and allow app_res, app_con, volt

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
+        # Change output for dtype
-        if stype == 'pdp':
+        if dtype == 'volt':
-            leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
 
-            leg = np.log10(abs(1/leg))
+            rho = np.hstack([rho,data])
 
-        elif stype == 'dpdp':
+        else:
-            leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
 
-            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

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()