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

implemented a function that generates all grids for FMTOMO

Browse Source
This commit is contained in:
Marcel Paffrath 2015-11-12 13:24:48 +01:00
parent d611b8606e
commit d53e4b7c0c
1 changed files with 66 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -381,6 +381,34 @@ class SeisArray(object):
return surface return surface
def generateFMTOMOinputFromArray(self, nRP, nThetaP, nPhiP, nRI, nThetaI, nPhiI,
Rbt, cushionfactor, interpolationMethod = 'linear',
customgrid = 'mygrid.in'):
print('\n------------------------------------------------------------')
print('Automatically generating input for FMTOMO from array size.')
print('Propgrid: nR = %s, nTheta = %s, nPhi = %s'%(nRP, nThetaP, nPhiP))
print('Interpolation Grid and Interfaces Grid: nR = %s, nTheta = %s, nPhi = %s'%(nRI, nThetaI, nPhiI))
print('Bottom and Top of model: (%s, %s)'%(Rbt[0], Rbt[1]))
print('Method: %s, customgrid = %s'%(interpolationMethod, customgrid))
print('------------------------------------------------------------')
def getZmin(surface):
z = []
for point in surface:
z.append(point[2])
return min(z)
self.generatePropgrid(nThetaP, nPhiP, nRP, Rbt, cushionpropgrid = 0.05)
surface = self.generateVgrid(nThetaI, nPhiI, nRI, Rbt, method = interpolationMethod,
cushionfactor = cushionfactor, infilename = customgrid,
returnTopo = True)
depthmax = abs(Rbt[0] - getZmin(surface)) - 1.0 # cushioning for the bottom interface
self.generateInterfaces(nThetaI, nPhiI, depthmax, cushionfactor = cushionfactor,
returnInterfaces = False, method = interpolationMethod)
def generateInterfaces(self, nTheta, nPhi, depthmax, cushionfactor = 0.1, def generateInterfaces(self, nTheta, nPhi, depthmax, cushionfactor = 0.1,
outfilename = 'interfaces.in', method = 'linear', outfilename = 'interfaces.in', method = 'linear',
returnInterfaces = False): returnInterfaces = False):
@ -398,6 +426,9 @@ class SeisArray(object):
:param: cushionfactor, add some extra space to the model (default: 0.1 = 10%) :param: cushionfactor, add some extra space to the model (default: 0.1 = 10%)
type: float type: float
''' '''
print('\n------------------------------------------------------------')
print('Generating interfaces...')
nInterfaces = 2 nInterfaces = 2
# generate dimensions of the grid from array # generate dimensions of the grid from array
@ -436,6 +467,9 @@ class SeisArray(object):
if returnInterfaces == True: if returnInterfaces == True:
return interface1, interface2 return interface1, interface2
print('Finished generating interfaces.')
print('------------------------------------------------------------')
def getThetaPhiFromArray(self, cushionfactor = 0.1): def getThetaPhiFromArray(self, cushionfactor = 0.1):
''' '''
Determine and returns PhiWE (tuple: (West, East)) and thetaSN (tuple (South, North)) from the SeisArray boundaries. Determine and returns PhiWE (tuple: (West, East)) and thetaSN (tuple (South, North)) from the SeisArray boundaries.
@ -466,7 +500,7 @@ class SeisArray(object):
:param: nR, number of points in R :param: nR, number of points in R
type: int type: int
:param: Rbt (bot, top) extensions of the model in km :param: Rbt (bot, top) extensions of the model in m
type: tuple type: tuple
:param: cushionpropogrid, cushionfactor for the propagationgrid (cushion direction :param: cushionpropogrid, cushionfactor for the propagationgrid (cushion direction
@ -478,6 +512,13 @@ class SeisArray(object):
''' '''
outfile = open(outfilename, 'w') outfile = open(outfilename, 'w')
print('\n------------------------------------------------------------')
print('Generating Propagation Grid for nTheta = %s, nPhi'
' = %s, nR = %s and a cushioning of %s'
%(nTheta, nPhi, nR, cushionpropgrid))
print('Bottom of the grid: %s, top of the grid %s'
%(Rbt[0], Rbt[1]))
thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor = 0) thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor = 0)
thetaS = thetaSN[0] + cushionpropgrid thetaS = thetaSN[0] + cushionpropgrid
@ -498,10 +539,13 @@ class SeisArray(object):
outfile.close() outfile.close()
def generateVgrid(self, nTheta = 80, nPhi = 80, nR = 120, print('Created Propagation Grid and saved it to %s' %outfilename)
Rbt = (-62.0, 6.0), thetaSN = None, print('------------------------------------------------------------')
phiWE = None, outfilename = 'vgrids.in',
method = 'linear', infilename = 'mygrid.in'): def generateVgrid(self, nTheta, nPhi, nR, Rbt, thetaSN = None,
phiWE = None, cushionfactor = 0.1,
outfilename = 'vgrids.in', method = 'linear',
infilename = 'mygrid.in', returnTopo = False):
''' '''
Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s]. Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s].
@ -520,7 +564,7 @@ class SeisArray(object):
:param: phiWE (W, E) extensions of the model in degree :param: phiWE (W, E) extensions of the model in degree
type: tuple type: tuple
:param: Rbt (bot, top) extensions of the model in km :param: Rbt (bot, top) extensions of the model in m
type: tuple type: tuple
:param: vbot, velocity at the bottom of the model :param: vbot, velocity at the bottom of the model
@ -529,6 +573,8 @@ class SeisArray(object):
:param: method, interpolation method for topography :param: method, interpolation method for topography
type: str type: str
''' '''
print('\n------------------------------------------------------------')
print('generateVgrid: Starting...')
# def getRad(angle): # def getRad(angle):
# PI = np.pi # PI = np.pi
@ -547,6 +593,7 @@ class SeisArray(object):
infile = open(filename, 'r') infile = open(filename, 'r')
nlayers = readMygridNlayers(filename) nlayers = readMygridNlayers(filename)
print('\nreadMygrid: Reading file %s.'%filename)
for index in range(nlayers): for index in range(nlayers):
line1 = infile.readline() line1 = infile.readline()
line2 = infile.readline() line2 = infile.readline()
@ -554,6 +601,10 @@ class SeisArray(object):
vtop.append(float(line1.split()[1])) vtop.append(float(line1.split()[1]))
zbot.append(float(line2.split()[0])) zbot.append(float(line2.split()[0]))
vbot.append(float(line2.split()[1])) vbot.append(float(line2.split()[1]))
print('Layer %s:\n[Top: v = %s [km/s], z = %s [m]]'
'\n[Bot: v = %s [km/s], z = %s [m]]'
%(index + 1, vtop[index], ztop[index],
vbot[index], zbot[index]))
if not ztop[0] == 0: if not ztop[0] == 0:
print('ERROR: there must be a velocity set for z = 0 in the file %s'%filename) print('ERROR: there must be a velocity set for z = 0 in the file %s'%filename)
@ -569,7 +620,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:
thetaSN, phiWE = self.getThetaPhiFromArray() thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor)
thetaS, thetaN = thetaSN thetaS, thetaN = thetaSN
phiW, phiE = phiWE phiW, phiE = phiWE
@ -597,15 +648,15 @@ class SeisArray(object):
surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method) surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method)
nlayers = readMygridNlayers(infilename)
ztop, zbot, vtop, vbot = readMygrid(infilename)
print("\nGenerating velocity grid for FMTOMO. " print("\nGenerating velocity grid for FMTOMO. "
"Output filename = %s, interpolation method = %s"%(outfilename, method)) "Output filename = %s, interpolation method = %s"%(outfilename, method))
print("nTheta = %s, nPhi = %s, nR = %s, " print("nTheta = %s, nPhi = %s, nR = %s, "
"thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt)) "thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt))
count = 0 count = 0
nlayers = readMygridNlayers(infilename)
ztop, zbot, vtop, vbot = readMygrid(infilename)
for radius in rGrid: for radius in rGrid:
for theta in thetaGrid: for theta in thetaGrid:
for phi in phiGrid: for phi in phiGrid:
@ -636,9 +687,13 @@ class SeisArray(object):
progress = float(count) / float(nTotal) * 100 progress = float(count) / float(nTotal) * 100
self._update_progress(progress) self._update_progress(progress)
print('Wrote %d points to file %s for %d layers'%(count, outfilename, nlayers)) print('\nWrote %d points to file %s for %d layers'%(count, outfilename, nlayers))
print('------------------------------------------------------------')
outfile.close() outfile.close()
if returnTopo == True:
return surface
def addCheckerboard(self, spacing = 20., pertubation = 0.1, inputfile = 'vgrids.in', outputfile = 'vgrids_cb.in'): def addCheckerboard(self, spacing = 20., pertubation = 0.1, inputfile = 'vgrids.in', outputfile = 'vgrids_cb.in'):
''' '''
Add a checkerboard to an existing vgrids.in velocity model. Add a checkerboard to an existing vgrids.in velocity model.