wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 20:23:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

return DCdata class?

Browse Source
This commit is contained in:
seogi_macbook
2016-02-11 15:29:16 -08:00
parent d9c94b9793
commit c1bd8edf77
2 changed files with 91 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/BaseDC.py
+1 -1
View File
@@ -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)
SimPEG/DCIP/DCIPUtils.py
+90 -78
View File
@@ -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):