From 376d1cc6f8f2e6162167851c577b9ab0ecd031b5 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Thu, 1 Oct 2015 13:14:44 +0200 Subject: [PATCH] name change --- pylot/core/active/seismicArrayPreparation.py | 665 +++++++++++++++++++ 1 file changed, 665 insertions(+) create mode 100644 pylot/core/active/seismicArrayPreparation.py diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py new file mode 100644 index 00000000..ae1b42ec --- /dev/null +++ b/pylot/core/active/seismicArrayPreparation.py @@ -0,0 +1,665 @@ +import sys +import numpy as np +from scipy.interpolate import griddata + +class SeisArray(object): + ''' + Can be used to interpolate missing values of a receiver grid, if only support points were measured. + Input file should contain in each line: ('traceID' 'receiverLineID' 'number of the geophone on recLine' 'X' 'Y' 'Z') + + Can be used to generate a velocity grid file (vgrids.in) for FMTOMO with a topography adapting gradient. + + Can be used to generate an interface file for FMTOMO (right now only interface.z used by grid3dg) for the topography. + + Supports vtk output for sources and receivers. + Note: Source and Receiver files for FMTOMO will be generated by the Survey object (because traveltimes will be added directly). + ''' + def __init__(self, recfile): + self._receiverlines = {} + self._receiverCoords = {} + self._measuredReceivers = {} + self._measuredTopo = {} + self._sourceLocs = {} + self._geophoneNumbers = {} + self._receiverlist = open(recfile, 'r').readlines() + self._generateReceiverlines() + self._setReceiverCoords() + self._setGeophoneNumbers() + + def _generateReceiverlines(self): + ''' + Connects the traceIDs to the lineIDs + for each receiverline in a dictionary. + ''' + for receiver in self._receiverlist: + traceID = int(receiver.split()[0]) + lineID = int(receiver.split()[1]) + if not lineID in self._receiverlines.keys(): + self._receiverlines[lineID] = [] + self._receiverlines[lineID].append(traceID) + + def _setReceiverCoords(self): + ''' + Fills the three x, y, z dictionaries with measured coordinates + ''' + for line in self._getReceiverlist(): + traceID = int(line.split()[0]) + x = float(line.split()[3]) + y = float(line.split()[4]) + z = float(line.split()[5]) + self._receiverCoords[traceID] = (x, y, z) + self._measuredReceivers[traceID] = (x, y, z) + + def _setGeophoneNumbers(self): + for line in self._getReceiverlist(): + traceID = int(line.split()[0]) + gphoneNum = float(line.split()[2]) + self._geophoneNumbers[traceID] = gphoneNum + + def _getReceiverlines(self): + return self._receiverlines + + def _getReceiverlist(self): + return self._receiverlist + + def getReceiverCoordinates(self): + return self._receiverCoords + + def _getXreceiver(self, traceID): + return self._receiverCoords[traceID][0] + + def _getYreceiver(self, traceID): + return self._receiverCoords[traceID][1] + + def _getZreceiver(self, traceID): + return self._receiverCoords[traceID][2] + + def _getXshot(self, shotnumber): + return self._sourceLocs[shotnumber][0] + + def _getYshot(self, shotnumber): + return self._sourceLocs[shotnumber][1] + + def _getZshot(self, shotnumber): + return self._sourceLocs[shotnumber][2] + + def _getReceiverValue(self, traceID, coordinate): + setCoordinate = {'X': self._getXreceiver, + 'Y': self._getYreceiver, + 'Z': self._getZreceiver} + return setCoordinate[coordinate](traceID) + + def _getGeophoneNumber(self, traceID): + return self._geophoneNumbers[traceID] + + def getMeasuredReceivers(self): + return self._measuredReceivers + + def getMeasuredTopo(self): + return self._measuredTopo + + def getSourceLocations(self): + return self._sourceLocs + + def _setXvalue(self, traceID, value): + self._checkKey(traceID) + self._receiverCoords[traceID][0] = value + + def _setYvalue(self, traceID, value): + self._checkKey(traceID) + self._receiverCoords[traceID][1] = value + + def _setZvalue(self, traceID, value): + self._checkKey(traceID) + self._receiverCoords[traceID][2] = value + + def _setValue(self, traceID, coordinate, value): + setCoordinate = {'X': self._setXvalue, + 'Y': self._setYvalue, + 'Z': self._setZvalue} + setCoordinate[coordinate](traceID, value) + + def _checkKey(self, traceID): + if not traceID in self._receiverCoords.keys(): + self._receiverCoords[traceID] = [None, None, None] + + def _checkTraceIDdirection(self, traceID1, traceID2): + if traceID2 > traceID1: + direction = +1 + return direction + if traceID2 < traceID1: + direction = -1 + return direction + print "Error: Same Value for traceID1 = %s and traceID2 = %s" %(traceID1, traceID2) + + def _checkCoordDirection(self, traceID1, traceID2, coordinate): + ''' + Checks whether the interpolation direction is positive or negative + ''' + if self._getReceiverValue(traceID1, coordinate) < self._getReceiverValue(traceID2, coordinate): + direction = +1 + return direction + if self._getReceiverValue(traceID1, coordinate) > self._getReceiverValue(traceID2, coordinate): + direction = -1 + return direction + print "Error: Same Value for traceID1 = %s and traceID2 = %s" %(traceID1, traceID2) + + def _interpolateMeanDistances(self, traceID1, traceID2, coordinate): + ''' + Returns the mean distance between two traceID's depending on the number of geophones in between + ''' + num_spaces = abs(self._getGeophoneNumber(traceID1) - self._getGeophoneNumber(traceID2)) + mean_distance = abs(self._getReceiverValue(traceID1, coordinate) - self._getReceiverValue(traceID2, coordinate))/num_spaces + return mean_distance + + def interpolateValues(self, coordinate): + ''' + Interpolates and sets all values (linear) for coordinate = 'X', 'Y' or 'Z' + ''' + for lineID in self._getReceiverlines().keys(): + number_measured = len(self._getReceiverlines()[lineID]) + for index, traceID1 in enumerate(self._getReceiverlines()[lineID]): + if index + 1 < number_measured: + traceID2 = self._getReceiverlines()[lineID][index + 1] + + traceID_dir = self._checkTraceIDdirection(traceID1, traceID2) + traceID_interp = traceID1 + traceID_dir + + coord_dir = self._checkCoordDirection(traceID1, traceID2, coordinate) + mean_distance = self._interpolateMeanDistances(traceID1, traceID2, coordinate) + + while (traceID_dir * traceID_interp) < (traceID_dir * traceID2): + self._setValue(traceID_interp, coordinate, + (self._getReceiverValue(traceID_interp - traceID_dir, coordinate) + + coord_dir * (mean_distance))) + traceID_interp += traceID_dir + + def addMeasuredTopographyPoints(self, filename): + ''' + Use more measured points for a higher precision of height interpolation. + Input file should contain in each line: ('point ID' 'X' 'Y' 'Z') + ''' + topolist = open(filename, 'r').readlines() + for line in topolist: + line = line.split() + pointID = int(line[0]) + x = float(line[1]) + y = float(line[2]) + z = float(line[3]) + self._measuredTopo[pointID] = (x, y, z) + + def addSourceLocations(self, filename): + ''' + Use source locations for a higher precision of height interpolation. + Input file should contain in each line: ('point ID' 'X' 'Y' 'Z') + + Source locations must be added before they can be written to vtk files. + ''' + topolist = open(filename, 'r').readlines() + for line in topolist: + line = line.split() + pointID = int(line[0]) + x = float(line[1]) + y = float(line[2]) + z = float(line[3]) + self._sourceLocs[pointID] = (x, y, z) + + def interpZcoords4rec(self, method = 'linear'): + ''' + Interpolates z values for all receivers. + ''' + measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists() + + for traceID in self.getReceiverCoordinates().keys(): + if type(self.getReceiverCoordinates()[traceID]) is not tuple: + z = griddata((measured_x, measured_y), measured_z, (self._getXreceiver(traceID), self._getYreceiver(traceID)), method = method) + self._setZvalue(traceID, float(z)) + + def _getAngle(self, 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 + + def _getDistance(self, angle): + ''' + Function returns the distance [km] on a Sphere of the radius R = 6371 [km] for an angle. + ''' + PI = np.pi + R = 6371. + distance = angle / 180 * (PI * R) + return distance + + def getMeasuredReceiverLists(self): + ''' + Returns a list of all measured receivers known to SeisArray. + ''' + x = []; y = []; z = [] + for traceID in self.getMeasuredReceivers().keys(): + x.append(self.getMeasuredReceivers()[traceID][0]) + y.append(self.getMeasuredReceivers()[traceID][1]) + z.append(self.getMeasuredReceivers()[traceID][2]) + return x, y, z + + def getMeasuredTopoLists(self): + ''' + Returns a list of all measured topography points known to the SeisArray. + ''' + x = []; y = []; z = [] + for pointID in self.getMeasuredTopo().keys(): + x.append(self.getMeasuredTopo()[pointID][0]) + y.append(self.getMeasuredTopo()[pointID][1]) + z.append(self.getMeasuredTopo()[pointID][2]) + return x, y, z + + def getSourceLocsLists(self): + ''' + Returns a list of all measured source locations known to SeisArray. + ''' + x = []; y = []; z = [] + for pointID in self.getSourceLocations().keys(): + x.append(self.getSourceLocations()[pointID][0]) + y.append(self.getSourceLocations()[pointID][1]) + z.append(self.getSourceLocations()[pointID][2]) + return x, y, z + + def getAllMeasuredPointsLists(self): + ''' + Returns a list of all measured points known to SeisArray. + ''' + mtopo_x, mtopo_y, mtopo_z = self.getMeasuredTopoLists() + msource_x, msource_y, msource_z = self.getSourceLocsLists() + mrec_x, mrec_y, mrec_z = self.getMeasuredReceiverLists() + + x = mtopo_x + mrec_x + msource_x + y = mtopo_y + mrec_y + msource_y + z = mtopo_z + mrec_z + msource_z + return x, y, z + + def getReceiverLists(self): + ''' + Returns a list of all receivers (measured and interpolated). + ''' + x = []; y =[]; z = [] + for traceID in self.getReceiverCoordinates().keys(): + x.append(self.getReceiverCoordinates()[traceID][0]) + y.append(self.getReceiverCoordinates()[traceID][1]) + z.append(self.getReceiverCoordinates()[traceID][2]) + return x, y, z + + def _interpolateXY4rec(self): + ''' + Interpolates the X and Y coordinates for all receivers. + ''' + for coordinate in ('X', 'Y'): + self.interpolateValues(coordinate) + + def interpolateAll(self): + self._interpolateXY4rec() + self.interpZcoords4rec() + + def interpolateTopography(self, nTheta, nPhi, thetaSN, phiWE, method = 'linear', filename = 'interface1.in'): + ''' + Interpolate Z values on a regular grid with cushion nodes to use it as FMTOMO topography interface. + Returns a surface in form of a list of points [[x1, y1, z1], [x2, y2, y2], ...] (cartesian). + + :param: nTheta, number of points in theta (NS) + type: integer + + :param: nPhi, number of points in phi (WE) + type: integer + + :param: thetaSN (S, N) extensions of the model in degree + type: tuple + + :param: phiWE (W, E) extensions of the model in degree + type: tuple + ''' + + surface = [] + elevation = 0.25 # elevate topography so that no source lies above the surface + + if filename is not None: + outfile = open(filename, 'w') + + print "Interpolating topography on regular grid with the dimensions:" + print "nTheta = %s, nPhi = %s, thetaSN = %s, phiWE = %s"%(nTheta, nPhi, thetaSN, phiWE) + print "method = %s, filename = %s" %(method, filename) + + thetaS, thetaN = thetaSN + phiW, phiE = phiWE + + measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists() + + # need to determine the delta to add two cushion nodes around the min/max values + thetaDelta = (thetaN - thetaS) / (nTheta - 1) + phiDelta = (phiE - phiW) / (nPhi - 1) + + thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes + phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes + + nTotal = len(thetaGrid) * len(phiGrid); count = 0 + for theta in thetaGrid: + for phi in phiGrid: + xval = self._getDistance(phi) + yval = self._getDistance(theta) + z = griddata((measured_x, measured_y), measured_z, (xval, yval), method = method) + # in case the point lies outside, nan will be returned. Find nearest: + if np.isnan(z) == True: + z = griddata((measured_x, measured_y), measured_z, (xval, yval), method = 'nearest') + z = float(z) + surface.append((xval, yval, z)) + count += 1 + progress = float(count) / float(nTotal) * 100 + self._update_progress(progress) + + if filename is not None: + outfile.writelines('%10s\n'%(z + elevation)) + + return surface + + def generateVgrid(self, nTheta = 80, nPhi = 80, nR = 120, + thetaSN = (-0.2, 1.2), phiWE = (-0.2, 1.2), + Rbt = (-62.0, 6.0), vbot = 5.5, filename = 'vgrids.in', + method = 'linear' ): + ''' + Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s]. + + :param: nTheta, number of points in theta (NS) + type: integer + + :param: nPhi, number of points in phi (WE) + type: integer + + :param: nR, number of points in depth + type: integer + + :param: thetaSN (S, N) extensions of the model in degree + type: tuple + + :param: phiWE (W, E) extensions of the model in degree + type: tuple + + :param: Rbt (bot, top) extensions of the model in km + type: tuple + + :param: vbot, velocity at the bottom of the model + type: real + ''' + + def getRad(angle): + PI = np.pi + rad = angle / 180 * PI + return rad + + def getZmax(surface): + z = [] + for point in surface: + z.append(point[2]) + return max(z) + + R = 6371 + vmin = 0.34 + decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3) + outfile = open(filename, 'w') + + thetaS, thetaN = thetaSN + phiW, phiE = phiWE + rbot = Rbt[0] + R + rtop = Rbt[1] + R + + # need to determine the delta to add two cushion nodes around the min/max values + thetaDelta = abs(thetaN - thetaS) / (nTheta - 1) + phiDelta = abs(phiE - phiW) / (nPhi - 1) + rDelta = abs(rbot - rtop) / (nR - 1) + + # create a regular grid including +2 cushion nodes in every direction + thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes + phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes + rGrid = np.linspace(rbot - rDelta, rtop + rDelta, num = nR + 2) # +2 cushion nodes + + nTotal = len(rGrid) * len(thetaGrid) * len(phiGrid) + print "Total number of grid nodes: %s"%nTotal + + # write header for velocity grid file (in RADIANS) + outfile.writelines('%10s %10s \n' %(1, 1)) + outfile.writelines('%10s %10s %10s\n' %(nR + 2, nTheta + 2, nPhi + 2)) + outfile.writelines('%10s %10s %10s\n' %(rDelta, getRad(thetaDelta), getRad(phiDelta))) + outfile.writelines('%10s %10s %10s\n' %(rbot - rDelta, getRad(thetaS - thetaDelta), getRad(phiW - phiDelta))) + + surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method, filename = None) + zmax = getZmax(surface) + + print "\nGenerating velocity grid for FMTOMO. Output filename = %s, interpolation method = %s"%(filename, method) + print "nTheta = %s, nPhi = %s, nR = %s, thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt) + count = 0 + for radius in rGrid: + for theta in thetaGrid: + for phi in phiGrid: + xval = self._getDistance(phi) + yval = self._getDistance(theta) + for point in surface: + if point[0] == xval and point[1] == yval: + z = point[2] + if radius > (R + z + 1): + vel = 0.0 + # elif radius > (R + z - 15): ########### TESTING + # vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * 1.0 + vmin ########################## + else: + vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * vbot + vmin ########################## + count += 1 + outfile.writelines('%10s %10s\n'%(vel, decm)) + + progress = float(count) / float(nTotal) * 100 + self._update_progress(progress) + + outfile.close() + + def exportAll(self, filename = 'interpolated_receivers.out'): + recfile_out = open(filename, 'w') + count = 0 + for traceID in self.getReceiverCoordinates().keys(): + count += 1 + x, y, z = self.getReceiverCoordinates()[traceID] + recfile_out.writelines('%5s %15s %15s %15s\n' %(traceID, x, y, z)) + print "Exported coordinates for %s traces to file > %s" %(count, filename) + recfile_out.close() + + def plotArray2D(self, plot_topo = False, highlight_measured = False, annotations = True): + import matplotlib.pyplot as plt + plt.interactive(True) + plt.figure() + xmt, ymt, zmt = self.getMeasuredTopoLists() + xsc, ysc, zsc = self.getSourceLocsLists() + xmr, ymr, zmr = self.getMeasuredReceiverLists() + xrc, yrc, zrc = self.getReceiverLists() + + plt.plot(xrc, yrc, 'k.', markersize = 10, label = 'all receivers') + plt.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations') + + if plot_topo == True: + plt.plot(xmt, ymt, 'b', markersize = 10, label = 'measured topo points') + if highlight_measured == True: + plt.plot(xmr, ymr, 'ro', label = 'measured receivers') + + plt.xlabel('X [m]') + plt.ylabel('Y [m]') + plt.legend() + if annotations == True: + for traceID in self.getReceiverCoordinates().keys(): + plt.annotate(str(traceID), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small') + + def plotArray3D(self, ax = None): + import matplotlib.pyplot as plt + from mpl_toolkits.mplot3d import Axes3D + plt.interactive(True) + + if ax == None: + fig = plt.figure() + ax = plt.axes(projection = '3d') + + xmt, ymt, zmt = self.getMeasuredTopoLists() + xmr, ymr, zmr = self.getMeasuredReceiverLists() + xin, yin, zin = self.getReceiverLists() + + ax.plot(xmt, ymt, zmt, 'b*', markersize = 10, label = 'measured topo points') + ax.plot(xin, yin, zin, 'k.', markersize = 10, label = 'interpolated receivers') + ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers') + ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation') + ax.legend() + + return ax + + + def plotSurface3D(self, ax = None, step = 0.5, method = 'linear'): + from matplotlib import cm + import matplotlib.pyplot as plt + from mpl_toolkits.mplot3d import Axes3D + plt.interactive(True) + + if ax == None: + fig = plt.figure() + ax = plt.axes(projection = '3d') + + xmt, ymt, zmt = self.getMeasuredTopoLists() + xmr, ymr, zmr = self.getMeasuredReceiverLists() + + x = xmt + xmr + y = ymt + ymr + z = zmt + zmr + + xaxis = np.arange(min(x)+1, max(x), step) + yaxis = np.arange(min(y)+1, max(y), step) + + xgrid, ygrid = np.meshgrid(xaxis, yaxis) + + zgrid = griddata((x, y), z, (xgrid, ygrid), method = method) + + ax.plot_surface(xgrid, ygrid, zgrid, linewidth = 0, cmap = cm.jet, vmin = min(z), vmax = max(z)) + + ax.set_zlim(-(max(x) - min(x)/2),(max(x) - min(x)/2)) + ax.set_aspect('equal') + + ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation') + ax.legend() + + return ax + + def _update_progress(self, progress): + sys.stdout.write("%d%% done \r" % (progress) ) + sys.stdout.flush() + + def receivers2VTK(self, filename = 'receivers.vtk'): + ''' + Generates vtk files from all receivers of the SeisArray object. + ''' + outfile = open(filename, 'w') + traceIDs = [] + + for traceID in self.getReceiverCoordinates(): + traceIDs.append(traceID) + + nPoints = len(traceIDs) + + # write header + print("Writing header for VTK file...") + outfile.writelines('# vtk DataFile Version 3.1\n') + outfile.writelines('Receivers with traceIDs\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 traceID in traceIDs: + x = self._getXreceiver(traceID) + y = self._getYreceiver(traceID) + z = self._getZreceiver(traceID) + + 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 traceIDs int %d\n' %(1)) + outfile.writelines('LOOKUP_TABLE default\n') + + # write traceIDs + print("Writing traceIDs to VTK file...") + for traceID in traceIDs: + outfile.writelines('%10d\n' %traceID) + + outfile.close() + print("Wrote receiver grid for %d points to file: %s" %(nPoints, filename)) + return + + def sources2VTK(self, filename = 'sources.vtk'): + ''' + Generates vtk-files for all source locations in the SeisArray object. + ''' + outfile = open(filename, 'w') + shotnumbers = [] + + for shotnumber in self.getSourceLocations(): + shotnumbers.append(shotnumber) + + nPoints = len(shotnumbers) + + # write header + print("Writing header for VTK file...") + outfile.writelines('# vtk DataFile Version 3.1\n') + outfile.writelines('Shots with shotnumbers\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 shotnumber in shotnumbers: + x = self._getXshot(shotnumber) + y = self._getYshot(shotnumber) + z = self._getZshot(shotnumber) + + 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 shotnumbers int %d\n' %(1)) + outfile.writelines('LOOKUP_TABLE default\n') + + # write shotnumber + print("Writing shotnumbers to VTK file...") + for shotnumber in shotnumbers: + outfile.writelines('%10d\n' %shotnumber) + + outfile.close() + print("Wrote receiver grid for %d points to file: %s" %(nPoints, filename)) + return + + + def saveSeisArray(self, filename = 'seisArray.pickle'): + import cPickle + outfile = open(filename, 'wb') + + cPickle.dump(self, outfile, -1) + print('saved SeisArray to file %s'%(filename)) + + @staticmethod + def from_pickle(filename): + import cPickle + infile = open(filename, 'rb') + return cPickle.load(infile)