From 0a0ab4bfa9b8be749761696014f2baa184bb6d18 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 16 Nov 2015 11:48:45 +0100 Subject: [PATCH] added checkerboard to fmtomo utils and added different checkerboard ampl. functions --- pylot/core/active/fmtomoUtils.py | 417 +++++++++++++++++++++++++++++++ 1 file changed, 417 insertions(+) create mode 100644 pylot/core/active/fmtomoUtils.py diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py new file mode 100644 index 00000000..20586b16 --- /dev/null +++ b/pylot/core/active/fmtomoUtils.py @@ -0,0 +1,417 @@ +# -*- coding: utf-8 -*- +import sys +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 + ''' + 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 + + :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)) + +def addCheckerboard(spacing = 10., pertubation = 0.1, inputfile = 'vgrids.in', + outputfile = 'vgrids_cb.in', ampmethod = 'linear', rect = (None, None)): + ''' + Add a checkerboard to an existing vgrids.in velocity model. + + :param: spacing, size of the tiles + type: float + + :param: pertubation, pertubation (default: 0.1 = 10%) + type: float + ''' + 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 + + def linearAmp(InCell): + decimal = InCell - np.floor(InCell) + return (-abs(decimal - 0.5) + 0.5) * 2 + + def rectAmp(InCell, rect): + decimal = InCell - np.floor(InCell) + r1, r2 = rect + if r1 <= decimal <= r2: + return 1 + else: + return 0 + + def ampFunc(InCell, method = 'linear', rect = None): + if method == 'linear': + return linearAmp(InCell) + if method == 'rect' and rect is not None: + return rectAmp(InCell, rect) + else: + print('ampFunc: Could not amplify cb pattern') + + + 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(_getAngle(spacing), dTheta, '[degree], Theta') + spacPhi = correctSpacing(_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. + # "InCell" points e.g. rInCell are floats with their integer number corresponding to the cell number and + # their decimal place (0 - 1) corresponding to the position inside the cell. + # The amplification factor ampFactor comes from a linear relationship and ranges between 0 (cell border) + # and 1 (cell middle) + for radius in rGrid: + rInCell = (radius - sR - dR/2) / spacR + ampR = ampFunc(rInCell, ampmethod, rect) + if np.floor(rInCell) % 2: + evenOddR = 1 + else: + evenOddR = -1 + for theta in thetaGrid: + thetaInCell = (theta - sTheta - dTheta/2) / spacTheta + ampTheta = ampFunc(thetaInCell, ampmethod, rect) + if np.floor(thetaInCell) % 2: + evenOddT = 1 + else: + evenOddT = -1 + for phi in phiGrid: + phiInCell = (phi - sPhi - dPhi/2) / spacPhi + ampPhi = ampFunc(phiInCell, ampmethod, rect) + if np.floor(phiInCell) % 2: + evenOddP = 1 + else: + evenOddP = -1 + velocity = vel[count] + ampFactor = (ampR + ampTheta + ampPhi) / 3 + evenOdd = evenOddR * evenOddT * evenOddP * ampFactor + velocity += evenOdd * pertubation * velocity + + outfile.writelines('%10s %10s\n'%(velocity, decm)) + count += 1 + + progress = float(count) / float(nPoints) * 100 + _update_progress(progress) + + print('Added checkerboard to the grid in file %s with a spacing of %s and a pertubation of %s %%. ' + 'Outputfile: %s.'%(inputfile, spacing, pertubation, outputfile)) + outfile.close() + +def _update_progress(progress): + sys.stdout.write("%d%% done \r" % (progress) ) + sys.stdout.flush() + +def _getAngle(distance): + ''' + Function returns the angle on a Sphere of the radius R = 6371 [km] for a distance [km]. + ''' + PI = np.pi + R = 6371. + angle = distance * 180. / (PI * R) + return angle +