From 46c152b7a139a1cf1e1e81bfb3ba5454d7d289b8 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Wed, 11 May 2016 12:02:09 +0200 Subject: [PATCH 1/3] improved speed on writing to file (write instead of writelines) --- pylot/core/active/fmtomoUtils.py | 438 +++++++++++++------ pylot/core/active/seismicArrayPreparation.py | 137 +++--- 2 files changed, 375 insertions(+), 200 deletions(-) diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index 92ef7e5d..7dfe928b 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -6,32 +6,152 @@ import datetime import numpy as np class Tomo3d(object): - def __init__(self, nproc): + def __init__(self, nproc, iterations): + self.setCWD() self.defParas() self.nproc = nproc + self.iter = iterations # number of iterations + self.citer = 0 # current iteration self.sources = self.readSrcFile() self.traces = self.readTraces() - + self.directories = [] + def defParas(self): - self.fmm = 'fm3d' + self.defFMMParas() + self.defInvParas() + + def defFMMParas(self): + self.fmm = '{0}/fm3d'.format(self.cwd) + self.frechgen = '{0}/frechgen'.format(self.cwd) self.cvg = 'vgrids.in' self.cig = 'interfaces.in' self.csl = 'sources.in' self.pg = 'propgrid.in' self.rec = 'receivers.in' - #self.ot = 'otimes.dat' self.frech = 'frechet.in' + self.ot = 'otimes.dat' + self.ttim = 'arrivals.dat' self.mode = 'mode_set.in' - self.cwd = subprocess.check_output(['pwd'])[0:-1] + '/' - self.folder = 'test_' + self.folder = '.proc_' + + def defInvParas(self): + # Name of program for performing inversion + self.inv = '{0}/invert3d'.format(self.cwd) + # Name of file containing current model traveltimes + self.mtrav = 'mtimes.dat' + # Name of file containing reference model traveltimes + self.rtrav = 'rtimes.dat' + # Name of file containing initial velocity grid + self.ivg = 'vgridsref.in' + # Name of file containing initial interface grid + self.iig = 'interfacesref.in' + # Name of file containing initial source locations + self.isl = 'sourcesref.in' + # Name of program for calculating traveltime residuals + self.resid = '{0}/residuals'.format(self.cwd) + # Name of output file for calculating traveltime residuals + 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)) + os.system('cp %s %s'%(self.iig, self.cig)) + os.system('cp %s %s'%(self.isl, self.csl)) + + def setCWD(self): + self.cwd = subprocess.check_output(['pwd'])[0:-1] + print('Working directory is pwd: %s'%self.cwd) + + def runFrech(self): + os.system(self.frechgen) def runFmm(self, directory, logfile, processes): os.chdir(directory) - processes.append(subprocess.Popen('fm3d', stdout = None)) -# os.system('%s > %s &'%(self.fmm, logfile)) + fout = open(logfile, 'w') + processes.append(subprocess.Popen(self.fmm, stdout = fout)) + fout.close() os.chdir(self.cwd) return processes + def calcRes(self): + resout = os.path.join(self.cwd, self.resout) + if self.citer == 0: + os.system('%s > %s'%(self.resid, resout)) + else: + os.system('%s >> %s'%(self.resid, resout)) + + with open(resout, 'r') as infile: + residuals = infile.readlines() + 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() + + 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): err = os.system('mkdir %s'%directory) if err is 256: @@ -40,6 +160,21 @@ class Tomo3d(object): print('Overwriting existing files.') else: sys.exit('Aborted') + + 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) + + def clearDirectories(self): + for directory in self.directories: + self.clearDIR(directory) + self.directories = [] def readNsrc(self): srcfile = open(self.csl, 'r') @@ -73,29 +208,6 @@ class Tomo3d(object): start = (srcPK + 1) * remain + srcPK * (proc - remain) + 1 return range(start, start + srcPK) - def startTomo(self): - starttime = datetime.datetime.now() - processes = [] - for procID in range(1, self.nproc + 1): - directory = self.getProcDir(procID) - log_out = self.cwd + 'fm3dlog_' + str(procID) + '.out' - - self.makeDIR(directory) # Problem bei Iteration - self.writeSrcFile(procID, directory) - self.writeTracesFile(procID, directory) - os.system('cp %s %s %s %s %s %s %s' - %(self.cvg, self.cig, self.frech, - self.fmm, self.mode, self.pg, directory)) - processes = self.runFmm(directory, log_out, processes) - - for p in processes: - p.wait() - - self.mergeOutput() - - tdelta = datetime.datetime.now() - starttime - print('Finished after %s'%tdelta) - def readSrcFile(self): nsrc = self.readNsrc() srcfile = open(self.csl, 'r') @@ -150,32 +262,32 @@ class Tomo3d(object): srcfile = open('%s/sources.in'%directory, 'w') sourceIDs = self.srcIDs4Kernel(procID) - srcfile.writelines('%s\n'%len(sourceIDs)) + srcfile.write('%s\n'%len(sourceIDs)) for sourceID in sourceIDs: source = self.sources[sourceID] coords = source['coords'] interactions = source['interactions'] - srcfile.writelines('%s\n'%source['teleflag']) - srcfile.writelines('%s %s %s\n'%(float(coords[0]), float(coords[1]), float(coords[2]))) - srcfile.writelines('%s\n'%source['numpaths']) - srcfile.writelines('%s\n'%source['steps']) - srcfile.writelines('%s %s\n'%(int(interactions[0]), int(interactions[1]))) - srcfile.writelines('%s\n'%source['veltype']) + srcfile.write('%s\n'%source['teleflag']) + srcfile.write('%s %s %s\n'%(float(coords[0]), float(coords[1]), float(coords[2]))) + srcfile.write('%s\n'%source['numpaths']) + srcfile.write('%s\n'%source['steps']) + srcfile.write('%s %s\n'%(int(interactions[0]), int(interactions[1]))) + srcfile.write('%s\n'%source['veltype']) def writeTracesFile(self, procID, directory): recfile = open('%s/receivers.in'%directory, 'w') sourceIDs = self.srcIDs4Kernel(procID) traceIDs = self.getTraceIDs4Sources(sourceIDs) - recfile.writelines('%s\n'%len(traceIDs)) + recfile.write('%s\n'%len(traceIDs)) for traceID in traceIDs: trace = self.traces[traceID] coords = trace['coords'] source = int(trace['source']) - sourceIDs[0] + 1 - recfile.writelines('%s %s %s\n'%(float(coords[0]), float(coords[1]), float(coords[2]))) - recfile.writelines('%s\n'%trace['paths']) - recfile.writelines('%s\n'%source) - recfile.writelines('%s\n'%trace['path']) + recfile.write('%s %s %s\n'%(float(coords[0]), float(coords[1]), float(coords[2]))) + recfile.write('%s\n'%trace['paths']) + recfile.write('%s\n'%source) + recfile.write('%s\n'%trace['path']) def getTraceIDs4Sources(self, sourceIDs): traceIDs = [] @@ -208,79 +320,139 @@ class Tomo3d(object): return arrivals - 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 mergeArrivals(self): - arrivalsOut = open(self.cwd + '/arrivals.dat', 'w') + arrivalsOut = open(os.path.join(self.cwd, self.ttim), 'w') print('Merging arrivals.dat...') for procID in range(1, self.nproc + 1): arrivals = self.readArrivals(procID) for line in arrivals: - arrivalsOut.writelines('%6s %6s %6s %6s %15s %5s %5s\n'%tuple(line)) + 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() + # 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) + + 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: - frech = frechet[traceID] - frechetOut.writelines('%6s %6s %6s %6s %6s\n'% - (traceID, - frech['sourceID'], - frech['raypath'], - frech['normal'], - frech['NPDEV'])) - for pdev in frech['PDEV']: - frechetOut.writelines(pdev) + 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): print('Merging rays.dat...') - filenames = [] - for procID in range(1, self.nproc + 1): - directory = self.getProcDir(procID) - filenames.append(directory + '/rays.dat') - with open(self.cwd + 'rays.dat', 'w') as outfile: - for fname in filenames: - with open(fname) as infile: - for line in infile: - outfile.write(line) + 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: + for procID in range(1, self.nproc + 1): + filename = self.getProcDir(procID) + '/frechet.dat' + with open(filename) as infile: + for sourceID in self.srcIDs4Kernel(procID): + for traceID in self.getTraceIDs4Source(sourceID): + recID_proc, srcID_proc, ray, normal, NPDEV = infile.readline().split() + outfile.write('%6s %6s %6s %6s %6s\n'%(traceID, sourceID, ray, normal, NPDEV)) + for index in range(int(NPDEV)): + outfile.write(infile.readline()) + def getProcDir(self, procID): - return self.cwd + self.folder + str(procID) + return os.path.join(self.cwd, self.folder) + str(procID) def mergeOutput(self): self.mergeArrivals() self.mergeFrechet() self.mergeRays() - - def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', inputfileref='vgridsref.in'): @@ -314,23 +486,23 @@ def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', i # write header print("Writing header for VTK file...") - outfile.writelines('# vtk DataFile Version 3.1\n') - outfile.writelines('Velocity on FMTOMO vgrids.in points\n') - outfile.writelines('ASCII\n') - outfile.writelines('DATASET STRUCTURED_POINTS\n') + outfile.write('# vtk DataFile Version 3.1\n') + outfile.write('Velocity on FMTOMO vgrids.in points\n') + outfile.write('ASCII\n') + outfile.write('DATASET STRUCTURED_POINTS\n') - outfile.writelines('DIMENSIONS %d %d %d\n' % (nX, nY, nZ)) - outfile.writelines('ORIGIN %f %f %f\n' % (sX, sY, sZ)) - outfile.writelines('SPACING %f %f %f\n' % (dX, dY, dZ)) + outfile.write('DIMENSIONS %d %d %d\n' % (nX, nY, nZ)) + outfile.write('ORIGIN %f %f %f\n' % (sX, sY, sZ)) + outfile.write('SPACING %f %f %f\n' % (dX, dY, dZ)) - outfile.writelines('POINT_DATA %15d\n' % (nPoints)) + outfile.write('POINT_DATA %15d\n' % (nPoints)) if absOrRel == 'abs': - outfile.writelines('SCALARS velocity float %d\n' %(1)) + outfile.write('SCALARS velocity float %d\n' %(1)) if absOrRel == 'relDepth': - outfile.writelines('SCALARS velocity2depthMean float %d\n' %(1)) + outfile.write('SCALARS velocity2depthMean float %d\n' %(1)) elif absOrRel == 'rel': - outfile.writelines('SCALARS velChangePercent float %d\n' % (1)) - outfile.writelines('LOOKUP_TABLE default\n') + outfile.write('SCALARS velChangePercent float %d\n' % (1)) + outfile.write('LOOKUP_TABLE default\n') pointsPerR = nTheta * nPhi @@ -338,7 +510,7 @@ def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', i if absOrRel == 'abs': print("Writing velocity values to VTK file...") for velocity in vel: - outfile.writelines('%10f\n' %velocity) + outfile.write('%10f\n' %velocity) elif absOrRel == 'relDepth': print("Writing velocity values to VTK file relative to mean of each depth...") index = 0; count = 0 @@ -350,7 +522,7 @@ def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', i velmean = np.mean(veldepth) #print velmean, count, count/pointsPerR for vel in veldepth: - outfile.writelines('%10f\n' %(vel - velmean)) + outfile.write('%10f\n' %(vel - velmean)) veldepth = [] elif absOrRel == 'rel': nref, dref, sref, velref = _readVgrid(inputfileref) @@ -372,7 +544,7 @@ def vgrids2VTK(inputfile='vgrids.in', outputfile='vgrids.vtk', absOrRel='abs', i return print("Writing velocity values to VTK file...") for velocity in velrel: - outfile.writelines('%10f\n' % velocity) + outfile.write('%10f\n' % velocity) print('Pertubations: min: %s %%, max: %s %%' % (min(velrel), max(velrel))) outfile.close() @@ -430,29 +602,29 @@ def rays2VTK(fnin, fdirout='./vtk_files/', nthPoint=50): # write header # print("Writing header for VTK file...") print("Writing shot %d to file %s" % (shotnumber, fnameout)) - outfile.writelines('# vtk DataFile Version 3.1\n') - outfile.writelines('FMTOMO rays\n') - outfile.writelines('ASCII\n') - outfile.writelines('DATASET POLYDATA\n') - outfile.writelines('POINTS %15d float\n' % (nPoints)) + outfile.write('# vtk DataFile Version 3.1\n') + outfile.write('FMTOMO rays\n') + outfile.write('ASCII\n') + outfile.write('DATASET POLYDATA\n') + outfile.write('POINTS %15d float\n' % (nPoints)) # write coordinates # print("Writing coordinates to VTK file...") for raynumber in rays[shotnumber].keys(): for raypoint in rays[shotnumber][raynumber]: - outfile.writelines('%10f %10f %10f \n' % (raypoint[0], raypoint[1], raypoint[2])) + outfile.write('%10f %10f %10f \n' % (raypoint[0], raypoint[1], raypoint[2])) - outfile.writelines('LINES %15d %15d\n' % (len(rays[shotnumber]), len(rays[shotnumber]) + nPoints)) + outfile.write('LINES %15d %15d\n' % (len(rays[shotnumber]), len(rays[shotnumber]) + nPoints)) # write indices # print("Writing indices to VTK file...") count = 0 for raynumber in rays[shotnumber].keys(): - outfile.writelines('%d ' % (len(rays[shotnumber][raynumber]))) + outfile.write('%d ' % (len(rays[shotnumber][raynumber]))) for index in range(len(rays[shotnumber][raynumber])): - outfile.writelines('%d ' % (count)) + outfile.write('%d ' % (count)) count += 1 - outfile.writelines('\n') + outfile.write('\n') def _readVgrid(filename): @@ -594,10 +766,10 @@ def addCheckerboard(spacing=10., pertubation=0.1, inputfile='vgrids.in', nPoints = nR * nTheta * nPhi # write header for velocity grid file (in RADIANS) - outfile.writelines('%10s %10s \n' % (1, 1)) - outfile.writelines('%10s %10s %10s\n' % (nR, nTheta, nPhi)) - outfile.writelines('%10s %10s %10s\n' % (dR, np.deg2rad(dTheta), np.deg2rad(dPhi))) - outfile.writelines('%10s %10s %10s\n' % (sR, np.deg2rad(sTheta), np.deg2rad(sPhi))) + outfile.write('%10s %10s \n' % (1, 1)) + outfile.write('%10s %10s %10s\n' % (nR, nTheta, nPhi)) + outfile.write('%10s %10s %10s\n' % (dR, np.deg2rad(dTheta), np.deg2rad(dPhi))) + outfile.write('%10s %10s %10s\n' % (sR, np.deg2rad(sTheta), np.deg2rad(sPhi))) spacR = correctSpacing(spacing, dR, '[meter], R') spacTheta = correctSpacing(_getAngle(spacing), dTheta, '[degree], Theta') @@ -642,7 +814,7 @@ def addCheckerboard(spacing=10., pertubation=0.1, inputfile='vgrids.in', evenOdd = evenOddR * evenOddT * evenOddP * ampFactor velocity += evenOdd * pertubation * velocity - outfile.writelines('%10s %10s\n' % (velocity, decm)) + outfile.write('%10s %10s\n' % (velocity, decm)) count += 1 progress = float(count) / float(nPoints) * 100 @@ -697,10 +869,10 @@ def addBox(x=(None, None), y=(None, None), z=(None, None), nPoints = nR * nTheta * nPhi # write header for velocity grid file (in RADIANS) - outfile.writelines('%10s %10s \n' % (1, 1)) - outfile.writelines('%10s %10s %10s\n' % (nR, nTheta, nPhi)) - outfile.writelines('%10s %10s %10s\n' % (dR, np.deg2rad(dTheta), np.deg2rad(dPhi))) - outfile.writelines('%10s %10s %10s\n' % (sR, np.deg2rad(sTheta), np.deg2rad(sPhi))) + outfile.write('%10s %10s \n' % (1, 1)) + outfile.write('%10s %10s %10s\n' % (nR, nTheta, nPhi)) + outfile.write('%10s %10s %10s\n' % (dR, np.deg2rad(dTheta), np.deg2rad(dPhi))) + outfile.write('%10s %10s %10s\n' % (sR, np.deg2rad(sTheta), np.deg2rad(sPhi))) count = 0 for radius in rGrid: @@ -722,7 +894,7 @@ def addBox(x=(None, None), y=(None, None), z=(None, None), if rFlag * thetaFlag * phiFlag is not 0: velocity = boxvelocity - outfile.writelines('%10s %10s\n' % (velocity, decm)) + outfile.write('%10s %10s\n' % (velocity, decm)) count += 1 progress = float(count) / float(nPoints) * 100 diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 275587a0..3fc10b9b 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -134,7 +134,8 @@ class SeisArray(object): if traceID2 < traceID1: direction = -1 return direction - print "Error: Same Value for traceID1 = %s and traceID2 = %s" % (traceID1, traceID2) + err_msg = "Same Value for traceID1 = %s and traceID2 = %s" % (traceID1, traceID2) + raise RuntimeError(err_msg) def _checkCoordDirection(self, traceID1, traceID2, coordinate): ''' @@ -146,7 +147,8 @@ class SeisArray(object): 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) + err_msg = "Same Value for traceID1 = %s and traceID2 = %s" % (traceID1, traceID2) + raise RuntimeError(err_msg) def _interpolateMeanDistances(self, traceID1, traceID2, coordinate): ''' @@ -455,7 +457,7 @@ class SeisArray(object): recx, recy, recz = self.getReceiverLists() nsrc = len(self.getSourceLocations()) - outfile.writelines('%s\n' % (len(zip(recx, recy, recz)) * nsrc)) + outfile.write('%s\n' % (len(zip(recx, recy, recz)) * nsrc)) for index in range(nsrc): for point in zip(recx, recy, recz): @@ -463,10 +465,10 @@ class SeisArray(object): rad = - rz lat = self._getAngle(ry) lon = self._getAngle(rx) - outfile.writelines('%15s %15s %15s\n' % (rad, lat, lon)) - outfile.writelines('%15s\n' % (1)) - outfile.writelines('%15s\n' % (index + 1)) - outfile.writelines('%15s\n' % (1)) + outfile.write('%15s %15s %15s\n' % (rad, lat, lon)) + outfile.write('%15s\n' % (1)) + outfile.write('%15s\n' % (index + 1)) + outfile.write('%15s\n' % (1)) outfile.close() @@ -506,22 +508,22 @@ class SeisArray(object): deltaPhi = abs(phiE - phiW) / float((nPhi - 1)) # write header for interfaces grid file (in RADIANS) - outfile.writelines('%10s\n' % (nInterfaces)) - outfile.writelines('%10s %10s\n' % (nTheta + 2, nPhi + 2)) # +2 cushion nodes - outfile.writelines('%10s %10s\n' % (np.deg2rad(deltaTheta), np.deg2rad(deltaPhi))) - outfile.writelines('%10s %10s\n' % (np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi))) + outfile.write('%10s\n' % (nInterfaces)) + outfile.write('%10s %10s\n' % (nTheta + 2, nPhi + 2)) # +2 cushion nodes + outfile.write('%10s %10s\n' % (np.deg2rad(deltaTheta), np.deg2rad(deltaPhi))) + outfile.write('%10s %10s\n' % (np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi))) interface1 = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method=method) interface2 = self.interpolateOnRegularGrid(nTheta, nPhi, thetaSN, phiWE, -depthmax, method=method) for point in interface1: z = point[2] - outfile.writelines('%10s\n' % (z + R)) + outfile.write('%10s\n' % (z + R)) - outfile.writelines('\n') + outfile.write('\n') for point in interface2: z = point[2] - outfile.writelines('%10s\n' % (z + R)) + outfile.write('%10s\n' % (z + R)) outfile.close() @@ -596,10 +598,10 @@ class SeisArray(object): deltaPhi = abs(phiE - phiW) / float(nPhi - 1) deltaR = abs(rbot - rtop) / float(nR - 1) - outfile.writelines('%10s %10s %10s\n' % (nR, nTheta, nPhi)) - outfile.writelines('%10s %10s %10s\n' % (deltaR, deltaTheta, deltaPhi)) - outfile.writelines('%10s %10s %10s\n' % (rtop, thetaS, phiW)) - outfile.writelines('%10s %10s\n' % refinement) + outfile.write('%10s %10s %10s\n' % (nR, nTheta, nPhi)) + outfile.write('%10s %10s %10s\n' % (deltaR, deltaTheta, deltaPhi)) + outfile.write('%10s %10s %10s\n' % (rtop, thetaS, phiW)) + outfile.write('%10s %10s\n' % refinement) outfile.close() @@ -708,10 +710,10 @@ class SeisArray(object): 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' % (deltaR, np.deg2rad(deltaTheta), np.deg2rad(deltaPhi))) - outfile.writelines( + outfile.write('%10s %10s \n' % (1, 1)) + outfile.write('%10s %10s %10s\n' % (nR + 2, nTheta + 2, nPhi + 2)) + outfile.write('%10s %10s %10s\n' % (deltaR, np.deg2rad(deltaTheta), np.deg2rad(deltaPhi))) + outfile.write( '%10s %10s %10s\n' % (rbot - deltaR, np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi))) surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method=method) @@ -746,15 +748,16 @@ class SeisArray(object): vbot[index] - vtop[index]) + vtop[index] break if not (ztop[index]) >= depth > (zbot[index]): - print( - 'ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile' % ( - depth - topo)) - return + err_msg = 'ERROR in grid inputfile, could not find velocity' + 'for a z-value of %s in the inputfile' % (depth - topo) + raise ValueError(err_msg) count += 1 - if vel < 0: - print( - 'ERROR, vel <0; z, topo, zbot, vbot, vtop:', depth, topo, zbot[index], vbot[index], vtop[index]) - outfile.writelines('%10s %10s\n' % (vel, decm)) + if not vel >= 0: + err_msg = 'vel < 0; z, topo, zbot, vbot, vtop:', + depth, topo, zbot[index], vbot[index], vtop[index] + raise ValueError(err_msg) + + outfile.write('%10s %10s\n' % (vel, decm)) progress = float(count) / float(nTotal) * 100 self._update_progress(progress) @@ -775,7 +778,7 @@ class SeisArray(object): 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)) + recfile_out.write('%5s %15s %15s %15s\n' % (traceID, x, y, z)) print "Exported coordinates for %s traces to file > %s" % (count, filename) recfile_out.close() @@ -895,11 +898,11 @@ class SeisArray(object): # write header print("Writing header for VTK file...") - outfile.writelines('# vtk DataFile Version 3.1\n') - outfile.writelines('Surface Points\n') - outfile.writelines('ASCII\n') - outfile.writelines('DATASET POLYDATA\n') - outfile.writelines('POINTS %15d float\n' % (nPoints)) + outfile.write('# vtk DataFile Version 3.1\n') + outfile.write('Surface Points\n') + outfile.write('ASCII\n') + outfile.write('DATASET POLYDATA\n') + outfile.write('POINTS %15d float\n' % (nPoints)) # write coordinates print("Writing coordinates to VTK file...") @@ -908,23 +911,23 @@ class SeisArray(object): y = point[1] z = point[2] - outfile.writelines('%10f %10f %10f \n' % (x, y, z)) + outfile.write('%10f %10f %10f \n' % (x, y, z)) - outfile.writelines('VERTICES %15d %15d\n' % (nPoints, 2 * nPoints)) + outfile.write('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.write('%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') + # outfile.write('POINT_DATA %15d\n' %(nPoints)) + # outfile.write('SCALARS traceIDs int %d\n' %(1)) + # outfile.write('LOOKUP_TABLE default\n') # # write traceIDs # print("Writing traceIDs to VTK file...") # for traceID in traceIDs: - # outfile.writelines('%10d\n' %traceID) + # outfile.write('%10d\n' %traceID) outfile.close() print("Wrote %d points to file: %s" % (nPoints, filename)) @@ -944,11 +947,11 @@ class SeisArray(object): # 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)) + outfile.write('# vtk DataFile Version 3.1\n') + outfile.write('Receivers with traceIDs\n') + outfile.write('ASCII\n') + outfile.write('DATASET POLYDATA\n') + outfile.write('POINTS %15d float\n' % (nPoints)) # write coordinates print("Writing coordinates to VTK file...") @@ -957,23 +960,23 @@ class SeisArray(object): y = self._getYreceiver(traceID) z = self._getZreceiver(traceID) - outfile.writelines('%10f %10f %10f \n' % (x, y, z)) + outfile.write('%10f %10f %10f \n' % (x, y, z)) - outfile.writelines('VERTICES %15d %15d\n' % (nPoints, 2 * nPoints)) + outfile.write('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.write('%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') + outfile.write('POINT_DATA %15d\n' % (nPoints)) + outfile.write('SCALARS traceIDs int %d\n' % (1)) + outfile.write('LOOKUP_TABLE default\n') # write traceIDs print("Writing traceIDs to VTK file...") for traceID in traceIDs: - outfile.writelines('%10d\n' % traceID) + outfile.write('%10d\n' % traceID) outfile.close() print("Wrote %d receiver for to file: %s" % (nPoints, filename)) @@ -993,11 +996,11 @@ class SeisArray(object): # 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)) + outfile.write('# vtk DataFile Version 3.1\n') + outfile.write('Shots with shotnumbers\n') + outfile.write('ASCII\n') + outfile.write('DATASET POLYDATA\n') + outfile.write('POINTS %15d float\n' % (nPoints)) # write coordinates print("Writing coordinates to VTK file...") @@ -1006,23 +1009,23 @@ class SeisArray(object): y = self._getYshot(shotnumber) z = self._getZshot(shotnumber) - outfile.writelines('%10f %10f %10f \n' % (x, y, z)) + outfile.write('%10f %10f %10f \n' % (x, y, z)) - outfile.writelines('VERTICES %15d %15d\n' % (nPoints, 2 * nPoints)) + outfile.write('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.write('%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') + outfile.write('POINT_DATA %15d\n' % (nPoints)) + outfile.write('SCALARS shotnumbers int %d\n' % (1)) + outfile.write('LOOKUP_TABLE default\n') # write shotnumber print("Writing shotnumbers to VTK file...") for shotnumber in shotnumbers: - outfile.writelines('%10d\n' % shotnumber) + outfile.write('%10d\n' % shotnumber) outfile.close() print("Wrote %d sources to file: %s" % (nPoints, filename)) From 1cb86d7b479fc31a7953e584050b8aee646665e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 11 May 2016 13:42:51 +0200 Subject: [PATCH 2/3] Cleaned input file, additional parameter apverbose --- inputs/autoPyLoT_regional.in | 38 ++++++++++++++++++++---------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/inputs/autoPyLoT_regional.in b/inputs/autoPyLoT_regional.in index 4248fd30..14972000 100644 --- a/inputs/autoPyLoT_regional.in +++ b/inputs/autoPyLoT_regional.in @@ -12,23 +12,27 @@ e1412.008.06 #eventID# %event ID for single event pr /DATA/Egelados/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure# %choose data structure 0 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything -AUTOPHASES_AIC_HOS4_ARH #phasefile# %name of autoPILOT output phase file -AUTOLOC_AIC_HOS4_ARH #locfile# %name of autoPILOT output location file -AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file containing polarities -HYPOSAT #locrt# %location routine used ("HYPOINVERSE" or "HYPOSAT") -6 #pmin# %minimum required P picks for location -4 #p0min# %minimum required P picks for location if at least - %3 excellent P picks are found -2 #smin# %minimum required S picks for location -/home/ludger/bin/run_HYPOSAT4autoPILOT.csh #cshellp# %path and name of c-shell script to run location routine -7.6 8.5 #blon# %longitude bounding for location map -49 49.4 #blat# %lattitude bounding for location map -#parameters for moment magnitude estimation# -5000 #vp# %average P-wave velocity -2800 #vs# %average S-wave velocity -2200 #rho# %rock density [kg/m^3] -300 #Qp# %quality factor for P waves -100 #Qs# %quality factor for S waves +True #apverbose# %choose 'True' or 'False' for terminal output +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +#NLLoc settings# +/home/ludger/NLLOC #nllocbin# %path to NLLoc executable +/home/ludger/NLLOC/Insheim #nllocroot# %root of NLLoc-processing directory +AUTOPHASES.obs #phasefile# %name of autoPyLoT-output phase file for NLLoc + %(in nllocroot/obs) +Insheim_min1d2015_auto.in #ctrfile# %name of autoPyLoT-output control file for NLLoc + %(in nllocroot/run) +ttime #ttpatter# %pattern of NLLoc ttimes from grid + %(in nllocroot/times) +AUTOLOC_nlloc #outpatter# %pattern of NLLoc-output file + %(returns 'eventID_outpatter') +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +#parameters for seismic moment estimation# +3530 #vp# %average P-wave velocity +2700 #rho# %average rock density [kg/m^3] +1000f**0.8 #Qp# %quality factor for P waves (Qp*f^a) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file containing derived polarities +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #common settings picker# 20 #pstart# %start time [s] for calculating CF for P-picking 100 #pstop# %end time [s] for calculating CF for P-picking From 2bb2ddbc7b5027897781ae29c39c9dfba2f1343a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 11 May 2016 13:44:26 +0200 Subject: [PATCH 3/3] Cleaned input file --- inputs/autoPyLoT_local.in | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/inputs/autoPyLoT_local.in b/inputs/autoPyLoT_local.in index 725a9dec..d01a9c0d 100644 --- a/inputs/autoPyLoT_local.in +++ b/inputs/autoPyLoT_local.in @@ -6,8 +6,8 @@ #main settings# /DATA/Insheim #rootpath# %project path EVENT_DATA/LOCAL #datapath# %data path -2015.11_Insheim #database# %name of data base - #eventID# %event ID for single event processing +2013.02_Insheim #database# %name of data base +e0019.048.13 #eventID# %event ID for single event processing /DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure#%choose data structure 0 #iplot# %flag for plotting: 0 none, 1 partly, >1 everything @@ -30,13 +30,7 @@ AUTOLOC_nlloc #outpatter# %pattern of NLLoc-output file 2500 #rho# %average rock density [kg/m^3] 300f**0.8 #Qp# %quality factor for P waves (Qp*f^a) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file containing polarities -6 #pmin# %minimum required P picks for location -4 #p0min# %minimum required P picks for location if at least - %3 excellent P picks are found -2 #smin# %minimum required S picks for location -7.6 8.5 #blon# %longitude bounding for location map -49 49.4 #blat# %lattitude bounding for location map +AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file containing derived polarities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #common settings picker# 15.0 #pstart# %start time [s] for calculating CF for P-picking @@ -62,7 +56,6 @@ HOS #algoP# %choose algorithm for P-onset 3 0.1 0.5 0.5 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s] 3.0 #pickwinP# %for initial AIC pick, length of P-pick window [s] 6.0 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick) -0 #peps4aic# %for HOS/AR, artificial uplift of samples of AIC-function (P) 0.2 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s] 0.1 #tsmoothP# %for HOS/AR, take average of samples for smoothing CF [s] 0.001 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P)