first implementation of a checkerboard modification for a given vgrids.in file

This commit is contained in:
Marcel Paffrath 2015-11-09 16:12:12 +01:00
parent 2a15ec75d2
commit 7aedc35c16

View File

@ -511,6 +511,136 @@ class SeisArray(object):
print('Wrote %d points to file %s for %d layers'%(count, outfilename, nlayers))
outfile.close()
def addCheckerboard(self, spacing = 20, pertubation = 1, inputfile = 'vgrids.in', outputfile = 'vgrids_cb.in'):
def readNumberOfPoints(filename):
fin = open(filename, 'r')
vglines = fin.readlines()
nR = int(vglines[1].split()[0])
nTheta = int(vglines[1].split()[1])
nPhi = int(vglines[1].split()[2])
print('readNumberOf Points: Awaiting %d grid points in %s'
%(nR*nTheta*nPhi, filename))
fin.close()
return nR, nTheta, nPhi
def readDelta(filename):
fin = open(filename, 'r')
vglines = fin.readlines()
dR = float(vglines[2].split()[0])
dTheta = float(vglines[2].split()[1])
dPhi = float(vglines[2].split()[2])
fin.close()
return dR, dTheta, dPhi
def readStartpoints(filename):
fin = open(filename, 'r')
vglines = fin.readlines()
sR = float(vglines[3].split()[0])
sTheta = float(vglines[3].split()[1])
sPhi = float(vglines[3].split()[2])
fin.close()
return sR, sTheta, sPhi
def readVelocity(filename):
'''
Reads in velocity from vgrids file and returns a list containing all values in the same order
'''
vel = []; count = 0
fin = open(filename, 'r')
vglines = fin.readlines()
for line in vglines:
count += 1
if count > 4:
vel.append(float(line.split()[0]))
print("Read %d points out of file: %s" %(count - 4, filename))
return vel
def correctSpacing(spacing, delta, disttype = None):
if spacing > delta:
spacing_corr = round(spacing / delta) * delta
elif spacing < delta:
spacing_corr = delta
print('The spacing of the checkerboard of %s (%s) was corrected to '
'a value of %s to fit the grid spacing of %s.' %(spacing, disttype, spacing_corr, delta))
return spacing_corr
R = 6371. # earth radius
decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)
outfile = open(outputfile, 'w')
# Theta, Phi in radians, R in km
nR, nTheta, nPhi = readNumberOfPoints(inputfile)
dR, dThetaRad, dPhiRad = readDelta(inputfile)
sR, sThetaRad, sPhiRad = readStartpoints(inputfile)
vel = readVelocity(inputfile)
dTheta, dPhi = np.rad2deg((dThetaRad, dPhiRad))
sTheta, sPhi = np.rad2deg((dThetaRad, dPhiRad))
eR = sR + (nR - 1) * dR
ePhi = sPhi + (nPhi - 1) * dPhi
eTheta = sTheta + (nTheta - 1) * dTheta
nPoints = nR * nTheta * nPhi
thetaGrid = np.linspace(sTheta, eTheta, num = nTheta)
phiGrid = np.linspace(sPhi, ePhi, num = nPhi)
rGrid = np.linspace(sR, eR, num = nR)
# write header for velocity grid file (in RADIANS)
outfile.writelines('%10s %10s \n' %(1, 1))
outfile.writelines('%10s %10s %10s\n' %(nR, nTheta, nPhi))
outfile.writelines('%10s %10s %10s\n' %(dR, dThetaRad, dPhiRad))
outfile.writelines('%10s %10s %10s\n' %(sR, sThetaRad, sPhiRad))
spacR = correctSpacing(spacing, dR, '[meter], R')
spacTheta = correctSpacing(self._getAngle(spacing), dTheta, '[degree], Theta')
spacPhi = correctSpacing(self._getAngle(spacing), dPhi, '[degree], Phi')
count = 0
evenOdd = 1
even = 0; odd = 0
# In the following loop it is checked whether the positive distance from the border of the model
# for a point on the grid divided by the spacing is even or odd and then pertubated.
# The position is also shifted by half of the delta so that the position is directly on the point and
# not on the border between two points.
for radius in rGrid:
if np.floor((radius - sR - dR/2) / spacR) % 2:
evenOddR = 1
else:
evenOddR = -1
for theta in thetaGrid:
if np.floor((theta - sTheta - dTheta/2) / spacTheta) % 2:
evenOddT = 1
else:
evenOddT = -1
for phi in phiGrid:
if np.floor((phi - sPhi - dPhi/2) / spacPhi) % 2:
evenOddP = 1
else:
evenOddP = -1
velocity = vel[count]
evenOdd = evenOddR * evenOddT * evenOddP
velocity += evenOdd * pertubation
outfile.writelines('%10s %10s\n'%(velocity, decm))
count += 1
progress = float(count) / float(nPoints) * 100
self._update_progress(progress)
outfile.close()
def exportAll(self, filename = 'interpolated_receivers.out'):
recfile_out = open(filename, 'w')
count = 0