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

Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

Browse Source
This commit is contained in:
Marcel Paffrath 2015-11-05 12:20:12 +01:00
commit 2a15ec75d2
1 changed files with 19 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -427,7 +427,7 @@ class SeisArray(object):
infile.close() infile.close()
return ztop, zbot, vtop, vbot return ztop, zbot, vtop, vbot
R = 6371 R = 6371.
vmin = 0.34 vmin = 0.34
cushionfactor = 0.1 # add some extra space to the model cushionfactor = 0.1 # add some extra space to the model
decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3) decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)
@ -435,7 +435,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:
x, y, z = self.getAllMeasuredPointsLists() x, y, _ = self.getAllMeasuredPointsLists()
phi_min, phi_max = (self._getAngle(min(x)), self._getAngle(max(x))) phi_min, phi_max = (self._getAngle(min(x)), self._getAngle(max(x)))
theta_min, theta_max = (self._getAngle(min(y)), self._getAngle(max(y))) theta_min, theta_max = (self._getAngle(min(y)), self._getAngle(max(y)))
cushionPhi = abs(phi_max - phi_min) * cushionfactor cushionPhi = abs(phi_max - phi_min) * cushionfactor
@ -459,7 +459,7 @@ class SeisArray(object):
rGrid = np.linspace(rbot - rDelta, rtop + rDelta, num = nR + 2) # +2 cushion nodes rGrid = np.linspace(rbot - rDelta, rtop + rDelta, 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)
# 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))
@ -469,8 +469,10 @@ class SeisArray(object):
surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method, filename = None) surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method, filename = None)
print "\nGenerating velocity grid for FMTOMO. Output filename = %s, interpolation method = %s"%(outfilename, method) print("\nGenerating velocity grid for FMTOMO. "
print "nTheta = %s, nPhi = %s, nR = %s, thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt) "Output filename = %s, interpolation method = %s"%(outfilename, method))
print("nTheta = %s, nPhi = %s, nR = %s, "
"thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt))
count = 0 count = 0
nlayers = readMygridNlayers(infilename) nlayers = readMygridNlayers(infilename)
@ -484,22 +486,23 @@ class SeisArray(object):
for point in surface: for point in surface:
if point[0] == xval and point[1] == yval: if point[0] == xval and point[1] == yval:
topo = point[2] topo = point[2]
z = -(R + topo - radius) break
if z > (topo + 1): depth = -(R + topo - radius)
if depth > 1:
vel = 0.0 vel = 0.0
elif (topo + 1) >= z > topo: # cushioning around topography elif 1 >= depth > 0: # cushioning around topography
vel = vtop[0] vel = vtop[0]
else: else:
for index in range(nlayers): for index in range(nlayers):
if (topo + ztop[index]) >= z > (topo + zbot[index]): if (ztop[index]) >= depth > (zbot[index]):
vel = (z - ztop[index] - topo) / (zbot[index] - ztop[index]) * (vbot[index] - vtop[index]) + vtop[index] vel = (depth - ztop[index]) / (zbot[index] - ztop[index]) * (vbot[index] - vtop[index]) + vtop[index]
break break
if not (topo + ztop[index]) >= z > (topo + zbot[index]): if not (ztop[index]) >= depth > (zbot[index]):
print('ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile'%(z - topo)) print('ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile'%(depth-topo))
return return
count += 1 count += 1
if vel < 0: if vel < 0:
print('ERROR, vel <0; z, topo, zbot, vbot, vtop:', z, topo, zbot[index], vbot[index], vtop[index]) print('ERROR, vel <0; z, topo, zbot, vbot, vtop:', depth, topo, zbot[index], vbot[index], vtop[index])
outfile.writelines('%10s %10s\n'%(vel, decm)) outfile.writelines('%10s %10s\n'%(vel, decm))
progress = float(count) / float(nTotal) * 100 progress = float(count) / float(nTotal) * 100