From c9017bd1801df7fe9e43f80f1f2e42e75146854f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 31 May 2016 08:42:52 -0700 Subject: [PATCH] Delete DCIP folder. --- SimPEG/DCIP/BaseDC.py | 292 ----------- SimPEG/DCIP/BaseIP.py | 182 ------- SimPEG/DCIP/DCIPUtils.py | 1043 -------------------------------------- SimPEG/DCIP/Utils.py | 38 -- SimPEG/DCIP/__init__.py | 4 - 5 files changed, 1559 deletions(-) delete mode 100644 SimPEG/DCIP/BaseDC.py delete mode 100644 SimPEG/DCIP/BaseIP.py delete mode 100644 SimPEG/DCIP/DCIPUtils.py delete mode 100644 SimPEG/DCIP/Utils.py delete mode 100644 SimPEG/DCIP/__init__.py diff --git a/SimPEG/DCIP/BaseDC.py b/SimPEG/DCIP/BaseDC.py deleted file mode 100644 index 475c621e..00000000 --- a/SimPEG/DCIP/BaseDC.py +++ /dev/null @@ -1,292 +0,0 @@ -from SimPEG import * - -class FieldsDC_CC(Problem.Fields): - knownFields = {'phi_sol':'CC'} - aliasFields = { - 'phi' : ['phi_sol','CC','_phi'], - 'e' : ['phi_sol','F','_e'], - 'j' : ['phi_sol','F','_j'] - } - - def __init__(self,mesh,survey,**kwargs): - super(FieldsDC_CC, self).__init__(mesh, survey, **kwargs) - - def startup(self): - self._cellGrad = self.survey.prob.mesh.cellGrad - self._Mfinv = self.survey.prob.mesh.getFaceInnerProduct(invMat=True) - - def _phi(self, phi_sol, srcList): - phi = phi_sol - # for i, src in enumerate(srcList): - # phi_p = src.phi_p(self.survey.prob) - # if phi_p is not None: - # phi[:,i] += phi_p - return phi - - def _e(self, phi_sol, srcList): - e = -self._cellGrad*phi_sol - # for i, src in enumerate(srcList): - # e_p = src.e_p(self.survey.prob) - # if e_p is not None: - # e[:,i] += e_p - return e - - def _j(self, phi_sol, srcList): - - j = -self._Mfinv*self.survey.prob.Msig*self._cellGrad*phi_sol - # for i, src in enumerate(srcList): - # j_p = src.j_p(self.survey.prob) - # if j_p is not None: - # j[:,i] += j_p - return j - - - -class SrcDipole(Survey.BaseSrc): - """A dipole source, locA and locB are moved to the closest cell-centers""" - - current = 1 - loc = None - # _rhsDict = None - - def __init__(self, rxList, locA, locB, **kwargs): - self.loc = (locA, locB) - super(SrcDipole, self).__init__(rxList, **kwargs) - - def eval(self, prob): - # Recompute rhs - # if getattr(self, '_rhsDict', None) is None: - # self._rhsDict = {} - # if mesh not in self._rhsDict: - pts = [self.loc[0], self.loc[1]] - inds = Utils.closestPoints(prob.mesh, pts) - q = np.zeros(prob.mesh.nC) - q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] ) - # self._rhsDict[mesh] = q - # return self._rhsDict[mesh] - return q - - -class RxDipole(Survey.BaseRx): - """A dipole source, locA and locB are moved to the closest cell-centers""" - def __init__(self, locsM, locsN, **kwargs): - locs = (locsM, locsN) - assert locsM.shape == locsN.shape, 'locs must be the same shape.' - super(RxDipole, self).__init__(locs, 'dipole', storeProjections=False, **kwargs) - - @property - def nD(self): - """Number of data in the receiver.""" - return self.locs[0].shape[0] - - def getP(self, mesh): - P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc) - P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc) - return P0 - P1 - - -class SurveyDC(Survey.BaseSurvey): - """ - **SurveyDC** - - Geophysical DC resistivity data. - - """ - uncert = None - def __init__(self, srcList, **kwargs): - self.srcList = srcList - Survey.BaseSurvey.__init__(self, **kwargs) - # self._rhsDict = {} - self._Ps = {} - - def eval(self, u): - """ - Predicted data. - - .. math:: - d_\\text{pred} = Pu(m) - """ - P = self.getP(self.prob.mesh) - return P*mkvc(u[self.srcList, 'phi_sol']) - - def getP(self, mesh): - if mesh in self._Ps: - return self._Ps[mesh] - - P_src = [sp.vstack([rx.getP(mesh) for rx in src.rxList]) for src in self.srcList] - - self._Ps[mesh] = sp.block_diag(P_src) - return self._Ps[mesh] - - -class ProblemDC_CC(Problem.BaseProblem): - """ - **ProblemDC** - - Geophysical DC resistivity problem. - - """ - - surveyPair = SurveyDC - Solver = Solver - fieldsPair = FieldsDC_CC - Ainv = None - - def __init__(self, mesh, **kwargs): - Problem.BaseProblem.__init__(self, mesh) - self.mesh.setCellGradBC('neumann') - Utils.setKwargs(self, **kwargs) - - - deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig'] - - @property - def Msig(self): - if getattr(self, '_Msig', None) is None: - sigma = self.curModel.transform - Av = self.mesh.aveF2CC - self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma))) - return self._Msig - - @property - def dMdsig(self): - if getattr(self, '_dMdsig', None) is None: - sigma = self.curModel.transform - Av = self.mesh.aveF2CC - dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2) - self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop - return self._dMdsig - - @property - def A(self): - """ - Makes the matrix A(m) for the DC resistivity problem. - - :param numpy.array m: model - :rtype: scipy.csc_matrix - :return: A(m) - - .. math:: - c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0 - - Where M() is the mass matrix and mT is the model transform. - """ - if getattr(self, '_A', None) is None: - D = self.mesh.faceDiv - G = self.mesh.cellGrad - self._A = D*self.Msig*G - # Remove the null space from the matrix. - self._A[0,0] /= self.mesh.vol[0] - self._A = self._A.tocsc() - return self._A - - def getRHS(self): - # if self.mesh not in self._rhsDict: - RHS = np.array([src.eval(self) for src in self.survey.srcList]).T - # self._rhsDict[mesh] = RHS - # return self._rhsDict[mesh] - return RHS - - def fields(self, m): - - F = self.fieldsPair(self.mesh, self.survey) - self.curModel = m - A = self.A - self.Ainv = self.Solver(A, **self.solverOpts) - RHS = self.getRHS() - Phi = self.Ainv * RHS - Srcs = self.survey.srcList - F[Srcs, 'phi_sol'] = Phi - - return F - - def Jvec(self, m, v, f=None): - """ - :param numpy.array m: model - :param numpy.array v: vector to multiply - :param Fields f: fields - :rtype: numpy.array - :return: Jv - - .. math:: - c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0 - - \\nabla_u (A(m)u - q) = A(m) - - \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) - - Where M() is the mass matrix and mT is the model transform. - - .. math:: - J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u) - - J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) ) - """ - # Set current model; clear dependent property $\mathbf{A(m)}$ - self.curModel = m - sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - if f is None: - # Run forward simulation if $u$ not provided - f = self.fields(self.curModel) - u = f[self.survey.srcList, 'phi_sol'] - - D = self.mesh.faceDiv - G = self.mesh.cellGrad - # Derivative of model transform, $\deriv{\sigma}{\m}$ - dsigdm_x_v = self.curModel.transformDeriv * v - - # Take derivative of $C(m,u)$ w.r.t. $m$ - dCdm_x_v = np.empty_like(u) - # loop over fields for each source - for i in range(self.survey.nSrc): - # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$ - dAdsig = D * self.dMdsig( G * u[:,i] ) - dCdm_x_v[:, i] = dAdsig * dsigdm_x_v - - # Take derivative of $C(m,u)$ w.r.t. $u$ - dA_du = self.A - # Solve for $\deriv{u}{m}$ - # dCdu_inv = self.Solver(dCdu, **self.solverOpts) - if self.Ainv is None: - self.Ainv = self.Solver(dA_du, **self.solverOpts) - - P = self.survey.getP(self.mesh) - Jv = - P * mkvc( self.Ainv * dCdm_x_v ) - return Jv - - def Jtvec(self, m, v, f=None): - - self.curModel = m - sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - if f is None: - # Run forward simulation if $f$ not provided - f = self.fields(self.curModel) - u = f[self.survey.srcList, 'phi_sol'] - - shp = u.shape - P = self.survey.getP(self.mesh) - PT_x_v = (P.T*v).reshape(shp, order='F') - - D = self.mesh.faceDiv - G = self.mesh.cellGrad - dA_du = self.A - mT_dm = self.mapping.deriv(m) - - # We probably always need this due to the linesearch .. (?) - self.Ainv = self.Solver(dA_du.T, **self.solverOpts) - # if self.Ainv is None: - # self.Ainv = self.Solver(dCdu, **self.solverOpts) - - w = self.Ainv * PT_x_v - - Jtv = 0 - for i, ui in enumerate(u.T): # loop over each column - Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] ) - - Jtv = - mT_dm.T * ( Jtv ) - return Jtv - - - - - diff --git a/SimPEG/DCIP/BaseIP.py b/SimPEG/DCIP/BaseIP.py deleted file mode 100644 index a18b5a47..00000000 --- a/SimPEG/DCIP/BaseIP.py +++ /dev/null @@ -1,182 +0,0 @@ -from SimPEG import * -from BaseDC import SurveyDC, FieldsDC_CC - -class SurveyIP(SurveyDC): - """ - **SurveyDC** - - Geophysical DC resistivity data. - - """ - - def __init__(self, srcList, **kwargs): - self.srcList = srcList - Survey.BaseSurvey.__init__(self, **kwargs) - self._Ps = {} - - def dpred(self, m, f=None): - """ - Predicted data. - - .. math:: - d_\\text{pred} = Pf(m) - """ - - return self.prob.forward(m) - - -class ProblemIP(Problem.BaseProblem): - """ - **ProblemIP** - - Geophysical IP resistivity problem. - - """ - - surveyPair = SurveyDC - Solver = Solver - sigma = None - Ainv = None - u = None - - def __init__(self, mesh, **kwargs): - Problem.BaseProblem.__init__(self, mesh) - self.mesh.setCellGradBC('neumann') - Utils.setKwargs(self, **kwargs) - - # deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig'] - - @property - def Msig(self): - if getattr(self, '_Msig', None) is None: - # sigma = self.curModel.transform - sigma = self.sigma - Av = self.mesh.aveF2CC - self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma))) - return self._Msig - - @property - def dMdsig(self): - if getattr(self, '_dMdsig', None) is None: - # sigma = self.curModel.transform - sigma = self.sigma - Av = self.mesh.aveF2CC - dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2) - self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop - return self._dMdsig - - @property - def A(self): - """ - Makes the matrix A(m) for the DC resistivity problem. - - :param numpy.array m: model - :rtype: scipy.csc_matrix - :return: A(m) - - .. math:: - c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0 - - Where M() is the mass matrix and mT is the model transform. - """ - if getattr(self, '_A', None) is None: - D = self.mesh.faceDiv - G = self.mesh.cellGrad - self._A = D*self.Msig*G - # Remove the null space from the matrix. - self._A[-1,-1] /= self.mesh.vol[-1] - self._A = self._A.tocsc() - return self._A - - def getRHS(self): - # if self.mesh not in self._rhsDict: - RHS = np.array([src.eval(self) for src in self.survey.srcList]).T - # self._rhsDict[mesh] = RHS - # return self._rhsDict[mesh] - return RHS - - def fields(self, m): - if self.u is None: - A = self.A - if self.Ainv == None: - self.Ainv = self.Solver(A, **self.solverOpts) - Q = self.getRHS() - self.u = self.Ainv * Q - return self.u - - def forward(self, m, u=None): - # Set current model; clear dependent property $\mathbf{A(m)}$ - self.curModel = m - # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - sigma = self.sigma - if self.u is None: - # Run forward simulation if $u$ not provided - u = self.fields(sigma) - - shp = (self.mesh.nC, self.survey.nSrc) - u = self.u.reshape(shp, order='F') - - D = self.mesh.faceDiv - G = self.mesh.cellGrad - # Derivative of model transform, $\deriv{\sigma}{\m}$ - # dsigdm_x_v = self.curModel.transformDeriv * v - - dsigdm_x_v = Utils.sdiag(sigma) * self.curModel.transformDeriv * m - - # Take derivative of $C(m,u)$ w.r.t. $m$ - dCdm_x_v = np.empty_like(u) - # loop over fields for each source - for i in range(self.survey.nSrc): - # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$ - dAdsig = D * self.dMdsig( G * u[:,i] ) - dCdm_x_v[:, i] = dAdsig * dsigdm_x_v - - # Take derivative of $C(m,u)$ w.r.t. $u$ - - if self.Ainv == None: - self.Ainv = self.Solver(A, **self.solverOpts) - - # dCdu = self.A - # Solve for $\deriv{u}{m}$ - # dCdu_inv = self.Solver(dCdu, **self.solverOpts) - P = self.survey.getP(self.mesh) - J_x_v = - P * mkvc( self.Ainv * dCdm_x_v ) - return -J_x_v - - def Jvec(self, m, v, f=None): - return self.forward(v) - - def Jtvec(self, m, v, f=None): - - self.curModel = m - # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - sigma = self.sigma - if self.u is None: - u = self.fields(sigma) - else: - u = self.u - shp = (self.mesh.nC, self.survey.nSrc) - u = u.reshape(shp, order='F') - P = self.survey.getP(self.mesh) - PT_x_v = (P.T*v).reshape(shp, order='F') - - D = self.mesh.faceDiv - G = self.mesh.cellGrad - A = self.A - mT_dm = Utils.sdiag(sigma)*self.mapping.deriv(m) - # mT_dm = self.mapping.deriv(m) - - # dCdu = A.T - # Ainv = self.Solver(dCdu, **self.solverOpts) - # if self.Ainv == None: - self.Ainv = self.Solver(A.T, **self.solverOpts) - - w = self.Ainv * PT_x_v - - Jtv = 0 - for i, ui in enumerate(u.T): # loop over each column - Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] ) - - Jtv = - mT_dm.T * ( Jtv ) - return -Jtv - diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py deleted file mode 100644 index 8e1ca5e9..00000000 --- a/SimPEG/DCIP/DCIPUtils.py +++ /dev/null @@ -1,1043 +0,0 @@ -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 - Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F') - Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2]) - XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy) - XY.shape - topo = Ftopo(XY) - actind = [] - for ixy in range(nCxy): - actind.append(topo[ixy] <= Zcc[ixy,:]) - else: - raise NotImplementedError("Only 3D is working") - - return Utils.mkvc(np.vstack(actind)) - -def gettopoCC(mesh, airind): -# 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') - ZC = zc.reshape((mesh.vnC[0]*mesh.vnC[1], mesh.vnC[2]), order='F') - topo = np.zeros(ZC.shape[0]) - topoCC = np.zeros(ZC.shape[0]) - for i in range(ZC.shape[0]): - ind = np.argmax(ZC[i,:][~AIRIND[i,:]]) - topo[i] = ZC[i,:][~AIRIND[i,:]].max() + mesh.hz[~AIRIND[i,:]][ind]*0.5 - topoCC[i] = ZC[i,:][~AIRIND[i,:]].max() - XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy) - return mesh2D, topoCC - -def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"): - """ - Seogi's personal readObs function. - - """ - text_file = open(filename, "r") - lines = text_file.readlines() - text_file.close() - SRC = [] - DATA = [] - srcLists = [] - isrc = 0 - # airind = getActiveindfromTopo(mesh, topo) - # mesh2D, topoCC = gettopoCC(mesh, airind) - - for line in lines: - if "!" in line.split(): continue - elif line == '\n': continue - elif line == ' \n': continue - temp = map(float, line.split()) - # Read a line for the current electrode - if len(temp) == 5: # SRC: Only X and Y are provided (assume no topography) - #TODO consider topography and assign the closest cell center in the earth - if isrc == 0: - DATA_temp = [] - else: - DATA.append(np.asarray(DATA_temp)) - DATA_temp = [] - indM = Utils.closestPoints(mesh2D, DATA[isrc-1][:,1:3]) - indN = Utils.closestPoints(mesh2D, DATA[isrc-1][:,3:5]) - rx = DCIP.RxDipole(np.c_[DATA[isrc-1][:,1:3], topoCC[indM]], np.c_[DATA[isrc-1][:,3:5], topoCC[indN]]) - temp = np.asarray(temp) - if [SRC[isrc-1][0], SRC[isrc-1][1]] == [SRC[isrc-1][2], SRC[isrc-1][3]]: - indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]]) - tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[mesh.vectorCCx.max(), mesh.vectorCCy.max(), topoCC[-1]]) - else: - indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]]) - indB = Utils.closestPoints(mesh2D, [SRC[isrc-1][2], SRC[isrc-1][3]]) - tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[SRC[isrc-1][2], SRC[isrc-1][3], topoCC[indB]]) - srcLists.append(tx) - SRC.append(temp) - isrc += 1 - elif len(temp) == 7: # SRC: X, Y and Z are provided - SRC.append(temp) - isrc += 1 - elif len(temp) == 6: # - DATA_temp.append(np.r_[isrc, np.asarray(temp)]) - elif len(temp) > 7: - DATA_temp.append(np.r_[isrc, np.asarray(temp)]) - - DATA.append(np.asarray(DATA_temp)) - DATA_temp = [] - indM = Utils.closestPoints(mesh2D, DATA[isrc-1][:,1:3]) - indN = Utils.closestPoints(mesh2D, DATA[isrc-1][:,3:5]) - rx = DCIP.RxDipole(np.c_[DATA[isrc-1][:,1:3], topoCC[indM]], np.c_[DATA[isrc-1][:,3:5], topoCC[indN]]) - temp = np.asarray(temp) - if [SRC[isrc-1][0], SRC[isrc-1][1]] == [SRC[isrc-1][2], SRC[isrc-1][3]]: - indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]]) - tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[mesh.vectorCCx.max(), mesh.vectorCCy.max(), topoCC[-1]]) - else: - indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]]) - indB = Utils.closestPoints(mesh2D, [SRC[isrc-1][2], SRC[isrc-1][3]]) - tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[SRC[isrc-1][2], SRC[isrc-1][3], topoCC[indB]]) - srcLists.append(tx) - text_file.close() - survey = DCIP.SurveyDC(srcLists) - - # Do we need this? - SRC = np.asarray(SRC) - DATA = np.vstack(DATA) - survey.dobs = np.vstack(DATA)[:,-2] - - - return {'DCsurvey':survey, 'airind':airind, 'topoCC':topoCC, 'SRC':SRC} - -def readUBC_DC2DModel(fileName): - """ - Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg - - :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='!') - - dim = np.array(obsfile[0].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) - model[:,ii] = mm - - model = model[:,::-1] - - else: - - if len(obsfile[1:])==1: - mm = np.array(obsfile[1:].split(), dtype=float) - - else: - mm = np.array(obsfile[1:], dtype=float) - - # Permute the second dimension to flip the order - model = mm.reshape(dim[1],dim[0]) - - model = model[::-1,:] - model = np.transpose(model, (1, 0)) - - model = mkvc(model) - - - return model - - -def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None): - """ - Read list of 2D tx-rx location and plot a speudo-section of apparent - resistivity. - - Assumes flat topo for now... - - :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 - from scipy.interpolate import griddata - import pylab as plt - - # Pre-allocate - midx = [] - midz = [] - rho = [] - 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 - MA = np.abs(Tx[0][0] - Rx[0][:,0]) - MB = np.abs(Tx[1][0] - Rx[0][:,0]) - NB = np.abs(Tx[1][0] - Rx[1][:,0]) - NA = np.abs(Tx[0][0] - Rx[1][:,0]) - MN = np.abs(Rx[1][:,0] - Rx[0][:,0]) - - # Create mid-point location - Cmid = (Tx[0][0] + Tx[1][0])/2 - Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 - - # Change output for unitType - if unitType == 'volt': - - rho = np.hstack([rho,data]) - - else: - - # Compute pant leg of apparent rho - if surveyType == 'pole-dipole': - - leg = data * 2*np.pi * MA * ( MA + MN ) / MN - - elif surveyType == 'dipole-dipole': - - leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA ) - - else: - print """unitType must be 'pole-dipole' | 'dipole-dipole' """ - break - - - if unitType == 'appConductivity': - - leg = np.log10(abs(1./leg)) - rho = np.hstack([rho,leg]) - - elif unitType == 'appResistivity': - - leg = np.log10(abs(leg)) - rho = np.hstack([rho,leg]) - - else: - print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """ - break - - midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) - midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ]) - - # 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') - - # Scale the color scheme - if clim == None: - vmin, vmax = rho.min(), rho.max() - else: - vmin, vmax = clim[0], clim[1] - - # Plot data - grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho) - - - - ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax) - plt.gca().tick_params(axis='both', which='major', labelsize=8) - - if contour is not None: - plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax) - - # Add scatter points - axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax) - - if colorbar: - - if unitType == 'volt': - cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal") - - else: - cbar = plt.colorbar(ph, ax = axs, 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 unitType == 'appConductivity': - cbar.set_label("App.Cond",size=12) - elif unitType == 'appResistivity': - cbar.set_label("App.Res.",size=12) - elif unitType == 'volt': - cbar.set_label("Potential (V)",size=12) - - - if not axlabel: - axs.set_xticklabels([]) - axs.set_yticklabels([]) - - plt.gca().set_aspect('equal', adjustable='box') - - - - return ph - -def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx): - """ - Load in endpoints and survey specifications to generate Tx, Rx location - stations. - - Assumes flat topo for now... - - :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 - - :rtype: DC.Survey, Src, Rx - :returns: DC survey, Source - - !! 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 / AM_sep ) - - # Compute discrete pole location along line - 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+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] - # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B] - Tx = [] - Rx = [] - SrcList = [] - - - if surveyType != 'gradient': - - for ii in range(0, int(nstn)-1): - - - if surveyType == 'dipole-dipole': - tx = np.c_[M[ii,:],N[ii,:]] - elif surveyType == 'pole-dipole': - 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 - 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*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+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) - Tx.append(tx) - if surveyType == 'dipole-dipole': - srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:]) - elif surveyType == 'pole-dipole': - srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:]) - SrcList.append(srcClass) - - 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 - - Tx.append(np.c_[M[0,:],N[-1,:]]) - - # Get the edge limit of survey area - 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 * 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 / AM_sep ) - - # Compute discrete pole location along line - 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 / AM_sep )) - 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]*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+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] - - Rx.append(rx) - rxClass = DC.RxDipole(rx[:,:3], rx[:,3:]) - srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:]) - SrcList.append(srcClass) - else: - print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """ - - survey = DC.SurveyDC(SrcList) - return survey, Tx, Rx - - -def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0): - """ - Write UBC GIF DCIP 2D or 3D observation file - - :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'" - assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" - - fid = open(fileName,'w') - fid.write('! ' + surveyType + ' FORMAT\n') - - 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 dim and surveyType - if dim=='2D': - - if surveyType == '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 surveyType == 'SURFACE': - - fid.writelines("%f " % ii for ii in mkvc(tx[0,:])) - M = M[:,0] - N = N[:,0] - - if surveyType == '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 dim=='3D': - - if surveyType == 'SURFACE': - - fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:])) - M = M[:,0:2] - N = N[:,0:2] - - if surveyType == '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() - -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) - In the 'local' system, station coordinates are referenced - to distance from the first srcLoc[0].loc[0] - - The Z value is preserved, but Y coordinates zeroed. - - :param DC.Survey survey3D: 3D simpeg DC survey - :rtype: DC.Survey - :return: survey2D - - """ - from SimPEG import np - - def stn_id(v0,v1,r): - """ - Compute station ID along line - """ - - dl = int(v0.dot(v1)) * r - - return dl - - srcLists = [] - - srcMat = getSrc_locs(DCsurvey) - - # Find all unique line id - uniqueID = np.unique(lineID) - - for jj in range(len(uniqueID)): - - indx = np.where(lineID==uniqueID[jj])[0] - - # Find origin of survey - r = 1e+8 # Initialize to some large number - - Tx = srcMat[indx] - - x0 = Tx[0][0,0:2] # Define station zero along line - - vecTx, r1 = r_unit(x0,Tx[-1][1,0:2]) - - for ii in range(len(indx)): - - # Get all receivers - Rx = DCsurvey.srcList[indx[ii]].rxList[0].locs - nrx = Rx[0].shape[0] - - if 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) - elif flag == 'Yloc': - """ Flip the XY axis locs""" - A = Tx[ii][0,1] - 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]]) ) ) - - - DCsurvey2D = DC.SurveyDC(srcLists) - - DCsurvey2D.dobs = np.asarray(DCsurvey.dobs) - DCsurvey2D.std = np.asarray(DCsurvey.std) - - return DCsurvey2D - -def readUBC_DC3Dobs(fileName, rtype = 'DC'): - """ - Read UBC GIF IP 3D observation file and generate survey - - :param string fileName:, path to the UBC GIF 3D obs file - :rtype: Survey - :return: DCIPsurvey - - """ - zflag = True # Flag for z value provided - - # Load file - if rtype == 'IP': - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE') - - elif rtype == 'DC': - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - - else: - print "rtype must be 'DC'(default) | 'IP'" - - # Pre-allocate - srcLists = [] - Rx = [] - d = [] - wd = [] - - - # Countdown for number of obs/tx - count = 0 - for ii in range(obsfile.shape[0]): - - # Skip if blank line - if not obsfile[ii]: - continue - - # First line or end of a transmitter block, read transmitter info - if count==0: - # Read the line - temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T) - count = int(temp[-1]) - - # Check if z value is provided, if False -> nan - if len(temp)==5: - tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan] - - zflag = False # Pass on the flag to the receiver loc - - else: - tx = temp[:-1] - - rx = [] - continue - - temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string - - # Filter out negative IP -# if temp[-2] < 0: -# count = count -1 -# print "Negative!" -# -# else: - - # If the Z-location is provided, otherwise put nan - if zflag: - - rx.append(temp[:-2]) - # Check if there is data with the location - if len(temp)==8: - d.append(temp[-2]) - wd.append(temp[-1]) - - else: - rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] ) - # Check if there is data with the location - if len(temp)==6: - d.append(temp[-2]) - wd.append(temp[-1]) - - count = count -1 - - # Reach the end of transmitter block, append the src, rx and continue - if count == 0: - rx = np.asarray(rx) - Rx = DC.RxDipole(rx[:,:3],rx[:,3:]) - srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) ) - - # Create survey class - survey = DC.SurveyDC(srcLists) - - survey.dobs = np.asarray(d) - survey.std = np.asarray(wd) - - return {'DCsurvey':survey} - -def readUBC_DC2Dobs(fileName): - """ - ------- NEEDS TO BE UPDATED ------ - Read UBC GIF 2D observation file and generate arrays for tx-rx location - - :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 - - # Load file - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - - # Check first line and figure out if 2D or 3D file format - line = np.array(obsfile[0].split(),dtype=float) - - tx_A = [] - tx_B = [] - rx_M = [] - rx_N = [] - d = [] - wd = [] - - for ii in range(obsfile.shape[0]): - - # If len==3, then simple format where tx-rx is listed on each line - if len(line) == 4: - - temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') - tx_A = np.hstack((tx_A,temp[0])) - tx_B = np.hstack((tx_B,temp[1])) - rx_M = np.hstack((rx_M,temp[2])) - rx_N = np.hstack((rx_N,temp[3])) - - - rx = np.transpose(np.array((rx_M,rx_N))) - tx = np.transpose(np.array((tx_A,tx_B))) - - return tx, rx, d, wd - -def readUBC_DC2Dpre(fileName): - """ - Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location - - Input: - :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 !! >> - - @author: dominiquef - - """ - - # Load file - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - - # Pre-allocate - srcLists = [] - Rx = [] - d = [] - zflag = True # Flag for z value provided - - for ii in range(obsfile.shape[0]): - - if not obsfile[ii]: - continue - - # First line is transmitter with number of receivers - - - temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T) - - - # Check if z value is provided, if False -> nan - if len(temp)==5: - tx = np.r_[temp[0],np.nan,np.nan,temp[1],np.nan,np.nan] - zflag = False - - 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]) - - - Rx = DC.RxDipole(rx[:,:3],rx[:,3:]) - srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) ) - - # Create survey class - survey = DC.SurveyDC(srcLists) - - survey.dobs = np.asarray(d) - - return {'DCsurvey':survey} - -def readUBC_DC2DMesh(fileName): - """ - Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg - - :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 - - @author: dominiquef - - """ - - from SimPEG import np - # Open file - fopen = open(fileName,'r') - - # Read down the file and unpack dx vector - def unpackdx(fid,nrows): - for ii in range(nrows): - - line = fid.readline() - var = np.array(line.split(),dtype=float) - - if ii==0: - x0= var[0] - xvec = np.ones(int(var[2])) * (var[1] - var[0]) / int(var[2]) - xend = var[1] - - else: - xvec = np.hstack((xvec,np.ones(int(var[1])) * (var[0] - xend) / int(var[1]))) - xend = var[0] - - return x0, xvec - - #%% Start with dx block - # First line specifies the number of rows for x-cells - line = fopen.readline() - nl = np.array(line.split(),dtype=float) - - [x0, dx] = unpackdx(fopen,nl) - - - #%% Move down the file until reaching the z-block - line = fopen.readline() - if not line: - line = fopen.readline() - - #%% End with dz block - # First line specifies the number of rows for z-cells - line = fopen.readline() - nl = np.array(line.split(),dtype=float) - - [z0, dz] = unpackdx(fopen,nl) - - # Flip z0 to be the bottom of the mesh for SimPEG - z0 = z0 - sum(dz) - dz = dz[::-1] - #%% Make the mesh using SimPEG - - from SimPEG import Mesh - tensMsh = Mesh.TensorMesh([dx,dz],(x0, z0)) - return tensMsh - -def xy_2_lineID(DCsurvey): - """ - Read DC survey class and append line ID. - Assumes that the locations are listed in the order - they were collected. May need to generalize for random - point locations, but will be more expensive - - :param numpy.array DCdict: Vectors of station location - :rtype: numpy.array - :return: LineID Vector of integers - - Created on Thu Feb 11, 2015 - - @author: dominiquef - - """ - - # Compute unit vector between two points - nstn = DCsurvey.nSrc - - # Pre-allocate space - lineID = np.zeros(nstn) - - linenum = 0 - indx = 0 - - for ii in range(nstn): - - if ii == 0: - - A = DCsurvey.srcList[ii].loc[0] - B = DCsurvey.srcList[ii].loc[1] - - xout = np.mean([A[0:2],B[0:2]], axis = 0) - - xy0 = A[:2] - xym = xout - - # Deal with replicate pole location - if np.all(xy0==xym): - - xym[0] = xym[0] + 1e-3 - - continue - - A = DCsurvey.srcList[ii].loc[0] - B = DCsurvey.srcList[ii].loc[1] - - xin = np.mean([A[0:2],B[0:2]], axis = 0) - - # Compute vector between neighbours - vec1, r1 = r_unit(xout,xin) - - # Compute vector between current stn and mid-point - vec2, r2 = r_unit(xym,xin) - - # Compute vector between current stn and start line - vec3, r3 = r_unit(xy0,xin) - - # Compute vector between mid-point and start line - vec4, r4 = r_unit(xym,xy0) - - # Compute dot product - ang1 = np.abs(vec1.dot(vec2)) - ang2 = np.abs(vec3.dot(vec4)) - - # If the angles are smaller then 45d, than next point is on a new line - if ((ang1 < np.cos(np.pi/4.)) | (ang2 < np.cos(np.pi/4.))) & (np.all(np.r_[r1,r2,r3,r4] > 0)): - - # Re-initiate start and mid-point location - xy0 = A[:2] - xym = xin - - # Deal with replicate pole location - if np.all(xy0==xym): - - xym[0] = xym[0] + 1e-3 - - linenum += 1 - indx = ii - - else: - xym = np.mean([xy0,xin], axis = 0) - - lineID[ii] = linenum - xout = xin - - return lineID - -def r_unit(p1,p2): - """ - r_unit(x,y) : Function computes the unit vector - between two points with coordinates p1(x1,y1) and p2(x2,y2) - - """ - - assert len(p1)==len(p2), 'locs must be the same shape.' - - dx = [] - for ii in range(len(p1)): - dx.append((p2[ii] - p1[ii])) - - # Compute length of vector - r = np.linalg.norm(np.asarray(dx)) - - - if r!=0: - vec = dx/r - - else: - vec = np.zeros(len(p1)) - - return vec, r - -def getSrc_locs(DCsurvey): - """ - - - """ - - srcMat = np.zeros((DCsurvey.nSrc,2,3)) - for ii in range(DCsurvey.nSrc): - srcMat[ii,:,:] = np.asarray(DCsurvey.srcList[ii].loc) - - return srcMat diff --git a/SimPEG/DCIP/Utils.py b/SimPEG/DCIP/Utils.py deleted file mode 100644 index 702a2333..00000000 --- a/SimPEG/DCIP/Utils.py +++ /dev/null @@ -1,38 +0,0 @@ -import numpy as np - -def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False): - - import SimPEG.DCIP as DC - - elocs = np.arange(0,aSpacing*nElecs,aSpacing) - elocs -= (nElecs*aSpacing - aSpacing)/2 - space = 1 - WENNER = np.zeros((0,),dtype=int) - for ii in range(nElecs): - for jj in range(nElecs): - test = np.r_[jj,jj+space,jj+space*2,jj+space*3] - if np.any(test >= nElecs): - break - WENNER = np.r_[WENNER, test] - space += 1 - WENNER = WENNER.reshape((-1,4)) - - - if plotIt: - for i, s in enumerate('rbkg'): - plt.plot(elocs[WENNER[:,i]],s+'.') - plt.show() - - # Create sources and receivers - i = 0 - if in2D: - getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0] - else: - getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0] - srcList = [] - for i in range(WENNER.shape[0]): - rx = DC.RxDipole(getLoc(i,1),getLoc(i,2)) - src = DC.SrcDipole([rx], getLoc(i,0),getLoc(i,3)) - srcList += [src] - - return srcList diff --git a/SimPEG/DCIP/__init__.py b/SimPEG/DCIP/__init__.py deleted file mode 100644 index 08d1fe12..00000000 --- a/SimPEG/DCIP/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -from BaseDC import * -from BaseIP import * -from DCIPUtils import * -import Utils