changed structure -> all FMTOMO to vtk functions to one module
This commit is contained in:
parent
5b8e2da59e
commit
3473ce732c
@ -192,116 +192,3 @@ def findTracesInRanges(shot_dict, distancebin, pickbin):
|
|||||||
shots_found[shot.getShotnumber()].append(traceID)
|
shots_found[shot.getShotnumber()].append(traceID)
|
||||||
|
|
||||||
return shots_found
|
return shots_found
|
||||||
|
|
||||||
def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'):
|
|
||||||
def getDistance(angle):
|
|
||||||
PI = np.pi
|
|
||||||
R = 6371.
|
|
||||||
distance = angle / 180 * (PI * R)
|
|
||||||
return distance
|
|
||||||
|
|
||||||
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])
|
|
||||||
|
|
||||||
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):
|
|
||||||
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
|
|
||||||
|
|
||||||
R = 6371 # earth radius
|
|
||||||
outfile = open(outputfile, 'w')
|
|
||||||
|
|
||||||
# Theta, Phi in radians, R in km
|
|
||||||
nR, nTheta, nPhi = readNumberOfPoints(inputfile)
|
|
||||||
dR, dTheta, dPhi = readDelta(inputfile)
|
|
||||||
sR, sTheta, sPhi = readStartpoints(inputfile)
|
|
||||||
vel = readVelocity(inputfile)
|
|
||||||
|
|
||||||
nX = nPhi; nY = nTheta; nZ = nR
|
|
||||||
|
|
||||||
sZ = sR - R
|
|
||||||
sX = getDistance(np.rad2deg(sPhi))
|
|
||||||
sY = getDistance(np.rad2deg(sTheta))
|
|
||||||
|
|
||||||
dX = getDistance(np.rad2deg(dPhi))
|
|
||||||
dY = getDistance(np.rad2deg(dTheta))
|
|
||||||
|
|
||||||
xGrid = np.linspace(sX, sX + (dX * nX), nX)
|
|
||||||
yGrid = np.linspace(sZ, sZ + (nY * dY), nY)
|
|
||||||
zGrid = np.linspace(sZ, sZ + (nR * dR), nR)
|
|
||||||
nPoints = len(xGrid) * len(yGrid) * len(zGrid)
|
|
||||||
|
|
||||||
# write header
|
|
||||||
print("Writing header for VTK file...")
|
|
||||||
outfile.writelines('# vtk DataFile Version 3.1\n')
|
|
||||||
outfile.writelines('Velocity on FMTOMO vgrids.in points\n')
|
|
||||||
outfile.writelines('ASCII\n')
|
|
||||||
outfile.writelines('DATASET POLYDATA\n')
|
|
||||||
outfile.writelines('POINTS %15d float\n' %(nPoints))
|
|
||||||
|
|
||||||
# write coordinates
|
|
||||||
print("Writing coordinates to VTK file...")
|
|
||||||
for z in zGrid:
|
|
||||||
for y in yGrid:
|
|
||||||
for x in xGrid:
|
|
||||||
outfile.writelines('%10f %10f %10f \n' %(x, y, z))
|
|
||||||
|
|
||||||
|
|
||||||
outfile.writelines('VERTICES %15d %15d\n' %(nPoints, 2 * nPoints))
|
|
||||||
|
|
||||||
# write indices
|
|
||||||
print("Writing indices to VTK file...")
|
|
||||||
for index in range(nPoints):
|
|
||||||
outfile.writelines('%10d %10d\n' %(1, index))
|
|
||||||
|
|
||||||
outfile.writelines('POINT_DATA %15d\n' %(nPoints))
|
|
||||||
outfile.writelines('SCALARS velocity float %d\n' %(1))
|
|
||||||
outfile.writelines('LOOKUP_TABLE default\n')
|
|
||||||
|
|
||||||
# write velocity
|
|
||||||
print("Writing velocity values to VTK file...")
|
|
||||||
for velocity in vel:
|
|
||||||
outfile.writelines('%10f\n' %velocity)
|
|
||||||
|
|
||||||
outfile.close()
|
|
||||||
print("Wrote velocity grid for %d points to file: %s" %(nPoints, outputfile))
|
|
||||||
return
|
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user