marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

finished generatePropgrid, changed getThetaPhiFromArray output to the right order

Browse Source
This commit is contained in:
Marcel Paffrath 2015-11-12 11:25:51 +01:00
parent 39a12bd1d1
commit d611b8606e
1 changed files with 34 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
pylot/core/active/seismicArrayPreparation.py
View File

@ -358,11 +358,11 @@ class SeisArray(object):
measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists() measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists()
# need to determine the delta to add two cushion nodes around the min/max values # need to determine the delta to add two cushion nodes around the min/max values
thetaDelta = (thetaN - thetaS) / (nTheta - 1) deltaTheta = (thetaN - thetaS) / (nTheta - 1)
phiDelta = (phiE - phiW) / (nPhi - 1) deltaPhi = (phiE - phiW) / (nPhi - 1)
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes thetaGrid = np.linspace(thetaS - deltaTheta, thetaN + deltaTheta, num = nTheta + 2) # +2 cushion nodes
phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes phiGrid = np.linspace(phiW - deltaPhi, phiE + deltaPhi, num = nPhi + 2) # +2 cushion nodes
nTotal = len(thetaGrid) * len(phiGrid); count = 0 nTotal = len(thetaGrid) * len(phiGrid); count = 0
for theta in thetaGrid: for theta in thetaGrid:
@ -401,7 +401,7 @@ class SeisArray(object):
nInterfaces = 2 nInterfaces = 2
# generate dimensions of the grid from array # generate dimensions of the grid from array
phiWE, thetaSN = self.getThetaPhiFromArray(cushionfactor) thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor)
thetaS, thetaN = thetaSN thetaS, thetaN = thetaSN
phiW, phiE = phiWE phiW, phiE = phiWE
@ -410,14 +410,14 @@ class SeisArray(object):
outfile = open(outfilename, 'w') outfile = open(outfilename, 'w')
# determine the deltas # determine the deltas
thetaDelta = abs(thetaN - thetaS) / float((nTheta - 1)) deltaTheta = abs(thetaN - thetaS) / float((nTheta - 1))
phiDelta = abs(phiE - phiW) / float((nPhi - 1)) deltaPhi = abs(phiE - phiW) / float((nPhi - 1))
# write header for interfaces grid file (in RADIANS) # write header for interfaces grid file (in RADIANS)
outfile.writelines('%10s\n' %(nInterfaces)) outfile.writelines('%10s\n' %(nInterfaces))
outfile.writelines('%10s %10s\n' %(nTheta + 2, nPhi + 2)) # +2 cushion nodes outfile.writelines('%10s %10s\n' %(nTheta + 2, nPhi + 2)) # +2 cushion nodes
outfile.writelines('%10s %10s\n' %(np.deg2rad(thetaDelta), np.deg2rad(phiDelta))) outfile.writelines('%10s %10s\n' %(np.deg2rad(deltaTheta), np.deg2rad(deltaPhi)))
outfile.writelines('%10s %10s\n' %(np.deg2rad(thetaS - thetaDelta), np.deg2rad(phiW - phiDelta))) outfile.writelines('%10s %10s\n' %(np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi)))
interface1 = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method) interface1 = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method)
interface2 = self.interpolateOnRegularGrid(nTheta, nPhi, thetaSN, phiWE, -depthmax, method = method) interface2 = self.interpolateOnRegularGrid(nTheta, nPhi, thetaSN, phiWE, -depthmax, method = method)
@ -450,9 +450,9 @@ class SeisArray(object):
cushionTheta = abs(theta_max - theta_min) * cushionfactor cushionTheta = abs(theta_max - theta_min) * cushionfactor
phiWE = (phi_min - cushionPhi, phi_max + cushionPhi) phiWE = (phi_min - cushionPhi, phi_max + cushionPhi)
thetaSN = (theta_min - cushionTheta, theta_max + cushionTheta) thetaSN = (theta_min - cushionTheta, theta_max + cushionTheta)
return phiWE, thetaSN return thetaSN, phiWE
def generatePropgrid(self, nTheta, nPhi, nR, Rbt, cushionfactor = 0.1, def generatePropgrid(self, nTheta, nPhi, nR, Rbt, cushionpropgrid = 0.05,
refinement = (5, 5), outfilename = 'propgrid.in'): refinement = (5, 5), outfilename = 'propgrid.in'):
''' '''
Create a propergation grid file for FMTOMO using SeisArray boundaries Create a propergation grid file for FMTOMO using SeisArray boundaries
@ -469,7 +469,8 @@ class SeisArray(object):
:param: Rbt (bot, top) extensions of the model in km :param: Rbt (bot, top) extensions of the model in km
type: tuple type: tuple
:param: cushionfactor, add some extra space to the model (default: 0.1 = 10%) :param: cushionpropogrid, cushionfactor for the propagationgrid (cushion direction
opposing to vgrids cushionfactor)
type: float type: float
:param: refinement, (refinement factor, number of local cells for refinement) used by FMTOMO :param: refinement, (refinement factor, number of local cells for refinement) used by FMTOMO
@ -477,22 +478,22 @@ class SeisArray(object):
''' '''
outfile = open(outfilename, 'w') outfile = open(outfilename, 'w')
R = 6371. thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor = 0)
phiWE, thetaSN = self.getThetaPhiFromArray(cushionfactor) thetaS = thetaSN[0] + cushionpropgrid
thetaN = thetaSN[1] - cushionpropgrid
phiW = phiWE[0] + cushionpropgrid
phiE = phiWE[1] - cushionpropgrid
rbot = Rbt[0]
rtop = Rbt[1]
thetaS, thetaN = thetaSN deltaTheta = abs(thetaN - thetaS) / float(nTheta - 1)
phiW, phiE = phiWE deltaPhi = abs(phiE - phiW) / float(nPhi - 1)
rbot = Rbt[0] + R deltaR = abs(rbot - rtop) / float(nR - 1)
rtop = Rbt[1] + R
thetaDelta = abs(thetaN - thetaS) / float((nTheta - 1))
phiDelta = abs(phiE - phiW) / float((nPhi - 1))
rDelta = abs(rbot - rtop) / float((nR - 1))
outfile.writelines('%10s %10s %10s\n' %(nR, nTheta, nPhi)) outfile.writelines('%10s %10s %10s\n' %(nR, nTheta, nPhi))
outfile.writelines('%10s %10s %10s\n' %(rDelta, thetaDelta, phiDelta)) outfile.writelines('%10s %10s %10s\n' %(deltaR, deltaTheta, deltaPhi))
outfile.writelines('%10s %10s %10s\n' %(Rbt[1], thetaS, phiW)) outfile.writelines('%10s %10s %10s\n' %(rtop, thetaS, phiW))
outfile.writelines('%10s %10s\n' %refinement) outfile.writelines('%10s %10s\n' %refinement)
outfile.close() outfile.close()
@ -568,7 +569,7 @@ class SeisArray(object):
# generate dimensions of the grid from array # generate dimensions of the grid from array
if thetaSN is None and phiWE is None: if thetaSN is None and phiWE is None:
phiWE, thetaSN = self.getThetaPhiFromArray() thetaSN, phiWE = self.getThetaPhiFromArray()
thetaS, thetaN = thetaSN thetaS, thetaN = thetaSN
phiW, phiE = phiWE phiW, phiE = phiWE
@ -576,14 +577,14 @@ class SeisArray(object):
rtop = Rbt[1] + R rtop = Rbt[1] + R
# need to determine the delta to add two cushion nodes around the min/max values # need to determine the delta to add two cushion nodes around the min/max values
thetaDelta = abs(thetaN - thetaS) / float((nTheta - 1)) deltaTheta = abs(thetaN - thetaS) / float((nTheta - 1))
phiDelta = abs(phiE - phiW) / float((nPhi - 1)) deltaPhi = abs(phiE - phiW) / float((nPhi - 1))
rDelta = abs(rbot - rtop) / float((nR - 1)) deltaR = abs(rbot - rtop) / float((nR - 1))
# create a regular grid including +2 cushion nodes in every direction # create a regular grid including +2 cushion nodes in every direction
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes thetaGrid = np.linspace(thetaS - deltaTheta, thetaN + deltaTheta, num = nTheta + 2) # +2 cushion nodes
phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes phiGrid = np.linspace(phiW - deltaPhi, phiE + deltaPhi, num = nPhi + 2) # +2 cushion nodes
rGrid = np.linspace(rbot - rDelta, rtop + rDelta, num = nR + 2) # +2 cushion nodes rGrid = np.linspace(rbot - deltaR, rtop + deltaR, num = nR + 2) # +2 cushion nodes
nTotal = len(rGrid) * len(thetaGrid) * len(phiGrid) nTotal = len(rGrid) * len(thetaGrid) * len(phiGrid)
print("Total number of grid nodes: %s"%nTotal) print("Total number of grid nodes: %s"%nTotal)
@ -591,8 +592,8 @@ class SeisArray(object):
# write header for velocity grid file (in RADIANS) # write header for velocity grid file (in RADIANS)
outfile.writelines('%10s %10s \n' %(1, 1)) outfile.writelines('%10s %10s \n' %(1, 1))
outfile.writelines('%10s %10s %10s\n' %(nR + 2, nTheta + 2, nPhi + 2)) outfile.writelines('%10s %10s %10s\n' %(nR + 2, nTheta + 2, nPhi + 2))
outfile.writelines('%10s %10s %10s\n' %(rDelta, np.deg2rad(thetaDelta), np.deg2rad(phiDelta))) outfile.writelines('%10s %10s %10s\n' %(deltaR, np.deg2rad(deltaTheta), np.deg2rad(deltaPhi)))
outfile.writelines('%10s %10s %10s\n' %(rbot - rDelta, np.deg2rad(thetaS - thetaDelta), np.deg2rad(phiW - phiDelta))) outfile.writelines('%10s %10s %10s\n' %(rbot - deltaR, np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi)))
surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method) surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method)