marcel/pylot
3
0
Fork 0
Code Issues 17 Pull Requests 1 Releases Wiki Activity

Merge branch 'develop'

Browse Source 
Conflicts:
	pylot/core/io/phases.py
This commit is contained in:
Sebastian Wehling-Benatelli 2016-05-18 13:02:51 +02:00
commit eb5f028a47
3 changed files with 348 additions and 271 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

122
pylot/core/active/activeSeismoPick.py
View File

@ -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')

482
pylot/core/active/fmtomoUtils.py
View File

@ -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.directories.append(directory)
def clearDIR(self, directory):
err = os.system('rm -r %s'%directory)
self.clearDir(directory)
self.directories.append(directory)
return
raise RuntimeError('Could not create directory: %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):
@ -208,6 +267,18 @@ class Tomo3d(object):
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()
srcfile = open(self.csl, 'r')
@ -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'])
with open(directory + '/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\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 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):
def mergeFrechet(self, directory):
print('Merging frechet.dat...')
with open(self.cwd + '/frechet.dat', 'w') as outfile:
frechfnout = os.path.join(directory, self.frechout)
with open(frechfnout, 'w') as outfile:
for procID in range(1, self.nproc + 1):
filename = self.getProcDir(procID) + '/frechet.dat'
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, directory):
self.mergeArrivals(directory)
self.mergeFrechet(directory)
self.mergeRays(directory)
def mergeOutput(self):
self.mergeArrivals()
self.mergeFrechet()
self.mergeRays()
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'):
'''

7
pylot/core/io/phases.py
View File

@ -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')