diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index afcfd200..2d02496d 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -4,11 +4,10 @@ import numpy as np from pylot.core.active import seismicshot from pylot.core.active.surveyUtils import cleanUp - class Survey(object): def __init__(self, path, sourcefile, receiverfile, useDefaultParas=False): ''' - The Survey Class contains all shots [type: seismicshot] of a survey + The Survey Class contains all shots [class: Seismicshot] of a survey as well as the aquisition geometry and the topography. It contains methods to pick all traces of all shots. @@ -24,7 +23,7 @@ class Survey(object): self._generateSurvey() self._initiateFilenames() if useDefaultParas == True: - self.setParametersForShots() + self.setParametersForAllShots() self._removeAllEmptyTraces() self._updateShots() @@ -51,32 +50,36 @@ class Survey(object): print ("Total number of traces: %d \n" % self.countAllTraces()) def _removeAllEmptyTraces(self): - filename = 'removeEmptyTraces.out' + ''' + Removes traces of the dataset that are not found in the input receiver files. + ''' + logfile = 'removeEmptyTraces.out' count = 0 for shot in self.data.values(): removed = shot.removeEmptyTraces() if removed is not None: - if count == 0: outfile = open(filename, 'w') + if count == 0: outfile = open(logfile, 'w') count += 1 outfile.writelines('shot: %s, removed empty traces: %s\n' % (shot.getShotnumber(), removed)) print ("\nremoveEmptyTraces: Finished! Removed %d traces" % count) if count > 0: print ("See %s for more information " - "on removed traces." % (filename)) + "on removed traces." % (logfile)) outfile.close() def _updateShots(self): ''' - Removes traces that do not exist in the dataset for any reason. + Removes traces that do not exist in the dataset for any reason, + but were set in the input files. ''' - filename = 'updateShots.out' + logfile = 'updateShots.out' count = 0; countTraces = 0 for shot in self.data.values(): del_traceIDs = shot.updateTraceList() if len(del_traceIDs) > 0: - if count == 0: outfile = open(filename, 'w') + if count == 0: outfile = open(logfile, 'w') count += 1 countTraces += len(del_traceIDs) outfile.writelines("shot: %s, removed traceID(s) %s because " @@ -87,36 +90,37 @@ class Survey(object): "%d traces" % (count, countTraces)) if count > 0: print ("See %s for more information " - "on removed traces." % (filename)) + "on removed traces." % (logfile)) outfile.close() def setArtificialPick(self, traceID, pick): ''' Sets an artificial pick for a traceID of all shots in the survey object. - (This can be used to create a pick with t = 0 at the source origin) + (Commonly used to generate a travel time t = 0 at the source origin) ''' for shot in self.data.values(): shot.setPick(traceID, pick) - def setParametersForShots(self, cutwindow=(0, 0.2), tmovwind=0.3, tsignal=0.03, tgap=0.0007): + def setParametersForAllShots(self, cutwindow=(0, 0.2), tmovwind=0.3, tsignal=0.03, tgap=0.0007): if (cutwindow == (0, 0.2) and tmovwind == 0.3 and tsignal == 0.03 and tgap == 0.0007): print ("Warning: Standard values used for " "setParamters. This might not be clever.") - # CHANGE this later. Parameters only needed for survey, not for each shot. for shot in self.data.values(): shot.setCut(cutwindow) shot.setTmovwind(tmovwind) shot.setTsignal(tsignal) shot.setTgap(tgap) shot.setOrder(order=4) - print ("setParametersForShots: Parameters set to:\n" + print ("setParametersForAllShots: Parameters set to:\n" "cutwindow = %s, tMovingWindow = %f, tsignal = %f, tgap = %f" % (cutwindow, tmovwind, tsignal, tgap)) def setManualPicksFromFiles(self, directory='picks'): ''' Read manual picks from *.pck files in a directory. + Can be used for comparison of automatic and manual picks. + The * must be identical with the shotnumber. ''' for shot in self.data.values(): @@ -125,6 +129,7 @@ class Survey(object): def getDiffsFromManual(self): ''' Returns a dictionary with the differences between manual and automatic pick for all shots. + Key: Seismicshot [object] ''' diffs = {} for shot in self.data.values(): @@ -136,6 +141,10 @@ class Survey(object): return diffs def plotDiffs(self): + ''' + Creates a plot of all Picks colored by the + difference between automatic and manual pick. + ''' import matplotlib.pyplot as plt diffs = []; dists = []; @@ -163,8 +172,12 @@ class Survey(object): ax.set_xlabel('Distance [m]') ax.set_ylabel('Time [s]') ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes, horizontalalignment='center') + plt.legend() def plotHist(self, nbins=20, ax=None): + ''' + Plot a histogram of the difference between automatic and manual picks. + ''' import matplotlib.pyplot as plt plt.interactive(True) diffs = [] @@ -180,9 +193,21 @@ class Survey(object): plt.xlabel('Difference in time (auto - manual) [s]') return diffs - def pickAllShots(self, windowsize, HosAic='hos', vmin=333, vmax=5500, folm=0.6): + def pickAllShots(self, vmin=333, vmax=5500, folm=0.6, HosAic='hos', aicwindow=(10, 0)): ''' Automatically pick all traces of all shots of the survey. + + :param: vmin, vmax, minimum (maximum) permitted apparent velocity on direct path between src and rec + :type: real + + :param: folm, fraction of local maximum for HOS pick (0.6 = 60% of the highest maximum) + :type: real + + :param: HosAic, pick with hos only ('hos') or use AIC ('aic') + :type: string + + :param: aicwindow, window around the initial pick to search for local AIC min (samples) + :type: tuple ''' from datetime import datetime starttime = datetime.now() @@ -207,7 +232,7 @@ class Survey(object): pickwin_used = (pwleft, pwright) shot.setPickwindow(traceID, pickwin_used) - shot.pickTraces(traceID, windowsize, folm, HosAic) # picker + shot.pickTraces(traceID, folm, HosAic, aicwindow) # picker shot.setSNR(traceID) # if shot.getSNR(traceID)[0] < snrthreshold: @@ -231,6 +256,9 @@ class Survey(object): % (pickedtraces, ntraces, float(pickedtraces) / float(ntraces) * 100.)) def cleanBySPE(self, maxSPE): + ''' + Sets all picks as invalid if they exceed a certain value of the symmetric pick error. + ''' for shot in self.data.values(): for traceID in shot.getTraceIDlist(): if shot.getPickFlag(traceID) == 1: @@ -238,6 +266,9 @@ class Survey(object): shot.setPickFlag(traceID, 0) def plotSPE(self): + ''' + Plots the symmetric pick error sorted by magnitude. + ''' import matplotlib.pyplot as plt spe = [] for shot in self.data.values(): @@ -251,7 +282,7 @@ class Survey(object): def recover(self): ''' - Recovers all (accidently) removed picks. Still regards SNR threshold. + Recovers all manually removed picks. Still regards SNR threshold. ''' print('Recovering survey...') numpicks = 0 @@ -266,18 +297,28 @@ class Survey(object): print('Recovered %d picks' % numpicks) def setArtificialPick(self, traceID, pick): + ''' + Sets an artificial pick for a certain receiver (traceID) for all shots. + ''' for shot in self.data.values(): shot.setPick(traceID, pick) shot.setPickwindow(traceID, shot.getCut()) def countAllTraces(self): + ''' + Returns the number of traces in total. + ''' numtraces = 0 for shot in self.getShotlist(): - for rec in self.getReceiverlist(): ### shot.getReceiverlist etc. + for rec in self.getReceiverlist(): numtraces += 1 + return numtraces def getShotlist(self): + ''' + Returns a list of all shotnumbers contained in the set Sourcefile. + ''' filename = self.getPath() + self.getSourcefile() srcfile = open(filename, 'r') shotlist = [] @@ -288,6 +329,9 @@ class Survey(object): return shotlist def getReceiverlist(self): + ''' + Returns a list of all trace IDs contained in the set Receiverfile. + ''' filename = self.getPath() + self.getReceiverfile() recfile = open(filename, 'r') reclist = [] @@ -313,6 +357,12 @@ class Survey(object): return self._obsdir def getStats(self): + ''' + Generates and returns a dictionary containing statistical information + of the survey. + + Key: shotnumber + ''' info_dict = {} for shot in self.data.values(): pickedTraces = 0 @@ -334,11 +384,17 @@ class Survey(object): return info_dict def getShotForShotnumber(self, shotnumber): + ''' + Returns Seismicshot [object] of a certain shotnumber if possible. + ''' for shot in self.data.values(): if shot.getShotnumber() == shotnumber: return shot def exportFMTOMO(self, directory='FMTOMO_export', sourcefile='input_sf.in', ttFileExtension='.tt'): + ''' + Exports all picks into a directory as travel time files readable by FMTOMO obsdata. + ''' def getAngle(distance): PI = np.pi R = 6371. @@ -354,13 +410,13 @@ class Survey(object): srcfile.writelines('%10s\n' % len(self.data)) # number of sources for shotnumber in self.getShotlist(): shot = self.getShotForShotnumber(shotnumber) - ttfilename = str(shotnumber) + ttFileExtension + ttfilename = str(shotnumber) + ttFileExtension # filename of travel time file for this shot (x, y, z) = shot.getSrcLoc() # getSrcLoc returns (x, y, z) - srcfile.writelines('%10s %10s %10s\n' % (getAngle(y), getAngle(x), (-1) * z)) # lat, lon, depth + srcfile.writelines('%10s %10s %10s\n' % (getAngle(y), getAngle(x), (-1) * z)) # transform to lat, lon, depth LatAll.append(getAngle(y)); LonAll.append(getAngle(x)); DepthAll.append((-1) * z) - srcfile.writelines('%10s\n' % 1) # + srcfile.writelines('%10s\n' % 1) srcfile.writelines('%10s %10s %10s\n' % (1, 1, ttfilename)) ttfile = open(directory + '/' + ttfilename, 'w') traceIDlist = shot.getTraceIDlist() @@ -393,6 +449,9 @@ class Survey(object): print(msg) def countPickedTraces(self, shot): + ''' + Counts all picked traces of a shot (type Seismicshot). + ''' count = 0 for traceID in shot.getTraceIDlist(): if shot.getPickFlag(traceID) is not 0: @@ -400,6 +459,9 @@ class Survey(object): return count def countAllPickedTraces(self): + ''' + Counts all picked traces of the survey. + ''' count = 0 for shot in self.data.values(): for traceID in shot.getTraceIDlist(): @@ -422,20 +484,9 @@ class Survey(object): figPerSubplot = columns * rows index = 1 - # shotnames = [] - # shotnumbers = [] - - # for shot in self.data.values(): - # shotnames.append(shot.getShotname()) - # shotnumbers.append(shot.getShotnumber()) - - # shotnumbers = [shotnumbers for (shotnumbers, shotnames) in sorted(zip(shotnumbers, shotnames))] for shotnumber in self.getShotlist(): if index <= figPerSubplot: - # ax = fig.add_subplot(3,3,i, projection = '3d', title = 'shot:' - # +str(shot_dict[shotnumber].getShotnumber()), xlabel = 'X', ylabel = 'Y', zlabel = 'traveltime') - # shot_dict[shotnumber].plot3dttc(ax = ax, plotpicks = True) ax = fig.add_subplot(rows, columns, index) if mode == '3d': self.getShot(shotnumber).matshow(ax=ax, colorbar=False, annotations=True, legend=False) @@ -516,6 +567,9 @@ class Survey(object): return ax def createPlot(self, dist, pick, inkByVal, label, ax=None, cbar=None): + ''' + Used by plotAllPicks. + ''' import matplotlib.pyplot as plt plt.interactive(True) cm = plt.cm.jet @@ -547,6 +601,10 @@ class Survey(object): sys.stdout.flush() def saveSurvey(self, filename='survey.pickle'): + ''' + Save Survey object to a file. + Can be loaded by using Survey.from_pickle(filename). + ''' import cPickle cleanUp(self) outfile = open(filename, 'wb') diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index 7dfe928b..f70b0b43 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -6,15 +6,28 @@ import datetime import numpy as np class Tomo3d(object): - def __init__(self, nproc, iterations): + def __init__(self, nproc, iterations, citer = 0, overwrite = False): + ''' + Class build from FMTOMO script tomo3d. Can be used to run several instances of FMM code in parallel. + + :param: nproc, number of parallel processes + :type: integer + + :param: iterations, number of iterations + :type: integer + + :param: citer, current iteration (default = 0: start new model) + :type: integer + ''' self.setCWD() self.defParas() self.nproc = nproc self.iter = iterations # number of iterations - self.citer = 0 # current iteration + self.citer = citer # current iteration self.sources = self.readSrcFile() self.traces = self.readTraces() self.directories = [] + self.overwrite = overwrite def defParas(self): self.defFMMParas() @@ -29,6 +42,7 @@ class Tomo3d(object): self.pg = 'propgrid.in' self.rec = 'receivers.in' self.frech = 'frechet.in' + self.frechout = 'frechet.dat' self.ot = 'otimes.dat' self.ttim = 'arrivals.dat' self.mode = 'mode_set.in' @@ -53,8 +67,8 @@ class Tomo3d(object): self.resout = 'residuals.dat' def copyRef(self): - # Attention: Copies reference grids to used grids (e.g. sourcesref.in to sources.in) - os.system('cp %s %s'%(self. ivg, self.cvg)) + # Copies reference grids to used grids (e.g. sourcesref.in to sources.in) + os.system('cp %s %s'%(self.ivg, self.cvg)) os.system('cp %s %s'%(self.iig, self.cig)) os.system('cp %s %s'%(self.isl, self.csl)) @@ -65,6 +79,30 @@ class Tomo3d(object): def runFrech(self): os.system(self.frechgen) + def runTOMO3D(self): + starttime = datetime.datetime.now() + print('Starting TOMO3D on %s parallel processes for %s iteration(s).' + %(self.nproc, self.iter)) + if self.citer == 0: + self.makeInvIterDir() + self.startForward(self.cInvIterDir) + self.raiseIter() + + while self.citer <= self.iter: + self.makeInvIterDir() + self.startInversion() + self.saveVgrid() + self.startForward(self.cInvIterDir) + self.raiseIter() + + if self.citer > self.iter: + self.removeDirectories() + self.unlink(os.path.join(self.cwd, self.frechout)) + self.unlink(os.path.join(self.cwd, self.ttim)) + + tdelta = datetime.datetime.now() - starttime + print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta)) + def runFmm(self, directory, logfile, processes): os.chdir(directory) fout = open(logfile, 'w') @@ -73,6 +111,47 @@ class Tomo3d(object): os.chdir(self.cwd) return processes + def startForward(self, logdir): + self._printLine() + print('Starting forward simulation for iteration %s.'%(self.citer)) + + if self.citer == 0: + self.copyRef() + self.runFrech() + self.makeDirectories() + + starttime = datetime.datetime.now() + processes = [] + + for procID in range(1, self.nproc + 1): + directory = self.getProcDir(procID) + logfn = 'fm3dlog_' + str(procID) + '.out' + log_out = os.path.join(logdir, logfn) + + self.writeSrcFile(procID, directory) + self.writeTracesFile(procID, directory) + os.system('cp {cvg} {cig} {mode} {pg} {frechin} {dest}' + .format(cvg=self.cvg, cig=self.cig, frechin=self.frech, + mode=self.mode, pg=self.pg, dest=directory)) + processes = self.runFmm(directory, log_out, processes) + + for p in processes: + p.wait() + + self.mergeOutput(self.cInvIterDir) + self.clearDirectories() + self.copyArrivals() + if self.citer == 0: + self.copyArrivals(self.rtrav) + + self.calcRes() + tdelta = datetime.datetime.now() - starttime + print('Finished Forward calculation after %s'%tdelta) + + def startInversion(self): + print('Calling %s...'%self.inv) + os.system(self.inv) + def calcRes(self): resout = os.path.join(self.cwd, self.resout) if self.citer == 0: @@ -85,111 +164,91 @@ class Tomo3d(object): RMS, var, chi2 = residuals[-1].split() print('Residuals: RMS = %s, var = %s, Chi^2 = %s.'%(RMS, var, chi2)) - def runTOMO3D(self): - starttime = datetime.datetime.now() - print('Starting TOMO3D on %s parallel processes for %s iteration(s).' - %(self.nproc, self.iter)) - if self.citer == 0: - self.startForward() - self.raiseIter() - - while self.citer <= self.iter: - self.startInversion() - self.startForward() - self.raiseIter() - - tdelta = datetime.datetime.now() - starttime - print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta)) - - def _printLine(self): - print('----------------------------------------') - def raiseIter(self): self.citer +=1 self._printLine() + invfile = open(self.cwd + '/inviter.in', 'w') + invfile.write('%s'%self.citer) + invfile.close() - def startInversion(self): - print('Calling %s...'%self.inv) - os.system(self.inv) - - def startForward(self): - self._printLine() - print('Starting forward simulation for iteration %s.'%(self.citer)) - self.makeDirectories() - starttime = datetime.datetime.now() - processes = [] - - if self.citer == 0: - self.copyRef() - self.runFrech() - - for procID in range(1, self.nproc + 1): - directory = self.getProcDir(procID) - log_out = self.cwd + '/fm3dlog_' + str(procID) + '.out' - - self.writeSrcFile(procID, directory) - self.writeTracesFile(procID, directory) - os.system('cp {cvg} {cig} {mode} {pg} {frechout} {dest}' - .format(cvg=self.cvg, cig=self.cig, frechout=self.frech, - mode=self.mode, pg=self.pg, dest=directory)) - processes = self.runFmm(directory, log_out, processes) - - for p in processes: - p.wait() - - self.mergeOutput() - self.clearDirectories() - self.copyArrivals() - if self.citer == 0: - self.copyArrivals(self.rtrav) - - self.calcRes() - tdelta = datetime.datetime.now() - starttime - print('Finished Forward calculation after %s'%tdelta) - - def copyArrivals(self, target = None): - if target == None: - target = self.mtrav - os.system('cp %s %s'%(self.ttim, target)) - - def makeDIR(self, directory): + def makeDir(self, directory): err = os.system('mkdir %s'%directory) + if err is 0: + self.directories.append(directory) + return if err is 256: - response = raw_input('Warning: Directory already existing. Continue (y/n)?\n') - if response == 'y': + if self.overwrite == True: print('Overwriting existing files.') - else: - sys.exit('Aborted') + self.clearDir(directory) + self.directories.append(directory) + return + raise RuntimeError('Could not create directory: %s'%directory) - self.directories.append(directory) - - def clearDIR(self, directory): - err = os.system('rm -r %s'%directory) - def makeDirectories(self): for procID in range(1, self.nproc + 1): directory = self.getProcDir(procID) - self.makeDIR(directory) + self.makeDir(directory) + def makeInvIterDir(self): + invIterDir = self.cwd + '/it_%s'%(self.citer) + err = os.system('mkdir %s'%invIterDir) + if err is 256: + if self.overwrite: + self.clearDir(invIterDir) + elif err is not 0: + raise RuntimeError('Could not create directory: %s'%invIterDir) + self.cInvIterDir = invIterDir + + def clearDir(self, directory): + print('Wiping directory %s...'%directory) + for filename in os.listdir(directory): + filenp = os.path.join(directory, filename) + os.remove(filenp) + def clearDirectories(self): for directory in self.directories: - self.clearDIR(directory) + self.clearDir(directory) + + def rmDir(self, directory): + print('Removing directory %s...'%directory) + return os.rmdir(directory) + + def removeDirectories(self): + for directory in self.directories: + self.rmDir(directory) self.directories = [] - def readNsrc(self): - srcfile = open(self.csl, 'r') - nsrc = int(srcfile.readline()) - srcfile.close() - return nsrc + def getProcDir(self, procID): + return os.path.join(self.cwd, self.folder) + str(procID) - def readNtraces(self): - recfile = open(self.rec, 'r') - nrec = int(recfile.readline()) - recfile.close() - return nrec + def getTraceIDs4Sources(self, sourceIDs): + traceIDs = [] + for traceID in self.traces.keys(): + if self.traces[traceID]['source'] in sourceIDs: + traceIDs.append(traceID) + return traceIDs + + def getTraceIDs4Source(self, sourceID): + traceIDs = [] + for traceID in self.traces.keys(): + if self.traces[traceID]['source'] == sourceID: + traceIDs.append(traceID) + return traceIDs + + def copyArrivals(self, target = None): + if target == None: + target = os.path.join(self.cwd, self.mtrav) + os.system('cp %s %s'%(os.path.join( + self.cInvIterDir, self.ttim), target)) + + def saveVgrid(self): + vgpath = os.path.join(self.cwd, 'vgrids.in') + os.system('cp %s %s'%(vgpath, self.cInvIterDir)) def calcSrcPerKernel(self): nsrc = self.readNsrc() + if self.nproc > nsrc: + raise ValueError('Number of spawned processes higher than number of sources') return nsrc/self.nproc, nsrc%self.nproc def srcIDs4Kernel(self, procID): @@ -207,6 +266,18 @@ class Tomo3d(object): elif proc > remain: start = (srcPK + 1) * remain + srcPK * (proc - remain) + 1 return range(start, start + srcPK) + + def readNsrc(self): + srcfile = open(self.csl, 'r') + nsrc = int(srcfile.readline()) + srcfile.close() + return nsrc + + def readNtraces(self): + recfile = open(self.rec, 'r') + nrec = int(recfile.readline()) + recfile.close() + return nrec def readSrcFile(self): nsrc = self.readNsrc() @@ -258,6 +329,62 @@ class Tomo3d(object): return traces + def readArrivals(self, procID): + directory = self.getProcDir(procID) + arrfile = open(directory + '/arrivals.dat', 'r') + sourceIDs = self.srcIDs4Kernel(procID) + + arrivals = [] + for sourceID in sourceIDs: + traceIDs = self.getTraceIDs4Source(sourceID) + for traceID in traceIDs: + line = arrfile.readline().split() + if line != []: + # recID and srcID for the individual processor will not be needed + recID_proc, srcID_proc, ray, normal, arrtime, diff, head = line + arrivals.append([traceID, sourceID, ray, normal, arrtime, diff, head]) + + return arrivals + + def readRays(self, procID): + directory = self.getProcDir(procID) + raysfile = open(directory + '/rays.dat', 'r') + sourceIDs = self.srcIDs4Kernel(procID) + + rays = {} + for sourceID in sourceIDs: + traceIDs = self.getTraceIDs4Source(sourceID) + for traceID in traceIDs: + line1 = raysfile.readline().split() + if line1 != []: + # recID and srcID for the individual processor will not be needed + recID_proc, srcID_proc, ray, normal, nsec = line1 + raysecs = {} + + for sec in range(int(nsec)): + line2 = raysfile.readline().split() + npoints, region, diff, head = line2 + raypoints = [] + + for j in range(int(npoints)): + raypoints.append(raysfile.readline() + '\n') + + raysecs[sec] = {'npoints': npoints, + 'region': region, + 'diff': diff, + 'head': head, + 'raypoints': raypoints + } + + + rays[traceID] = {'sourceID': sourceID, + 'raypath': ray, + 'normal': normal, + 'nsec': nsec, + 'raysections': raysecs + } + return rays + def writeSrcFile(self, procID, directory): srcfile = open('%s/sources.in'%directory, 'w') sourceIDs = self.srcIDs4Kernel(procID) @@ -289,154 +416,44 @@ class Tomo3d(object): recfile.write('%s\n'%source) recfile.write('%s\n'%trace['path']) - def getTraceIDs4Sources(self, sourceIDs): - traceIDs = [] - for traceID in self.traces.keys(): - if self.traces[traceID]['source'] in sourceIDs: - traceIDs.append(traceID) - return traceIDs - - def getTraceIDs4Source(self, sourceID): - traceIDs = [] - for traceID in self.traces.keys(): - if self.traces[traceID]['source'] == sourceID: - traceIDs.append(traceID) - return traceIDs - - def readArrivals(self, procID): - directory = self.getProcDir(procID) - arrfile = open(directory + '/arrivals.dat', 'r') - sourceIDs = self.srcIDs4Kernel(procID) - - arrivals = [] - for sourceID in sourceIDs: - traceIDs = self.getTraceIDs4Source(sourceID) - for traceID in traceIDs: - line = arrfile.readline().split() - if line != []: - # recID and srcID for the individual processor will not be needed - recID_proc, srcID_proc, ray, normal, arrtime, diff, head = line - arrivals.append([traceID, sourceID, ray, normal, arrtime, diff, head]) - - return arrivals - - def mergeArrivals(self): - arrivalsOut = open(os.path.join(self.cwd, self.ttim), 'w') + def mergeArrivals(self, directory): + arrfn = os.path.join(directory, self.ttim) + arrivalsOut = open(arrfn, 'w') print('Merging arrivals.dat...') for procID in range(1, self.nproc + 1): arrivals = self.readArrivals(procID) for line in arrivals: arrivalsOut.write('%6s %6s %6s %6s %15s %5s %5s\n'%tuple(line)) - # def mergeFrechet(self): - # print('Merging frechet.dat...') - # frechetOut = open(self.cwd + '/frechet.dat', 'w') - # for procID in range(1, self.nproc + 1): - # frechet = self.readFrechet(procID) - # traceIDs = frechet.keys() - # traceIDs.sort() - # for traceID in traceIDs: - # frech = frechet[traceID] - # frechetOut.write('%6s %6s %6s %6s %6s\n'% - # (traceID, - # frech['sourceID'], - # frech['raypath'], - # frech['normal'], - # frech['NPDEV'])) - # for pdev in frech['PDEV']: - # frechetOut.writelines(pdev) + os.system('ln -fs %s %s'%(arrfn, os.path.join(self.cwd, self.ttim))) - def readRays(self, procID): - directory = self.getProcDir(procID) - raysfile = open(directory + '/rays.dat', 'r') - sourceIDs = self.srcIDs4Kernel(procID) - - rays = {} - for sourceID in sourceIDs: - traceIDs = self.getTraceIDs4Source(sourceID) - for traceID in traceIDs: - line1 = raysfile.readline().split() - if line1 != []: - # recID and srcID for the individual processor will not be needed - recID_proc, srcID_proc, ray, normal, nsec = line1 - raysecs = {} - - for sec in range(int(nsec)): - line2 = raysfile.readline.split() - npoints, region, diff, head = line2 - raypoints = [] - - for j in range(int(npoints)): - raypoints.append(raysfile.readline() + '\n') - - raysecs[sec] = {'npoints': npoints, - 'region': region, - 'diff': diff, - 'head': head, - 'raypoints': raypoints - } - - - rays[traceID] = {'sourceID': sourceID, - 'raypath': ray, - 'normal': normal, - 'nsec': nsec, - 'raysections': raysecs - } - return rays - - def mergeRays(self): + def mergeRays(self, directory): print('Merging rays.dat...') - with open(self.cwd + 'rays.dat', 'w') as outfile: - for procID in range(1, self.nproc + 1): - rays = self.readRays(procID) - for traceID in rays: - ray = rays[traceID] - outfile.write('%6s %6s %6s %6s %6s'%(traceID, - ray['sourceID'], - ray['raypath'], - ray['normal'], - ray['nsec'])) - for sec in range(int(ray['nsec'])): - raysec = ray['raysections'][sec] - outfile.write('%6s %6s %6s %6s'%(raysec['npoints'], - raysec['region'], - raysec['diff'], - raysec['head'])) - outfile.writelines(raysec['raypoints']) - - - # def readFrechet(self, procID): - # directory = self.getProcDir(procID) - # frechfile = open(directory + '/frechet.dat', 'r') - # sourceIDs = self.srcIDs4Kernel(procID) - - # frechet = {} - # for sourceID in sourceIDs: - # traceIDs = self.getTraceIDs4Source(sourceID) - # for traceID in traceIDs: - # line = frechfile.readline().split() - # if line != []: - # # recID and srcID for the individual processor will not be needed - # PDEV = [] - # recID_proc, srcID_proc, ray, normal, NPDEV = line - # for i in range(int(NPDEV)): - # PDEV.append(frechfile.readline()) - - # frechet[traceID] = {'sourceID': sourceID, - # 'raypath': ray, - # 'normal': normal, - # 'NPDEV': NPDEV, - # 'PDEV': PDEV - # } - # return frechet - - def mergeFrechet(self): - print('Merging frechet.dat...') - - with open(self.cwd + '/frechet.dat', 'w') as outfile: + with open(directory + '/rays.dat', 'w') as outfile: for procID in range(1, self.nproc + 1): - filename = self.getProcDir(procID) + '/frechet.dat' + rays = self.readRays(procID) + for traceID in rays: + ray = rays[traceID] + outfile.write('%6s %6s %6s %6s %6s\n'%(traceID, + ray['sourceID'], + ray['raypath'], + ray['normal'], + ray['nsec'])) + for sec in range(int(ray['nsec'])): + raysec = ray['raysections'][sec] + outfile.write('%6s %6s %6s %6s\n'%(raysec['npoints'], + raysec['region'], + raysec['diff'], + raysec['head'])) + outfile.writelines(raysec['raypoints']) + + def mergeFrechet(self, directory): + print('Merging frechet.dat...') + + frechfnout = os.path.join(directory, self.frechout) + with open(frechfnout, 'w') as outfile: + for procID in range(1, self.nproc + 1): + filename = os.path.join(self.getProcDir(procID), self.frechout) with open(filename) as infile: for sourceID in self.srcIDs4Kernel(procID): for traceID in self.getTraceIDs4Source(sourceID): @@ -445,15 +462,18 @@ class Tomo3d(object): for index in range(int(NPDEV)): outfile.write(infile.readline()) + os.system('ln -fs %s %s'%(frechfnout, os.path.join(self.cwd, self.frechout))) - def getProcDir(self, procID): - return os.path.join(self.cwd, self.folder) + str(procID) - - def mergeOutput(self): - self.mergeArrivals() - self.mergeFrechet() - self.mergeRays() + def mergeOutput(self, directory): + self.mergeArrivals(directory) + self.mergeFrechet(directory) + self.mergeRays(directory) + def unlink(self, filepath): + return os.system('unlink %s'%filepath) + + def _printLine(self): + print('----------------------------------------') def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', inputfileref='vgridsref.in'): ''' diff --git a/pylot/core/io/phases.py b/pylot/core/io/phases.py index b5770385..bf2ba499 100644 --- a/pylot/core/io/phases.py +++ b/pylot/core/io/phases.py @@ -199,7 +199,7 @@ def picksdict_from_picks(evt): :param evt: Event object contain all available information :type evt: `~obspy.core.event.Event` :return: pick dictionary - ''' + """ picks = {} for pick in evt.picks: phase = {} @@ -345,7 +345,7 @@ def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbos def writephases(arrivals, fformat, filename): - ''' + """ Function of methods to write phases to the following standard file formats used for locating earthquakes: @@ -363,7 +363,7 @@ def writephases(arrivals, fformat, filename): :param: filename, full path and name of phase file :type: string - ''' + """ if fformat == 'NLLoc': print ("Writing phases to %s for NLLoc" % filename) @@ -425,7 +425,6 @@ def writephases(arrivals, fformat, filename): sweight)) fid.close() - elif fformat == 'HYPO71': print ("Writing phases to %s for HYPO71" % filename) fid = open("%s" % filename, 'w')