diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py index 64d81a05..0c766166 100644 --- a/pylot/core/active/surveyUtils.py +++ b/pylot/core/active/surveyUtils.py @@ -192,116 +192,3 @@ def findTracesInRanges(shot_dict, distancebin, pickbin): shots_found[shot.getShotnumber()].append(traceID) 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 -