From 30680a782021715cb79671f925ab49fca0a9acf4 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Thu, 22 Oct 2015 10:54:24 +0200 Subject: [PATCH 01/10] cosmetics for array plots (title, legend) --- pylot/core/active/seismicArrayPreparation.py | 41 +++++++++++++------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 7e7dd09a..4077a3d6 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -16,13 +16,14 @@ class SeisArray(object): Note: Source and Receiver files for FMTOMO will be generated by the Survey object (because traveltimes will be added directly). ''' def __init__(self, recfile): + self.recfile = recfile self._receiverlines = {} self._receiverCoords = {} self._measuredReceivers = {} self._measuredTopo = {} self._sourceLocs = {} self._geophoneNumbers = {} - self._receiverlist = open(recfile, 'r').readlines() + self._receiverlist = open(self.recfile, 'r').readlines() self._generateReceiverlines() self._setReceiverCoords() self._setGeophoneNumbers() @@ -472,28 +473,33 @@ class SeisArray(object): def plotArray2D(self, plot_topo = False, highlight_measured = False, annotations = True): import matplotlib.pyplot as plt plt.interactive(True) - plt.figure() + fig = plt.figure() + ax = plt.axes() xmt, ymt, zmt = self.getMeasuredTopoLists() xsc, ysc, zsc = self.getSourceLocsLists() xmr, ymr, zmr = self.getMeasuredReceiverLists() xrc, yrc, zrc = self.getReceiverLists() - plt.plot(xrc, yrc, 'k.', markersize = 10, label = 'all receivers') - plt.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations') + if len(xrc) > 0: + ax.plot(xrc, yrc, 'k.', markersize = 10, label = 'all receivers') + if len(xsc) > 0: + ax.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations') if plot_topo == True: - plt.plot(xmt, ymt, 'b', markersize = 10, label = 'measured topo points') + ax.plot(xmt, ymt, 'b', markersize = 10, label = 'measured topo points') if highlight_measured == True: - plt.plot(xmr, ymr, 'ro', label = 'measured receivers') + ax.plot(xmr, ymr, 'ro', label = 'measured receivers') - plt.xlabel('X [m]') - plt.ylabel('Y [m]') + plt.title('2D plot of seismic array %s'%self.recfile) + ax.set_xlabel('X [m]') + ax.set_ylabel('Y [m]') + ax.set_aspect('equal') plt.legend() if annotations == True: for traceID in self.getReceiverCoordinates().keys(): - plt.annotate((' ' + str(traceID)), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small', color = 'k') + ax.annotate((' ' + str(traceID)), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small', color = 'k') for shotnumber in self.getSourceLocations().keys(): - plt.annotate((' ' + str(shotnumber)), xy = (self._getXshot(shotnumber), self._getYshot(shotnumber)), fontsize = 'x-small', color = 'b') + ax.annotate((' ' + str(shotnumber)), xy = (self._getXshot(shotnumber), self._getYshot(shotnumber)), fontsize = 'x-small', color = 'b') @@ -508,11 +514,18 @@ class SeisArray(object): xmt, ymt, zmt = self.getMeasuredTopoLists() xmr, ymr, zmr = self.getMeasuredReceiverLists() - xin, yin, zin = self.getReceiverLists() + xrc, yrc, zrc = self.getReceiverLists() + xsc, ysc, zsc = self.getSourceLocsLists() - ax.plot(xmt, ymt, zmt, 'b*', markersize = 10, label = 'measured topo points') - ax.plot(xin, yin, zin, 'k.', markersize = 10, label = 'interpolated receivers') - ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers') + plt.title('3D plot of seismic array %s'%self.recfile) + if len(xmt) > 0: + ax.plot(xmt, ymt, zmt, 'b.', markersize = 10, label = 'measured topo points') + if len(xrc) > 0: + ax.plot(xrc, yrc, zrc, 'k.', markersize = 10, label = 'all receivers') + if len(xmr) > 0: + ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers') + if len(xsc) > 0: + ax.plot(xsc, ysc, zsc, 'b*', label = 'shot locations') ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation') ax.legend() From 84f1639e59bfea015b40bab6d9c89716dd25e02c Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Fri, 30 Oct 2015 11:56:45 +0100 Subject: [PATCH 02/10] added transformation to vgrids vtk with relative values using vgrids.in and vgridsref.in from FMTOMO --- pylot/core/active/fmtomo2vtk.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/pylot/core/active/fmtomo2vtk.py b/pylot/core/active/fmtomo2vtk.py index 95b90705..63990852 100644 --- a/pylot/core/active/fmtomo2vtk.py +++ b/pylot/core/active/fmtomo2vtk.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- import numpy as np -def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'): +def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk', absOrRel = 'abs', inputfileref = 'vgridsref.in'): ''' Generate a vtk-file readable by e.g. paraview from FMTOMO output vgrids.in ''' @@ -19,6 +19,8 @@ def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'): nTheta = int(vglines[1].split()[1]) nPhi = int(vglines[1].split()[2]) + print('readNumberOf Points: Awaiting %d grid points in %s' + %(nR*nTheta*nPhi, filename)) fin.close() return nR, nTheta, nPhi @@ -95,9 +97,24 @@ def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'): outfile.writelines('LOOKUP_TABLE default\n') # write velocity - print("Writing velocity values to VTK file...") - for velocity in vel: - outfile.writelines('%10f\n' %velocity) + if absOrRel == 'abs': + print("Writing velocity values to VTK file...") + for velocity in vel: + outfile.writelines('%10f\n' %velocity) + elif absOrRel == 'rel': + velref = readVelocity(inputfileref) + if not len(velref) == len(vel): + print('ERROR: Number of gridpoints mismatch for %s and %s'%(inputfile, inputfileref)) + return + velrel = [vel - velref for vel, velref in zip(vel, velref)] + nR_ref, nTheta_ref, nPhi_ref = readNumberOfPoints(inputfileref) + if not nR_ref == nR and nTheta_ref == nTheta and nPhi_ref == nPhi: + print('ERROR: Dimension mismatch of grids %s and %s'%(inputfile, inputfileref)) + return + print("Writing velocity values to VTK file...") + for velocity in velrel: + outfile.writelines('%10f\n' %velocity) + print('Pertubations: min: %s, max: %s'%(min(velrel), max(velrel))) outfile.close() print("Wrote velocity grid for %d points to file: %s" %(nPoints, outputfile)) From 422a76012e3be7b16461485cddb6f49c46adba79 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Fri, 30 Oct 2015 12:00:03 +0100 Subject: [PATCH 03/10] vert. exag for surface plot --- pylot/core/active/seismicArrayPreparation.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 4077a3d6..212f8495 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -470,7 +470,7 @@ class SeisArray(object): print "Exported coordinates for %s traces to file > %s" %(count, filename) recfile_out.close() - def plotArray2D(self, plot_topo = False, highlight_measured = False, annotations = True): + def plotArray2D(self, plot_topo = False, highlight_measured = False, annotations = True, pointsize = 10): import matplotlib.pyplot as plt plt.interactive(True) fig = plt.figure() @@ -481,14 +481,14 @@ class SeisArray(object): xrc, yrc, zrc = self.getReceiverLists() if len(xrc) > 0: - ax.plot(xrc, yrc, 'k.', markersize = 10, label = 'all receivers') + ax.plot(xrc, yrc, 'k.', markersize = pointsize, label = 'all receivers') if len(xsc) > 0: - ax.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations') + ax.plot(xsc, ysc, 'b*', markersize = pointsize, label = 'shot locations') if plot_topo == True: - ax.plot(xmt, ymt, 'b', markersize = 10, label = 'measured topo points') + ax.plot(xmt, ymt, 'b.', markersize = pointsize, label = 'measured topo points') if highlight_measured == True: - ax.plot(xmr, ymr, 'ro', label = 'measured receivers') + ax.plot(xmr, ymr, 'r.', markersize = pointsize, label = 'measured receivers') plt.title('2D plot of seismic array %s'%self.recfile) ax.set_xlabel('X [m]') @@ -532,7 +532,7 @@ class SeisArray(object): return ax - def plotSurface3D(self, ax = None, step = 0.5, method = 'linear'): + def plotSurface3D(self, ax = None, step = 0.5, method = 'linear', exag = False): from matplotlib import cm import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D @@ -558,7 +558,8 @@ class SeisArray(object): ax.plot_surface(xgrid, ygrid, zgrid, linewidth = 0, cmap = cm.jet, vmin = min(z), vmax = max(z)) - ax.set_zlim(-(max(x) - min(x)/2),(max(x) - min(x)/2)) + if exag == False: + ax.set_zlim(-(max(x) - min(x)/2),(max(x) - min(x)/2)) ax.set_aspect('equal') ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation') From 1b95ed0da7ff3e00d3fa264335f7ef74464c4569 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Fri, 30 Oct 2015 12:37:41 +0100 Subject: [PATCH 04/10] relative velocities in percent --- pylot/core/active/fmtomo2vtk.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pylot/core/active/fmtomo2vtk.py b/pylot/core/active/fmtomo2vtk.py index 63990852..10cbe879 100644 --- a/pylot/core/active/fmtomo2vtk.py +++ b/pylot/core/active/fmtomo2vtk.py @@ -106,7 +106,15 @@ def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk', absOrRel = 'a if not len(velref) == len(vel): print('ERROR: Number of gridpoints mismatch for %s and %s'%(inputfile, inputfileref)) return - velrel = [vel - velref for vel, velref in zip(vel, velref)] + #velrel = [((vel - velref) / velref * 100) for vel, velref in zip(vel, velref)] + velrel = [] + for velocities in zip(vel, velref): + v, vref = velocities + if not vref == 0: + velrel.append((v - vref) / vref * 100) + else: + velrel.append(0) + nR_ref, nTheta_ref, nPhi_ref = readNumberOfPoints(inputfileref) if not nR_ref == nR and nTheta_ref == nTheta and nPhi_ref == nPhi: print('ERROR: Dimension mismatch of grids %s and %s'%(inputfile, inputfileref)) From bcc6c8a73d599606e8fae71beffdfd016eb5522c Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Tue, 3 Nov 2015 10:25:52 +0100 Subject: [PATCH 05/10] enabled use of parameter file "mygrid.in" for generation of a starting model, prepared generation of vgrid model from array dimensions --- pylot/core/active/seismicArrayPreparation.py | 80 +++++++++++++++----- 1 file changed, 63 insertions(+), 17 deletions(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 212f8495..21a5c625 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -364,9 +364,9 @@ class SeisArray(object): return surface def generateVgrid(self, nTheta = 80, nPhi = 80, nR = 120, - thetaSN = (-0.2, 1.2), phiWE = (-0.2, 1.2), - Rbt = (-62.0, 6.0), vbot = 5.5, filename = 'vgrids.in', - method = 'linear' ): + Rbt = (-62.0, 6.0), thetaSN = None, + phiWE = None, outfilename = 'vgrids.in', + method = 'linear', infilename = 'mygrid.in'): ''' Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s]. @@ -387,9 +387,12 @@ class SeisArray(object): :param: Rbt (bot, top) extensions of the model in km type: tuple - + :param: vbot, velocity at the bottom of the model type: real + + :param: method, interpolation method for topography + type: str ''' def getRad(angle): @@ -397,16 +400,48 @@ class SeisArray(object): rad = angle / 180 * PI return rad - def getZmax(surface): - z = [] - for point in surface: - z.append(point[2]) - return max(z) + def readMygridNlayers(filename): + infile = open(filename, 'r') + nlayers = len(infile.readlines()) / 2 + infile.close() + + return nlayers + + def readMygrid(filename): + ztop = []; zbot = []; vtop = []; vbot = [] + infile = open(filename, 'r') + nlayers = readMygridNlayers(filename) + + for index in range(nlayers): + line1 = infile.readline() + line2 = infile.readline() + ztop.append(float(line1.split()[0])) + vtop.append(float(line1.split()[1])) + zbot.append(float(line2.split()[0])) + vbot.append(float(line2.split()[1])) + + if not ztop[0] == 0: + print('ERROR: there must be a velocity set for z = 0 in the file %s'%filename) + print('e.g.:\n0 0.33\n-5 1.0\netc.') + + infile.close() + return ztop, zbot, vtop, vbot R = 6371 vmin = 0.34 + cushionfactor = 0.1 # add some extra space to the model decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3) - outfile = open(filename, 'w') + outfile = open(outfilename, 'w') + + # generate dimensions of the grid from array + if thetaSN is None and phiWE is None: + x, y, z = self.getAllMeasuredPointsLists() + phi_min, phi_max = (self._getAngle(min(x)), self._getAngle(max(x))) + theta_min, theta_max = (self._getAngle(min(y)), self._getAngle(max(y))) + cushionPhi = abs(phi_max - phi_min) * cushionfactor + cushionTheta = abs(theta_max - theta_min) * cushionfactor + phiWE = (phi_min - cushionPhi, phi_max + cushionPhi) + thetaSN = (theta_min - cushionTheta, theta_max + cushionTheta) thetaS, thetaN = thetaSN phiW, phiE = phiWE @@ -433,11 +468,14 @@ class SeisArray(object): outfile.writelines('%10s %10s %10s\n' %(rbot - rDelta, getRad(thetaS - thetaDelta), getRad(phiW - phiDelta))) surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method, filename = None) - zmax = getZmax(surface) - print "\nGenerating velocity grid for FMTOMO. Output filename = %s, interpolation method = %s"%(filename, method) + print "\nGenerating velocity grid for FMTOMO. Output filename = %s, interpolation method = %s"%(outfilename, method) print "nTheta = %s, nPhi = %s, nR = %s, thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt) count = 0 + + nlayers = readMygridNlayers(infilename) + ztop, zbot, vtop, vbot = readMygrid(infilename) + for radius in rGrid: for theta in thetaGrid: for phi in phiGrid: @@ -445,19 +483,27 @@ class SeisArray(object): yval = self._getDistance(theta) for point in surface: if point[0] == xval and point[1] == yval: - z = point[2] - if radius > (R + z + 1): + topo = point[2] + z = -(R + topo - radius) + if z > (topo + 1): vel = 0.0 - # elif radius > (R + z - 15): ########### TESTING - # vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * 1.0 + vmin ########################## + elif (topo + 1) >= z > (topo): # cushioning around topography + vel = vtop[0] else: - vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * vbot + vmin ########################## + for index in range(nlayers): + if (topo + ztop[index]) >= z > (topo + zbot[index]): + vel = (z - topo) / (zbot[index] - topo) * (vbot[index] - vtop[index]) + vtop[index] + break + if not (topo + ztop[index]) >= z > (topo + zbot[index]): + print('ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile'%(z - topo)) + return count += 1 outfile.writelines('%10s %10s\n'%(vel, decm)) progress = float(count) / float(nTotal) * 100 self._update_progress(progress) + print('Wrote %d points to file %s for %d layers'%(count, outfilename, nlayers)) outfile.close() def exportAll(self, filename = 'interpolated_receivers.out'): From 1c08152e65de55e41067fdaf5ad1a8cdf9a393fe Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Tue, 3 Nov 2015 10:27:06 +0100 Subject: [PATCH 06/10] vkt file "rel" changed to changes in percent --- pylot/core/active/fmtomo2vtk.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pylot/core/active/fmtomo2vtk.py b/pylot/core/active/fmtomo2vtk.py index 10cbe879..59c0968e 100644 --- a/pylot/core/active/fmtomo2vtk.py +++ b/pylot/core/active/fmtomo2vtk.py @@ -93,7 +93,10 @@ def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk', absOrRel = 'a outfile.writelines('SPACING %f %f %f\n' %(dX, dY, dZ)) outfile.writelines('POINT_DATA %15d\n' %(nPoints)) - outfile.writelines('SCALARS velocity float %d\n' %(1)) + if absOrRel == 'abs': + outfile.writelines('SCALARS velocity float %d\n' %(1)) + elif absOrRel == 'rel': + outfile.writelines('SCALARS velChangePercent float %d\n' %(1)) outfile.writelines('LOOKUP_TABLE default\n') # write velocity From c28641d62902cbf3f34a80d3a01a51f29aef2ed1 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Tue, 3 Nov 2015 10:29:09 +0100 Subject: [PATCH 07/10] changed plot traces for presentation (SNR, _drawStream etc.) --- pylot/core/active/seismicshot.py | 56 +++++++++++++++++++++++++------- 1 file changed, 44 insertions(+), 12 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 07cd2b43..58267dd5 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -588,6 +588,23 @@ class SeismicShot(object): # plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_lowerb, ':k') # plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_upperb, ':k') + def plotTrace(self, traceID, plotSNR = True, lw = 1): + fig = plt.figure() + ax = fig.add_subplot(111) + ax = self._drawStream(traceID, ax = ax) + + tgap = self.getTgap() + tsignal = self.getTsignal() + pick = self.getPick(traceID) + tnoise = pick - tgap + snr, snrdb, noiselevel = self.getSNR(traceID) + + ax.plot([0, tnoise], [noiselevel, noiselevel], 'm', linewidth = lw, label = 'noise level') + ax.plot([tnoise, pick], [noiselevel, noiselevel], 'g:', linewidth = lw, label = 'gap') + ax.plot([tnoise + tgap, pick + tsignal], [noiselevel * snr, noiselevel * snr], 'b', linewidth = lw, label = 'signal level') + ax.legend() + ax.text(0.05, 0.95, 'SNR: %s' %snr, transform = ax.transAxes) + def plot_traces(self, traceID, folm = 0.6): ########## 2D, muss noch mehr verbessert werden ########## from matplotlib.widgets import Button @@ -621,7 +638,7 @@ class SeismicShot(object): self._drawStream(traceID) self._drawCFs(traceID, folm) - def _drawStream(self, traceID, refresh = False): + def _drawStream(self, traceID, refresh = False, ax = None): from pylot.core.util.utils import getGlobalTimes from pylot.core.util.utils import prepTimeAxis @@ -630,7 +647,8 @@ class SeismicShot(object): timeaxis = prepTimeAxis(stime, stream[0]) timeaxis -= stime - ax = self.traces4plot[traceID]['ax1'] + if ax is None: + ax = self.traces4plot[traceID]['ax1'] if refresh == True: xlim, ylim = ax.get_xlim(), ax.get_ylim() @@ -645,8 +663,9 @@ class SeismicShot(object): ax.plot([self.getPick(traceID), self.getPick(traceID)], [min(stream[0].data), max(stream[0].data)], - 'r', label = 'mostlikely') + 'r', label = 'most likely') ax.legend() + return ax def _drawCFs(self, traceID, folm, refresh = False): hoscf = self.getHOScf(traceID) @@ -665,7 +684,7 @@ class SeismicShot(object): ax.plot([self.getPick(traceID), self.getPick(traceID)], [min(np.minimum(hoscf.getCF(), aiccf.getCF())), max(np.maximum(hoscf.getCF(), aiccf.getCF()))], - 'r', label = 'mostlikely') + 'r', label = 'most likely') ax.plot([0, self.getPick(traceID)], [folm * max(hoscf.getCF()), folm * max(hoscf.getCF())], 'm:', label = 'folm = %s' %folm) @@ -745,11 +764,12 @@ class SeismicShot(object): :type: 'logical' ''' from scipy.interpolate import griddata + from matplotlib import cm + cmap = cm.jet x = []; xcut = [] y = []; ycut = [] z = []; zcut = [] - tmin, tmax = self.getCut() for traceID in self.pick.keys(): if self.getFlag(traceID) != 0: @@ -761,6 +781,9 @@ class SeismicShot(object): ycut.append(self.getRecLoc(traceID)[1]) zcut.append(self.getPickIncludeRemoved(traceID)) + tmin = 0.8 * min(z) # 20% cushion for colorbar + tmax = 1.2 * max(z) + xaxis = np.arange(min(x), max(x), step) yaxis = np.arange(min(y), max(y), step) xgrid, ygrid = np.meshgrid(xaxis, yaxis) @@ -770,19 +793,28 @@ class SeismicShot(object): fig = plt.figure() ax = plt.axes() - ax.matshow(zgrid, extent = [min(x), max(x), min(y), max(y)], origin = 'lower') + count = 0 + ax.imshow(zgrid, extent = [min(x), max(x), min(y), max(y)], vmin = tmin, vmax = tmax, cmap = cmap, origin = 'lower', alpha = 0.85) plt.text(0.45, 0.9, 'shot: %s' %self.getShotnumber(), transform = ax.transAxes) - sc = ax.scatter(x, y, c = z, s = 30, label = 'picked shots', vmin = tmin, vmax = tmax, linewidths = 1.5) - sccut = ax.scatter(xcut, ycut, c = zcut, s = 30, edgecolor = 'm', label = 'cut out shots', vmin = tmin, vmax = tmax, linewidths = 1.5) + sc = ax.scatter(x, y, c = z, s = 30, label = 'picked shots', vmin = tmin, vmax = tmax, cmap = cmap, linewidths = 1.5) + for xyz in zip(xcut, ycut, zcut): + x, y, z = xyz + label = None + if z > tmax: + count += 1 + z = 'w' + if count == 1: + label = 'cut out shots' + ax.scatter(x, y, c = z, s = 30, edgecolor = 'm', label = label, vmin = tmin, vmax = tmax, cmap = cmap, linewidths = 1.5) if colorbar == True: - plt.colorbar(sc) + cbar = plt.colorbar(sc) + cbar.set_label('Time [s]') + + ax.legend() ax.set_xlabel('X') ax.set_ylabel('Y') ax.plot(self.getSrcLoc()[0], self.getSrcLoc()[1],'*k', markersize = 15) # plot source location - if plotRec == True: - ax.scatter(x, y, c = z, s = 30) - if annotations == True: for traceID in self.getTraceIDlist(): if self.getFlag(traceID) is not 0: From dc4d19ba884202e912f66b33cd4131fb7d17e90a Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Wed, 4 Nov 2015 10:04:05 +0100 Subject: [PATCH 08/10] bugfixes: formula for gradients in multiple layers for vgrid corrected, fixed issue with integer division that would lead to a wrong spacing --- pylot/core/active/seismicArrayPreparation.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 21a5c625..c0c1d99b 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -449,9 +449,9 @@ class SeisArray(object): rtop = Rbt[1] + R # need to determine the delta to add two cushion nodes around the min/max values - thetaDelta = abs(thetaN - thetaS) / (nTheta - 1) - phiDelta = abs(phiE - phiW) / (nPhi - 1) - rDelta = abs(rbot - rtop) / (nR - 1) + thetaDelta = abs(thetaN - thetaS) / float((nTheta - 1)) + phiDelta = abs(phiE - phiW) / float((nPhi - 1)) + rDelta = abs(rbot - rtop) / float((nR - 1)) # create a regular grid including +2 cushion nodes in every direction thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes @@ -487,17 +487,19 @@ class SeisArray(object): z = -(R + topo - radius) if z > (topo + 1): vel = 0.0 - elif (topo + 1) >= z > (topo): # cushioning around topography + elif (topo + 1) >= z > topo: # cushioning around topography vel = vtop[0] else: for index in range(nlayers): if (topo + ztop[index]) >= z > (topo + zbot[index]): - vel = (z - topo) / (zbot[index] - topo) * (vbot[index] - vtop[index]) + vtop[index] + vel = (z - ztop[index]) / (zbot[index] - ztop[index]) * (vbot[index] - vtop[index]) + vtop[index] break if not (topo + ztop[index]) >= z > (topo + zbot[index]): print('ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile'%(z - topo)) return count += 1 + if vel < 0: + print('ERROR, vel <0; z, topo, zbot, vbot, vtop:', z, topo, zbot[index], vbot[index], vtop[index]) outfile.writelines('%10s %10s\n'%(vel, decm)) progress = float(count) / float(nTotal) * 100 From d8a59ac81dace573f002c7532909566f2c082ae5 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Wed, 4 Nov 2015 10:13:52 +0100 Subject: [PATCH 09/10] bugfix: not yet done: problem with cusioning around topography --- pylot/core/active/seismicArrayPreparation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index c0c1d99b..58d021be 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -486,7 +486,7 @@ class SeisArray(object): topo = point[2] z = -(R + topo - radius) if z > (topo + 1): - vel = 0.0 + vel = vtop[0] elif (topo + 1) >= z > topo: # cushioning around topography vel = vtop[0] else: From 223902f2d4a9fe992730bb1ba925bfa432f26550 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 4 Nov 2015 13:11:48 +0100 Subject: [PATCH 10/10] reformatted the file to avoid usage of tabs and spaces and made a more pythonic splash message than multiple calls to print --- autoPyLoT.py | 42 ++++++++++++++++++++---------------------- 1 file changed, 20 insertions(+), 22 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 1e60e54d..e7f79057 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -32,19 +32,17 @@ def autoPyLoT(inputfile): .. rubric:: Example """ - print '************************************' - print '*********autoPyLoT starting*********' - print 'The Python picking and Location Tool' - print ' Version ', _getVersionString(), '2015' - print ' ' - print 'Authors:' - print 'S. Wehling-Benatelli' - print ' Ruhr-Universität Bochum' - print 'L. Küperkoch' - print ' BESTEC GmbH, Landau (Pfalz)' - print 'K. Olbert' - print ' Christian-Albrechts Universität Kiel' - print '************************************' + splash = '''************************************\n + *********autoPyLoT starting*********\n + The Python picking and Location Tool\n + Version {version} 2015\n + \n + Authors:\n + S. Wehling-Benatelli (Ruhr-Universität Bochum)\n + L. Küperkoch (BESTEC GmbH, Landau i. d. Pfalz)\n + K. Olbert (Christian-Albrechts Universität zu Kiel)\n + ***********************************'''.format(version=_getVersionString()) + print(splash) # reading parameter file @@ -77,7 +75,7 @@ def autoPyLoT(inputfile): # get path to NLLoc executable nllocbin = parameter.getParam('nllocbin') nlloccall = '%s/NLLoc' % nllocbin - # get name of phase file + # get name of phase file phasef = parameter.getParam('phasefile') phasefile = '%s/obs/%s' % (nllocroot, phasef) # get name of NLLoc-control file @@ -99,8 +97,8 @@ def autoPyLoT(inputfile): if not parameter.hasParam('eventID'): for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]: data.setWFData(glob.glob(os.path.join(datapath, event, '*'))) - print 'Working on event %s' % event - print data + print('Working on event %s' % event) + print(data) wfdat = data.getWFData() # all available streams ########################################################## @@ -148,7 +146,7 @@ def autoPyLoT(inputfile): ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) - + ########################################################## # locating if locflag == 1: @@ -166,15 +164,15 @@ def autoPyLoT(inputfile): filedata = nllfile.read() if filedata.find(locfiles) < 0: # replace old command - filedata = filedata.replace('LOCFILES', locfiles) - nllfile = open(locfile, 'w') - nllfile.write(filedata) - nllfile.close() + filedata = filedata.replace('LOCFILES', locfiles) + nllfile = open(locfile, 'w') + nllfile.write(filedata) + nllfile.close() # locate the event subprocess.call([nlloccall, locfile]) ########################################################## - + print '------------------------------------------'