Merge branch 'develop' of ariadne.geophysik.rub.de:/data/git/pylot into develop

Conflicts:
	autoPyLoT.py

autoPyLoT.py

@@ -32,19 +32,17 @@ def autoPyLoT(inputfile):
     .. rubric:: Example
 
     """
-    print '************************************'
+    splash = '''************************************\n
-    print '*********autoPyLoT starting*********'
+                *********autoPyLoT starting*********\n
-    print 'The Python picking and Location Tool'
+                The Python picking and Location Tool\n
-    print '      Version ', _getVersionString(), '2015'
+                Version {version} 2015\n
-    print ' '
+                \n
-    print 'Authors:'
+                Authors:\n
-    print 'S. Wehling-Benatelli'
+                S. Wehling-Benatelli (Ruhr-Universität Bochum)\n
-    print '   Ruhr-Universität Bochum'
+                L. Küperkoch (BESTEC GmbH, Landau i. d. Pfalz)\n
-    print 'L. Küperkoch'
+                K. Olbert (Christian-Albrechts Universität zu Kiel)\n
-    print '   BESTEC GmbH, Landau (Pfalz)'
+                ***********************************'''.format(version=_getVersionString())
-    print 'K. Olbert'
+    print(splash)
-    print '   Christian-Albrechts Universität Kiel'
-    print '************************************'
 
     # 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('Working on event %s' % event)
-                print data
+                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)
+                    filedata = filedata.replace('LOCFILES', locfiles)
-            	    nllfile = open(locfile, 'w')
+                    nllfile = open(locfile, 'w')
-            	    nllfile.write(filedata)
+                    nllfile.write(filedata)
-            	    nllfile.close()
+                    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')

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
+PILOT				    #datastructure#%choose data structure
-0                                   #iplot#        %flag for plotting: 0 none, 1, partly, >1 everything
+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)

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#

pylot/RELEASE-VERSION

@@ -1 +1 @@
-ac7d-dirty
+a31e-dirty

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...")
+    if absOrRel == 'abs':
-    for velocity in vel:
+        print("Writing velocity values to VTK file...")
-        outfile.writelines('%10f\n' %velocity)
+        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))

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), thetaSN = None,
-                      Rbt = (-62.0, 6.0), vbot = 5.5, filename = 'vgrids.in',
+                      phiWE = None, outfilename = 'vgrids.in',
-                      method = 'linear' ):
+                      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):
+        def readMygridNlayers(filename):
-            z = []
+            infile = open(filename, 'r')
-            for point in surface:
+            nlayers = len(infile.readlines()) / 2
-                z.append(point[2])
+            infile.close()
-            return max(z)
+
+            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)
+        thetaDelta = abs(thetaN - thetaS) / float((nTheta - 1))
-        phiDelta = abs(phiE - phiW) / (nPhi - 1)
+        phiDelta = abs(phiE - phiW) / float((nPhi - 1))
-        rDelta = abs(rbot - rtop) / (nR - 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]
+                            topo = point[2]
-                    if radius > (R + z + 1):
+                    z = -(R + topo - radius)
-                        vel = 0.0
+                    if z > (topo + 1): 
-                    # elif radius > (R + z - 15): ########### TESTING
+                        vel = vtop[0]
-                    #     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 - 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')
+        if len(xrc) > 0:
-        plt.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations')
+            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.title('2D plot of seismic array %s'%self.recfile)
-        plt.ylabel('Y [m]')
+        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')
+        plt.title('3D plot of seismic array %s'%self.recfile)
-        ax.plot(xin, yin, zin, 'k.', markersize = 10, label = 'interpolated receivers')
+        if len(xmt) > 0:
-        ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers')
+            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')

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)
+        sc = ax.scatter(x, y, c = z, s = 30, label = 'picked shots', vmin = tmin, vmax = tmax, cmap = cmap, 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)
+        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:

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()