diff --git a/autoPyLoT.py b/autoPyLoT.py index 06a3bdef..1d6fb6cf 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,20 +164,24 @@ 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]) ########################################################## +<<<<<<< HEAD # write phase files for various location routines # HYPO71 hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID')) writephases(picks, 'HYPO71', hypo71file) +======= + +>>>>>>> 223902f2d4a9fe992730bb1ba925bfa432f26550 print '------------------------------------------' print '-------Finished event %s!-------' % parameter.getParam('eventID') diff --git a/autoPyLoT_local.in b/autoPyLoT_local.in index ae4556c0..baabf0ce 100644 --- a/autoPyLoT_local.in +++ b/autoPyLoT_local.in @@ -9,15 +9,15 @@ EVENT_DATA/LOCAL #datapath# %data path 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 +PILOT #datastructure#%choose data structure +0 #iplot# %flag for plotting: 0 none, 1 partly, >1 everything %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -#NLLoc settings +#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.in #locfile# %name of autoPyLoT-output control file for NLLoc +Insheim_min1d2015_auto.in #locfile# %name of autoPyLoT-output control file for NLLoc %(in nllocroot/run) ttime #ttpatter# %pattern of NLLoc ttimes from grid %(in nllocroot/times) diff --git a/autoPyLoT_regional.in b/autoPyLoT_regional.in index 52f5c465..4248fd30 100644 --- a/autoPyLoT_regional.in +++ b/autoPyLoT_regional.in @@ -88,7 +88,7 @@ ARH #algoS# %choose algorithm for S-onset 2.5 #noisefactor# %noiselevel*noisefactor=threshold 60 #minpercent# %required percentage of samples higher than threshold #check for spuriously picked S-onsets# -1.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude +0.5 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude #check statistics of P onsets# 45 #mdttolerance# %maximum allowed deviation of P picks from median [s] #wadati check# diff --git a/pylot/RELEASE-VERSION b/pylot/RELEASE-VERSION index 01eb3d3d..070ed7aa 100644 --- a/pylot/RELEASE-VERSION +++ b/pylot/RELEASE-VERSION @@ -1 +1 @@ -ac7d-dirty +a31e-dirty diff --git a/pylot/core/active/fmtomo2vtk.py b/pylot/core/active/fmtomo2vtk.py index 95b90705..59c0968e 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 @@ -91,13 +93,39 @@ def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'): 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 - 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) / 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)) + 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)) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 7e7dd09a..58d021be 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() @@ -363,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]. @@ -386,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): @@ -396,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 @@ -413,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 @@ -432,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: @@ -444,19 +483,29 @@ 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): - vel = 0.0 - # elif radius > (R + z - 15): ########### TESTING - # vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * 1.0 + vmin ########################## + topo = point[2] + z = -(R + topo - radius) + if z > (topo + 1): + vel = vtop[0] + 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 - 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 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'): @@ -469,31 +518,36 @@ 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) - 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 = pointsize, label = 'all receivers') + if len(xsc) > 0: + ax.plot(xsc, ysc, 'b*', markersize = pointsize, label = 'shot locations') if plot_topo == True: - plt.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: - plt.plot(xmr, ymr, 'ro', label = 'measured receivers') + ax.plot(xmr, ymr, 'r.', markersize = pointsize, 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,18 +562,25 @@ 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() 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 @@ -545,7 +606,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') 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: diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 19d33b02..5a5cc4e5 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -335,7 +335,8 @@ def autopickstation(wfstream, pickparam): "no zero crossings derived!") print ("Cannot calculate source spectrum!") else: - calcwin = (zc[3] - zc[0]) * z_copy[0].stats.delta + index = min([3, len(zc) - 1]) + calcwin = (zc[index] - zc[0]) * z_copy[0].stats.delta # calculate source spectrum and get w0 and fc specpara = DCfc(z_copy, mpickP, calcwin, iplot) w0 = specpara.getw0()