diff --git a/pylot/core/active/picking_script.py b/pylot/core/active/picking_script.py deleted file mode 100644 index 6a92b61b..00000000 --- a/pylot/core/active/picking_script.py +++ /dev/null @@ -1,107 +0,0 @@ -# -*- coding: utf-8 -*- -import sys -from obspy import read -from obspy import Stream -from obspy import Trace -from datetime import datetime -import numpy as np - -from pylot.core.active import surveyUtils -from pylot.core.active import seismicshot -import activeSeismoPick -reload(seismicshot) -reload(surveyUtils) - -##################################################################################### -# parameter definitions:############################################################# -#traceslist = list(np.arange(1, 49, 1)) # traces (1-48) -#shotlist = list(np.arange(302, 352, 1)) # shot-files(.dat) (Paffrath: 302-352) (Hauburg: 353-401) (arange+1!) -cutwindow = (0, 0.2) # cut out a part of the trace [seconds] -tmovwind = 0.3 # size of the moving window -windowsize = (5, 0) # windowsize for AIC picker (indices around HOS picks!!) -pickwindow = cutwindow # for local max and threshold picker: fraction of the seismogram used (0...1) TO BE DONE: depends on cutwindow!!!! -folm = 0.6 - -rockeskyll = False -if rockeskyll == True: - receiverfile = "Geophone_interpoliert_rockes" - sourcefile = "Schusspunkte_rockes" - obsdir = "/rscratch/minos22/marcel/flachseismik/rockeskyll_200615_270615/" - filename = 'survey_rockes.pickle' -else: - receiverfile = "Geophone_interpoliert_GZB" - sourcefile = "Schusspunkte_GZB" - obsdir = "/rscratch/minos22/marcel/flachseismik/GZB_26_06_15_01/" - filename = 'survey_GZB.pickle' - -# SNR -tsignal = 0.03 -tgap = 0.0007 -snrthreshold = 1 -###################################################################################### - -vmin = 333 -vmax = 5500 -distBinSize = 2 - -########################################### -############# Settings: ################### -thresholdpick=True -localmaxpick=False - -if thresholdpick == True: pickmethod = "threshold" -if localmaxpick == True: pickmethod = "localmax" - -HosAic = 'hos' # can be 'hos' or 'aic' -########################################### - -starttime = datetime.now() - -print '\n--------------------Starting Script at %s -------------------\n' %starttime.time() -if thresholdpick == True: print 'Using treshold pickmethod!\n' -elif localmaxpick == True: print 'Using local maximum pick method!\n' -if HosAic == 'aic': print 'picking with AIC\n' -if HosAic == 'hos': print 'picking with HOS\n' - -survey = activeSeismoPick.Survey(obsdir, sourcefile, receiverfile, True) -surveyUtils.setFittedSNR(survey.getShotDict()) -print '\nDone after %s seconds!\n------------------------------------------------------------------------------\n' % (datetime.now() - starttime).seconds - -count = 0; tpicksum = starttime - starttime - -for shot in survey.data.values(): - tstartpick = datetime.now(); count += 1 - for traceID in shot.getTraceIDlist(): - distance = shot.getDistance(traceID) # receive distance - - pickwin_used = pickwindow # use pickwindow set in the parameter section - # for higher distances use a linear vmin/vmax to cut out late/early regions with high noise - if distance > 5.: - pwleft = distance/vmax ################## TEST - pwright = distance/vmin - if pwright > cutwindow[1]: - pwright = cutwindow[1] - pickwin_used = (pwleft, pwright) - - shot.setPickwindow(traceID, pickwin_used) - shot.pickTraces(traceID, windowsize, folm, HosAic) # picker - #shot.setManualPicks(traceID, picklist) # set manual picks if given (yet used on 2D only) - - # ++ TEST: set and check SNR before adding to distance bin ############################ - shot.setSNR(traceID) - #if shot.getSNR(traceID)[0] < snrthreshold: - if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID): - shot.removePick(traceID) - # -- TEST: set and check SNR before adding to distance bin ############################ - - if shot.getPick(traceID) is not None: - shot.setEarllatepick(traceID) - - tpicksum += (datetime.now() - tstartpick); tpick = tpicksum/count - tremain = (tpick * (len(survey.getShotDict()) - count)) - tend = datetime.now() + tremain - print 'shot: %s, est. time to be finished is %s:%s:%s' % (shot.getShotname(), tend.hour, tend.minute, tend.second) - -survey.saveSurvey(filename) -print '\n--- Finished script ---' -print 'Elapsed time:', datetime.now()-starttime diff --git a/pylot/core/active/seismicArrayPreperation.py b/pylot/core/active/seismicArrayPreperation.py deleted file mode 100644 index 0bd3b1ad..00000000 --- a/pylot/core/active/seismicArrayPreperation.py +++ /dev/null @@ -1,666 +0,0 @@ -# -*- coding: utf-8 -*- -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)