From a0174e4f3043216cbc302be4fd775cda2674cf95 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 29 Apr 2016 12:52:45 -0700 Subject: [PATCH 1/7] kwarg name updates --- SimPEG/DCIP/DCIPUtils.py | 84 ++++++++++----------- SimPEG/Examples/DC_Forward_PseudoSection.py | 32 ++++---- 2 files changed, 55 insertions(+), 61 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index e94b930f..91eab0b3 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='dpdp', dtype="appc", clim=None): +def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None): """ Read list of 2D tx-rx location and plot a speudo-section of apparent resistivity. @@ -177,16 +177,12 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): 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) + :param DCsurvey, z0 + :switch surveyType -> Either 'pole-dipole' | 'dipole-dipole' + :switch unitType=-> Either 'appResistivity' | 'appConductivity' | 'volt' Output: :figure scatter plot overlayed on image - Edited Feb 17th, 2016 - - @author: dominiquef - """ from SimPEG import np from scipy.interpolate import griddata @@ -221,39 +217,39 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): Cmid = (Tx[0][0] + Tx[1][0])/2 Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 - # Change output for dtype - if dtype == 'volt': + # Change output for unitType + if unitType == 'volt': rho = np.hstack([rho,data]) else: # Compute pant leg of apparent rho - if stype == 'pdp': + if surveyType == 'pole-dipole': leg = data * 2*np.pi * MA * ( MA + MN ) / MN - elif stype == 'dpdp': + elif surveyType == 'dipole-dipole': leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA ) else: - print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """ + print """unitType must be 'pole-dipole' | 'dipole-dipole' """ break - if dtype == 'appc': + if unitType == 'appConductivity': leg = np.log10(abs(1./leg)) rho = np.hstack([rho,leg]) - elif dtype == 'appr': + elif unitType == 'appResistivity': leg = np.log10(abs(leg)) rho = np.hstack([rho,leg]) else: - print """dtype must be 'appr' | 'appc' | 'volt' """ + print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """ break midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) @@ -278,12 +274,12 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): ticks = np.linspace(cmin,cmax,3) cbar.set_ticks(ticks) cbar.ax.tick_params(labelsize=10) - - if dtype == 'appc': + + if unitType == 'appConductivity': cbar.set_label("App.Cond",size=12) - elif dtype == 'appr': + elif unitType == 'appResistivity': cbar.set_label("App.Res.",size=12) - elif dtype == 'volt': + elif unitType == 'volt': cbar.set_label("Potential (V)",size=12) # Plot apparent resistivity @@ -298,7 +294,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): return ph -def gen_DCIPsurvey(endl, mesh, stype, a, b, n): +def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): """ Load in endpoints and survey specifications to generate Tx, Rx location stations. @@ -308,7 +304,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): 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' + :switch surveyType -> 'dipole-dipole' | 'pole-dipole' | 'gradient' : param a, n -> pole seperation, number of rx dipoles per tx Output: @@ -354,14 +350,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): SrcList = [] - if stype != 'gradient': + if surveyType != 'gradient': for ii in range(0, int(nstn)-1): - if stype == 'dpdp': + if surveyType == 'dipole-dipole': tx = np.c_[M[ii,:],N[ii,:]] - elif stype == 'pdp': + elif surveyType == 'pole-dipole': tx = np.c_[M[ii,:],M[ii,:]] # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]]) @@ -390,13 +386,13 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): Rx.append(np.c_[P1,P2]) rxClass = DC.RxDipole(P1, P2) Tx.append(tx) - if stype == 'dpdp': + if surveyType == 'dipole-dipole': srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:]) - elif stype == 'pdp': + elif surveyType == 'pole-dipole': srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:]) SrcList.append(srcClass) - elif stype == 'gradient': + elif surveyType == '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 @@ -443,20 +439,20 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:]) SrcList.append(srcClass) else: - print """stype must be either 'pdp', 'dpdp' or 'gradient'. """ + print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """ survey = DC.SurveyDC(SrcList) return survey, Tx, Rx -def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): +def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): """ 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' + :string dim -> either '2D' | '3D' + :string surveyType -> either 'SURFACE' | 'GENERAL' Output: :param UBC2D-Data file @@ -469,11 +465,11 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): """ 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'" + assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'" + assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" fid = open(fileName,'w') - fid.write('! ' + stype + ' FORMAT\n') + fid.write('! ' + surveyType + ' FORMAT\n') count = 0 @@ -488,10 +484,10 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): M = rx[0] N = rx[1] - # Adapt source-receiver location for dtype and stype - if dtype=='2D': + # Adapt source-receiver location for dim and surveyType + if dim=='2D': - if stype == 'SIMPLE': + if surveyType == 'SIMPLE': #fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) A = np.repeat(tx[0,0],M.shape[0],axis=0) @@ -504,13 +500,13 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): else: - if stype == 'SURFACE': + if surveyType == 'SURFACE': fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) M = M[:,0] N = N[:,0] - if stype == 'GENERAL': + if surveyType == 'GENERAL': fid.writelines("%e " % ii for ii in mkvc(tx[::2,:])) M = M[:,0::2] @@ -519,15 +515,15 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): 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') - if dtype=='3D': + if dim=='3D': - if stype == 'SURFACE': + if surveyType == 'SURFACE': fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:])) M = M[:,0:2] N = N[:,0:2] - if stype == 'GENERAL': + if surveyType == 'GENERAL': fid.writelines("%e " % ii for ii in mkvc(tx)) @@ -538,7 +534,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): fid.close() -def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): +def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'): """ Read DC survey and projects the coordinate system according to the flag = 'Xloc' | 'Yloc' | 'local' (default) diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index 53467826..240cbb33 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -2,7 +2,7 @@ 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', dtype='appc', plotIt=True): +def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True): """ DC Forward Simulation ===================== @@ -15,14 +15,14 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p 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) + surveyType = survey type 'pole-dipole' or 'dipole-dipole' + unitType = Data type "appResistivity" | "appConductivity" | "volt" Created by @fourndo """ - 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)" + assert surveyType in ['pole-dipole', 'dipole-dipole'], "Source type (surveyType) must be pdp or dpdp (pole dipole or dipole dipole)" + assert unitType in ['appResistivity', 'appConductivity', 'volt'], "Unit type (unitType) must be appResistivity or appConductivity or volt (potential)" if loc is None: loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]] @@ -73,8 +73,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]] # We will handle the geometry of the survey for you and create all the combination of tx-rx along line - # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) - survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) + # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2]) + survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2]) # Define some global geometry dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) @@ -118,8 +118,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p rxloc_N = np.asarray(Rx[ii][:,3:]) - # For usual cases "dpdp" or "gradient" - if stype == 'pdp': + # For usual cases 'dipole-dipole' or "gradient" + if surveyType == 'pole-dipole': # Create an "inifinity" pole tx = np.squeeze(Tx[ii][:,0:1]) tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2 @@ -157,12 +157,12 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p fig = plt.figure(figsize=(7,7)) ax = plt.subplot(2,1,1, aspect='equal') # 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) + 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', + 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') @@ -188,15 +188,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p ax2 = plt.subplot(2,1,2, aspect='equal') # 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) + 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 - dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype) - - # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v') + dat = DC.plot_pseudoSection(survey2D, ax2, surveyType=surveyType, unitType=unitType) # 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') From 4257ea77b357bdab1811ac8f9adfb82845a6cf55 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 29 Apr 2016 15:09:04 -0700 Subject: [PATCH 2/7] remove InjectActiveCellsTopo. you should use InjectActiveCells --- SimPEG/Maps.py | 77 -------------------------------------------------- 1 file changed, 77 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index f17d2314..3e97499f 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -533,83 +533,6 @@ class ActiveCells(InjectActiveCells): FutureWarning) InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC) -class InjectActiveCellsTopo(IdentityMap): - """ - Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh) - - """ - - indActive = None #: Active Cells - valInactive = None #: Values of inactive Cells - nC = None #: Number of cells in the full model - - def __init__(self, mesh, indActive, nC=None): - self.mesh = mesh - - self.nC = nC or mesh.nC - - if indActive.dtype is not bool: - z = np.zeros(self.nC,dtype=bool) - z[indActive] = True - indActive = z - self.indActive = indActive - - self.indInactive = np.logical_not(indActive) - inds = np.nonzero(self.indActive)[0] - self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP)) - - @property - def shape(self): - return (self.nC, self.nP) - - @property - def nP(self): - """Number of parameters in the model.""" - return self.indActive.sum() - - def _transform(self, m): - val_temp = np.zeros(self.mesh.nC) - val_temp[self.indActive] = m - valInactive = np.zeros(self.mesh.nC) - #1D - if self.mesh.dim == 1: - z_temp = self.mesh.gridCC - val_temp[~self.indActive] = val_temp[np.argmax(z_temp[self.indActive])] - #2D - elif self.mesh.dim == 2: - act_temp = self.indActive.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - val_temp = val_temp.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - y_temp = self.mesh.gridCC[:,1].reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - for i in range(self.mesh.nCx): - act_tempx = act_temp[i,:] == 1 - val_temp[i,~act_tempx] = val_temp[i,np.argmax(y_temp[i,act_tempx])] - valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive] - #3D - elif self.mesh.dim == 3: - act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - for i in range(self.mesh.nCx*self.mesh.nCy): - act_tempxy = act_temp[i,:] == 1 - val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])] - valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive] - - self.valInactive = valInactive - - return self.P*m + self.valInactive - - def inverse(self, D): - return self.P.T*D - - def deriv(self, m): - return self.P - -class ActiveCellsTopo(InjectActiveCellsTopo): - def __init__(self, mesh, indActive, valInactive, nC=None): - warnings.warn( - "`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead", - FutureWarning) - InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC) class Weighting(IdentityMap): """ From ba8f270b3adb7d9aa6553b3ffac739ed7fef0cba Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 1 May 2016 13:17:16 -0700 Subject: [PATCH 3/7] start of surface2ind_topo --- SimPEG/Examples/Utils_surface2ind_topo.py | 42 +++++++++++++++ SimPEG/Utils/__init__.py | 1 + SimPEG/Utils/modelutils.py | 63 +++++++++++++++++++++++ 3 files changed, 106 insertions(+) create mode 100644 SimPEG/Examples/Utils_surface2ind_topo.py create mode 100644 SimPEG/Utils/modelutils.py diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py new file mode 100644 index 00000000..61a0f937 --- /dev/null +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -0,0 +1,42 @@ +from SimPEG import * +from SimPEG.Utils import surface2ind_topo + + +def run(plotIt=False, nx = 5, ny = 5): + """ + Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below + a topographic surface. + + """ + + mesh = Mesh.TensorMesh([nx,ny], x0='CC') # 2D mesh + xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max()) + topo = 0.4*np.sin(xtopo*5) # define a topographic surface + + Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array + + indcc = surface2ind_topo(mesh, Topo,'CC') + + if plotIt: + from matplotlib.pylab import plt + from scipy.interpolate import interp1d + fig, ax = plt.subplots(1,1,figsize=(6,6)) + mesh.plotGrid(ax=ax, nodes=True, centers=True) + ax.plot(xtopo,topo,'k',linewidth=1) + # ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3) + ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3) + + + aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T + a = aveN2CC * indcc + a[a > 0] = 1. + a[a < 0.25] = np.nan + a = a.reshape(mesh.vnN, order='F') + masked_array = np.ma.array(a, mask=np.isnan(a)) + ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2) + plt.show() + + + +if __name__ == '__main__': + run(plotIt=True) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 18c1994f..c7597fee 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,3 +7,4 @@ from CounterUtils import * import ModelBuilder import SolverUtils from coordutils import * +from modelutils import * diff --git a/SimPEG/Utils/modelutils.py b/SimPEG/Utils/modelutils.py new file mode 100644 index 00000000..dad92fae --- /dev/null +++ b/SimPEG/Utils/modelutils.py @@ -0,0 +1,63 @@ +from matutils import mkvc, ndgrid +import numpy as np + +def surface2ind_topo(mesh, topo, gridLoc='CC'): +# def genActiveindfromTopo(mesh, topo): + """ + Get active indices from topography + """ + + + if mesh.dim == 3: + from scipy.interpolate import NearestNDInterpolator + Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2]) + + if gridLoc == 'CC': + XY = ndgrid(mesh.vectorCCx, mesh.vectorCCy) + Zcc = mesh.gridCC[:,2].reshape((np.prod(mesh.vnC[:2]), mesh.nCz), order='F') + + gridTopo = Ftopo(XY) + actind = [gridTopo[ixy] <= Zcc[ixy,:] for ixy in range(np.prod(mesh.vnC[0]))] + actind = np.hstack(actind) + + elif gridLoc == 'N': + + XY = ndgrid(mesh.vectorNx, mesh.vectorNy) + gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F') + + if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: + raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) + + Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes + actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') + + for ii in range(mesh.nCx): + for jj in range(mesh.nCy): + actind[ii,jj,:] = [np.all(gridTopo[ii:ii+2, jj:jj+2] >= Nz[kk]) for kk in range(len(Nz)) ] + + elif mesh.dim == 2: + from scipy.interpolate import interp1d + Ftopo = interp1d(topo[:,0], topo[:,1]) + + if gridLoc == 'CC': + gridTopo = Ftopo(mesh.gridCC[:,0]) + actind = mesh.gridCC[:,1] <= gridTopo + + elif gridLoc == 'N': + + gridTopo = Ftopo(mesh.vectorNx) + if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: + raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) + + Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes + actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') + + for ii in range(mesh.nCx): + actind[ii,:] = [np.all(gridTopo[ii:ii+2] > Ny[kk]) for kk in range(len(Ny)) ] + + else: + raise NotImplementedError('surface2ind_topo not implemented for 1D mesh') + + return mkvc(actind) + + From dd45a6a0855a90bc1ef382293661a5240f700463 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 2 May 2016 11:40:02 -0700 Subject: [PATCH 4/7] name updates in DC_Forward_PseudoSection, DC_Utils, example for Utils_surface2ind_topo --- SimPEG/DCIP/DCIPUtils.py | 113 ++++++++++----------- SimPEG/EM/Base.py | 3 +- SimPEG/Examples/Utils_surface2ind_topo.py | 1 - SimPEG/Examples/__init__.py | 3 +- docs/examples/DC_Forward_PseudoSection.rst | 6 +- docs/examples/Utils_surface2ind_topo.rst | 24 +++++ 6 files changed, 82 insertions(+), 68 deletions(-) create mode 100644 docs/examples/Utils_surface2ind_topo.rst diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 91eab0b3..800b478c 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -1,4 +1,4 @@ -from SimPEG import np +from SimPEG import np, Utils import BaseDC as DC import BaseDC as IP @@ -118,34 +118,27 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"): def readUBC_DC2DModel(fileName): """ - Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg + Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg - Input: - :param fileName, path to the UBC GIF 2D model file - - Output: - :param SimPEG TensorMesh 2D object - :return - - Created on Thu Nov 12 13:14:10 2015 - - @author: dominiquef + :param string fileName: path to the UBC GIF 2D model file + :rtype: TensorMesh + :return: SimPEG TensorMesh 2D object """ from SimPEG import np, mkvc # Open fileand skip header... assume that we know the mesh already - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + obsfile = np.genfromtxt(fileName, delimiter=' \n', dtype=np.str, comments='!') - dim = np.array(obsfile[0].split(),dtype=float) + dim = np.array(obsfile[0].split(), dtype=float) - temp = np.array(obsfile[1].split(),dtype=float) + temp = np.array(obsfile[1].split(), dtype=float) if len(temp) > 1: model = np.zeros(dim) for ii in range(len(obsfile)-1): - mm = np.array(obsfile[ii+1].split(),dtype=float) + mm = np.array(obsfile[ii+1].split(), dtype=float) model[:,ii] = mm model = model[:,::-1] @@ -153,10 +146,10 @@ def readUBC_DC2DModel(fileName): else: if len(obsfile[1:])==1: - mm = np.array(obsfile[1:].split(),dtype=float) + mm = np.array(obsfile[1:].split(), dtype=float) else: - mm = np.array(obsfile[1:],dtype=float) + mm = np.array(obsfile[1:], dtype=float) # Permute the second dimension to flip the order model = mm.reshape(dim[1],dim[0]) @@ -171,17 +164,16 @@ def readUBC_DC2DModel(fileName): def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None): """ - Read list of 2D tx-rx location and plot a speudo-section of apparent - resistivity. + Read list of 2D tx-rx location and plot a speudo-section of apparent + resistivity. - Assumes flat topo for now... + Assumes flat topo for now... - Input: - :param DCsurvey, z0 - :switch surveyType -> Either 'pole-dipole' | 'dipole-dipole' - :switch unitType=-> Either 'appResistivity' | 'appConductivity' | 'volt' - Output: - :figure scatter plot overlayed on image + :param SurveyDC DCsurvey: + :param string surveyType: Either 'pole-dipole' | 'dipole-dipole' + :param string unitType: Either 'appResistivity' | 'appConductivity' | 'volt' + :rtype: matplotlib.plt + :return: figure scatter plot overlayed on image """ from SimPEG import np @@ -294,27 +286,24 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt return ph -def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): +def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx): """ - Load in endpoints and survey specifications to generate Tx, Rx location - stations. + Load in endpoints and survey specifications to generate Tx, Rx location + stations. - Assumes flat topo for now... + Assumes flat topo for now... - Input: - :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] - :object mesh -> SimPEG mesh object - :switch surveyType -> 'dipole-dipole' | 'pole-dipole' | 'gradient' - : param a, n -> pole seperation, number of rx dipoles per tx + :param numpy.array endl: input endpoints [[x1, y1] , [x2, y2]] + :param Mesh mesh: SimPEG mesh object + :param string surveyType: 'dipole-dipole' | 'pole-dipole' | 'gradient' + :param float AM_sep: transmitter (A) - receiver (M) seperation + :param float b: receiver dipole seperation + :param float nrx: 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 + :rtype: DC.Survey, Src, Rx + :returns: DC survey, Source - Created on Wed December 9th, 2015 - - @author: dominiquef - !! Require clean up to deal with DCsurvey + !! Require clean up to deal with DCsurvey """ from SimPEG import np @@ -330,17 +319,17 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): 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 ) + nstn = np.floor( dl_len / AM_sep ) # 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 + stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep # Create line of P1 locations M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] # Create line of P2 locations - N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + N = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] ## Build list of Tx-Rx locations depending on survey type # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn] @@ -366,22 +355,22 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): 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]) + nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx]) # 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 + stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep # Create receiver poles # Create line of P1 locations P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] # Create line of P2 locations - P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + P2 = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] Rx.append(np.c_[P1,P2]) rxClass = DC.RxDipole(P1, P2) @@ -400,23 +389,23 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): Tx.append(np.c_[M[0,:],N[-1,:]]) # Get the edge limit of survey area - min_x = endl[0,0] + dl_x * b - min_y = endl[0,1] + dl_y * b + min_x = endl[0,0] + dl_x * MN_sep + min_y = endl[0,1] + dl_y * MN_sep - max_x = endl[1,0] - dl_x * b - max_y = endl[1,1] - dl_y * b + max_x = endl[1,0] - dl_x * MN_sep + max_y = endl[1,1] - dl_y * MN_sep 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 ) + nstn = np.floor( box_l / AM_sep ) # 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 + stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep # Define number of cross lines - nlin = int(np.floor( box_w / a )) + nlin = int(np.floor( box_w / AM_sep )) lind = range(-nlin,nlin+1) ngrad = nstn * len(lind) @@ -425,12 +414,12 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): 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 + lxx = stn_x - lind[ii]*AM_sep*dl_y + lyy = stn_y + lind[ii]*AM_sep*dl_x M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]] - N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + N = np.c_[ lxx+AM_sep*dl_x, lyy+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N] diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index a16cdb91..032b4429 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -165,7 +165,8 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRho` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv + # (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @property diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py index 61a0f937..4ed1cdaf 100644 --- a/SimPEG/Examples/Utils_surface2ind_topo.py +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -37,6 +37,5 @@ def run(plotIt=False, nx = 5, ny = 5): plt.show() - if __name__ == '__main__': run(plotIt=True) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index c67652a2..24828466 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -20,8 +20,9 @@ import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation import MT_1D_ForwardAndInversion import MT_3D_Foward +import Utils_surface2ind_topo -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"] ##### AUTOIMPORTS ##### diff --git a/docs/examples/DC_Forward_PseudoSection.rst b/docs/examples/DC_Forward_PseudoSection.rst index 80cf0307..4231e944 100644 --- a/docs/examples/DC_Forward_PseudoSection.rst +++ b/docs/examples/DC_Forward_PseudoSection.rst @@ -20,9 +20,9 @@ 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 +surveyType = survey type 'pole-dipole' or 'dipole-dipole' +unitType = Data type "appResistivity" | "appConductivity" | "volt" +Created by @fourndo diff --git a/docs/examples/Utils_surface2ind_topo.rst b/docs/examples/Utils_surface2ind_topo.rst new file mode 100644 index 00000000..04b27d5d --- /dev/null +++ b/docs/examples/Utils_surface2ind_topo.rst @@ -0,0 +1,24 @@ +.. _examples_Utils_surface2ind_topo: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below +a topographic surface. + + + +.. plot:: + + from SimPEG import Examples + Examples.Utils_surface2ind_topo.run() + +.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py + :language: python + :linenos: From dbdcc3cefb28cfaa7802b28468d5f0cc60f70e9f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 22:06:32 -0700 Subject: [PATCH 5/7] use sigma in MfRhoDeriv - due to propmap bug --- SimPEG/EM/Base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 032b4429..a16cdb91 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -165,8 +165,7 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRho` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv - # (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @property From 66440b0478b434adaa9bbb07c13ce056fd3f45f3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 22:14:41 -0700 Subject: [PATCH 6/7] add depreciation warnings to DCIP utils for activeind from topo --- SimPEG/DCIP/DCIPUtils.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 800b478c..f4e643eb 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -1,12 +1,16 @@ from SimPEG import np, Utils import BaseDC as DC import BaseDC as IP +import warnings def getActiveindfromTopo(mesh, topo): # def genActiveindfromTopo(mesh, topo): """ Get active indices from topography """ + warnings.warn( + "`getActiveindfromTopo` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead", + FutureWarning) from scipy.interpolate import NearestNDInterpolator if mesh.dim==3: nCxy = mesh.nCx*mesh.nCy @@ -28,6 +32,9 @@ def gettopoCC(mesh, airind): """ Get topography from active indices of mesh. """ + warnings.warn( + "`gettopoCC` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead", + FutureWarning) mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2]) zc = mesh.gridCC[:,2] AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F') @@ -630,16 +637,9 @@ def readUBC_DC3Dobs(fileName): """ Read UBC GIF DCIP 3D observation file and generate survey - Input: - :param fileName, path to the UBC GIF 3D obs file - - Output: - :param DCIPsurvey - :return - - Created on Mon April 6th, 2015 - - @author: dominiquef + :param string fileName:, path to the UBC GIF 3D obs file + :rtype: Survey + :return: DCIPsurvey """ From 0379df2bf22bed20151ed29ece5e41e95e6e98cf Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 22:27:02 -0700 Subject: [PATCH 7/7] attempt to clean up docs in DCIP utils --- SimPEG/DCIP/DCIPUtils.py | 73 ++++++++++++---------------------------- 1 file changed, 22 insertions(+), 51 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index f4e643eb..49f9ba2b 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -444,21 +444,14 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): """ 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 dim -> either '2D' | '3D' - :string surveyType -> either 'SURFACE' | 'GENERAL' - - Output: - :param UBC2D-Data file - :return - - Last edit: February 16th, 2016 - - @author: dominiquef - + :param string fileName: including path where the file is written out + :param Survey DCsurvey: DC survey class object + :param string dim: either '2D' | '3D' + :param string surveyType: either 'SURFACE' | 'GENERAL' + :rtype: file + :return: UBC2D-Data file """ + from SimPEG import mkvc assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'" @@ -539,15 +532,9 @@ def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'): The Z value is preserved, but Y coordinates zeroed. - Input: - :param survey3D - - Output: - :figure survey2D - - Edited April 6th, 2016 - - @author: dominiquef + :param DC.Survey survey3D: 3D simpeg DC survey + :rtype: DC.Survey + :return: survey2D """ from SimPEG import np @@ -715,17 +702,9 @@ def readUBC_DC2Dobs(fileName): ------- NEEDS TO BE UPDATED ------ Read UBC GIF 2D observation file and generate arrays for tx-rx location - Input: - :param fileName, path to the UBC GIF 2D model file - - Output: - :param rx, tx - :return - - Created on Thu Nov 12 13:14:10 2015 - - @author: dominiquef - + :param string fileName: path to the UBC GIF 2D model file + :rtype: (DC.Src, DC.Rx, ??, ??) + :return: source_locs, rx_locs, ??, ?? """ from SimPEG import np @@ -765,11 +744,9 @@ def readUBC_DC2Dpre(fileName): 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 - - Output: - DCsurvey - :return + :param string fileName: path to the UBC GIF 3D obs file + :rtype: DC.Survey + :return: DCsurvey Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >> @@ -831,12 +808,9 @@ def readUBC_DC2DMesh(fileName): """ Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg - Input: - :param fileName, path to the UBC GIF mesh file - - Output: - :param SimPEG TensorMesh 2D object - :return + :param string fileName: path to the UBC GIF mesh file + :rtype: Mesh.TensorMesh + :return: SimPEG TensorMesh 2D object Created on Thu Nov 12 13:14:10 2015 @@ -902,12 +876,9 @@ def xy_2_lineID(DCsurvey): they were collected. May need to generalize for random point locations, but will be more expensive - Input: - :param DCdict Vectors of station location - - Output: - :param LineID Vector of integers - :return + :param numpy.array DCdict: Vectors of station location + :rtype: numpy.array + :return: LineID Vector of integers Created on Thu Feb 11, 2015