diff --git a/SimPEG/DCIP/BaseDC.py b/SimPEG/DCIP/BaseDC.py index 5fa73265..3761663c 100644 --- a/SimPEG/DCIP/BaseDC.py +++ b/SimPEG/DCIP/BaseDC.py @@ -92,7 +92,7 @@ class SurveyDC(Survey.BaseSurvey): Geophysical DC resistivity data. """ - + uncert = None def __init__(self, srcList, **kwargs): self.srcList = srcList Survey.BaseSurvey.__init__(self, **kwargs) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 15073000..39825ef2 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -2,7 +2,12 @@ from SimPEG import np import BaseDC as DC import BaseDC as IP -def genActiveindfromTopo(mesh, topo): +def getActiveindfromTopo(mesh, topo): +# def genActiveindfromTopo(mesh, topo): + """ + Get active indices from topography + """ + from scipy.interpolate import NearestNDInterpolator if mesh.dim==3: nCxy = mesh.nCx*mesh.nCy Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F') @@ -16,90 +21,97 @@ def genActiveindfromTopo(mesh, topo): else: raise NotImplementedError("Only 3D is working") - return actind + return Utils.mkvc(np.vstack(actind)) def gettopoCC(mesh, airind): - 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 topoCC +# def gettopoCC(mesh, airind): + """ + Get topography from active indices of mesh. + """ + 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_DC3Dobs(filename,mesh,topo,probType="CC"): - text_file = open("./mill_dc.dat", "r") - lines = text_file.readlines() - text_file.close() - SRC = [] - DATA = [] - srcLists = [] - isrc = 0 - airind = genActiveindfromTopo(mesh, topo) - topoCC, topoCCind = gettopoCC(mesh, airind) +def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"): + 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 = DC.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 = DC.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 = DC.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)]) + 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 = DC.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 = DC.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 = DC.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 = DC.SurveyDC(srcLists) + 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) + # 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} + # DCdata = Survey.Data(surveytest, surveytest.dobs) + # DCdata[src0, src0.rxList[0]] + return {'DCsurvey':survey, 'airind':airind, 'topoCC':topoCC, 'SRC':SRC} def readUBC_DC2DModel(fileName):