Merged branch develop into master

This commit is contained in:
Marcel Paffrath 2016-05-02 11:31:17 +02:00
commit dbe47690eb
2 changed files with 25 additions and 243 deletions

View File

@ -1,243 +0,0 @@
# -*- coding: utf-8 -*-
import numpy as np
def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', inputfileref='vgridsref.in'):
'''
Generate a vtk-file readable by e.g. paraview from FMTOMO output vgrids.in
:param: absOrRel, can be "abs" or "rel" for absolute or relative velocities. if "rel" inputfileref must be given
:type: str
'''
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])
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
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))
nPoints = nX * nY * nZ
dZ = dR
# 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 STRUCTURED_POINTS\n')
outfile.writelines('DIMENSIONS %d %d %d\n' % (nX, nY, nZ))
outfile.writelines('ORIGIN %f %f %f\n' % (sX, sY, sZ))
outfile.writelines('SPACING %f %f %f\n' % (dX, dY, dZ))
outfile.writelines('POINT_DATA %15d\n' % (nPoints))
if absOrRel == 'abs':
outfile.writelines('SCALARS velocity float %d\n' % (1))
elif absOrRel == 'rel':
outfile.writelines('SCALARS velChangePercent float %d\n' % (1))
outfile.writelines('LOOKUP_TABLE default\n')
# write velocity
if absOrRel == 'abs':
print("Writing velocity values to VTK file...")
for velocity in vel:
outfile.writelines('%10f\n' % velocity)
elif absOrRel == 'rel':
velref = readVelocity(inputfileref)
if not len(velref) == len(vel):
print('ERROR: Number of gridpoints mismatch for %s and %s' % (inputfile, inputfileref))
return
# velrel = [((vel - velref) / velref * 100) for vel, velref in zip(vel, velref)]
velrel = []
for velocities in zip(vel, velref):
v, vref = velocities
if not vref == 0:
velrel.append((v - vref) / vref * 100)
else:
velrel.append(0)
nR_ref, nTheta_ref, nPhi_ref = readNumberOfPoints(inputfileref)
if not nR_ref == nR and nTheta_ref == nTheta and nPhi_ref == nPhi:
print('ERROR: Dimension mismatch of grids %s and %s' % (inputfile, inputfileref))
return
print("Writing velocity values to VTK file...")
for velocity in velrel:
outfile.writelines('%10f\n' % velocity)
print('Pertubations: min: %s, max: %s' % (min(velrel), max(velrel)))
outfile.close()
print("Wrote velocity grid for %d points to file: %s" % (nPoints, outputfile))
return
def rays2VTK(fnin, fdirout='./vtk_files/', nthPoint=50):
'''
Writes VTK file(s) for FMTOMO rays from rays.dat. There is one file created for each ray.
:param: fdirout, output directory, must exist before
:type: str
:param: nthPoint, plot every nth point of the ray
:type: integer
'''
def getDistance(angle):
PI = np.pi
R = 6371.
distance = angle / 180 * (PI * R)
return distance
infile = open(fnin, 'r')
R = 6371
rays = {}
raynumber = 0
nPoints = 0
### NOTE: rays.dat seems to be in km and radians
while True:
raynumber += 1
firstline = infile.readline()
if firstline == '': break # break at EOF
raynumber = int(firstline.split()[0])
shotnumber = int(firstline.split()[1])
rayValid = int(firstline.split()[4]) # is zero if the ray is invalid
if rayValid == 0:
print('Invalid ray number %d for shot number %d' % (raynumber, shotnumber))
continue
nRayPoints = int(infile.readline().split()[0])
if not shotnumber in rays.keys():
rays[shotnumber] = {}
rays[shotnumber][raynumber] = []
for index in range(nRayPoints):
if index % nthPoint is 0 or index == (nRayPoints - 1):
rad, lat, lon = infile.readline().split()
rays[shotnumber][raynumber].append(
[getDistance(np.rad2deg(float(lon))), getDistance(np.rad2deg(float(lat))), float(rad) - R])
else:
dummy = infile.readline()
infile.close()
for shotnumber in rays.keys():
fnameout = fdirout + 'rays%03d.vtk' % (shotnumber)
outfile = open(fnameout, 'w')
nPoints = 0
for raynumber in rays[shotnumber]:
for ray in rays[shotnumber][raynumber]:
nPoints += 1
# write header
# print("Writing header for VTK file...")
print("Writing shot %d to file %s" % (shotnumber, fnameout))
outfile.writelines('# vtk DataFile Version 3.1\n')
outfile.writelines('FMTOMO rays\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 raynumber in rays[shotnumber].keys():
for raypoint in rays[shotnumber][raynumber]:
outfile.writelines('%10f %10f %10f \n' % (raypoint[0], raypoint[1], raypoint[2]))
outfile.writelines('LINES %15d %15d\n' % (len(rays[shotnumber]), len(rays[shotnumber]) + nPoints))
# write indices
# print("Writing indices to VTK file...")
count = 0
for raynumber in rays[shotnumber].keys():
outfile.writelines('%d ' % (len(rays[shotnumber][raynumber])))
for index in range(len(rays[shotnumber][raynumber])):
outfile.writelines('%d ' % (count))
count += 1
outfile.writelines('\n')
# outfile.writelines('POINT_DATA %15d\n' %(nPoints))
# outfile.writelines('SCALARS rays 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))

View File

@ -45,16 +45,41 @@ def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', i
outfile.writelines('POINT_DATA %15d\n' % (nPoints)) outfile.writelines('POINT_DATA %15d\n' % (nPoints))
if absOrRel == 'abs': if absOrRel == 'abs':
<<<<<<< HEAD
outfile.writelines('SCALARS velocity float %d\n' %(1)) outfile.writelines('SCALARS velocity float %d\n' %(1))
if absOrRel == 'relDepth':
outfile.writelines('SCALARS velocity2depthMean float %d\n' %(1))
=======
outfile.writelines('SCALARS velocity float %d\n' % (1))
>>>>>>> 37f9292c39246b327d3630995ca2521725c6cdd7
elif absOrRel == 'rel': elif absOrRel == 'rel':
outfile.writelines('SCALARS velChangePercent float %d\n' % (1)) outfile.writelines('SCALARS velChangePercent float %d\n' % (1))
outfile.writelines('LOOKUP_TABLE default\n') outfile.writelines('LOOKUP_TABLE default\n')
pointsPerR = nTheta * nPhi
# write velocity # write velocity
if absOrRel == 'abs': if absOrRel == 'abs':
print("Writing velocity values to VTK file...") print("Writing velocity values to VTK file...")
for velocity in vel: for velocity in vel:
<<<<<<< HEAD
outfile.writelines('%10f\n' %velocity) outfile.writelines('%10f\n' %velocity)
elif absOrRel == 'relDepth':
print("Writing velocity values to VTK file relative to mean of each depth...")
index = 0; count = 0
veldepth = []
for velocity in vel:
count += 1
veldepth.append(velocity)
if count%pointsPerR == 0:
velmean = np.mean(veldepth)
#print velmean, count, count/pointsPerR
for vel in veldepth:
outfile.writelines('%10f\n' %(vel - velmean))
veldepth = []
=======
outfile.writelines('%10f\n' % velocity)
>>>>>>> 37f9292c39246b327d3630995ca2521725c6cdd7
elif absOrRel == 'rel': elif absOrRel == 'rel':
nref, dref, sref, velref = _readVgrid(inputfileref) nref, dref, sref, velref = _readVgrid(inputfileref)
nR_ref, nTheta_ref, nPhi_ref = nref nR_ref, nTheta_ref, nPhi_ref = nref