Marginal changes.

pylot/core/pick/utils.py

@@ -13,8 +13,6 @@ import scipy as sc
 import matplotlib.pyplot as plt
 from obspy.core import Stream, UTCDateTime
 import warnings
-import pdb
-
 
 def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
     '''
@@ -257,6 +255,8 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
             elif P1[0] > 0 and P2[0] <= 0:
                 FM = '+'
            
+        print 'fmpicker: Found polarity %s' % FM
+
     if iplot > 1:
         plt.figure(iplot)
         plt.subplot(2, 1, 1)
@@ -301,48 +301,34 @@ def crossings_nonzero_all(data):
     return ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0]
 
 
-def getSNR(st, TSNR, t0):
+def getSNR(X, TSNR, t1):
-    """
+    '''
-    returns the maximum signal to noise ratio SNR (also in dB) and the
+    Function to calculate SNR of certain part of seismogram relative to
-    corresponding noise level for a given data stream ST ,initial time T0 and
+    given time (onset) out of given noise and signal windows. A safety gap
-    time window parameter tuple TSNR
+    between noise and signal part can be set. Returns SNR and SNR [dB] and
+    noiselevel.
 
-    :param: st, time series (seismogram)
+    :param: X, time series (seismogram)
     :type:  `~obspy.core.stream.Stream`
-    :param: TSNR, length of time windows [s] around t0 (onset) used to determine
+
-     SNR
+    :param: TSNR, length of time windows [s] around t1 (onset) used to determine SNR
     :type: tuple (T_noise, T_gap, T_signal)
-    :param: t0, initial time (onset) from which noise and signal windows are calculated
+
+    :param: t1, initial time (onset) from which noise and signal windows are calculated
     :type: float
-    :return: SNR, SNRdB, noiselevel
+    '''
 
-    ..examples:
+    assert isinstance(X, Stream), "%s is not a stream object" % str(X)
 
-    >>> from obspy.core import read
+    x = X[0].data
-    >>> st = read()
+    t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
-    >>> result = getSNR(st, (6., .3, 3.), 4.67)
+                  X[0].stats.delta)
-    >>> print result
-    (5.1267717641040758, 7.0984398375666435, 132.89370192191919)
-    >>> result = getSNR(st, (8., .2, 5.), 4.67)
-    >>> print result
-    (4.645441835797703, 6.6702702677384131, 133.03562794665109)
-    """
-
-    assert isinstance(st, Stream), "%s is not a stream object" % str(st)
-
-    SNR = None
-    noiselevel = None
-
-    for tr in st:
-        x = tr.data
-        t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate,
-                      tr.stats.delta)
 
     # get noise window
-        inoise = getnoisewin(t, t0, TSNR[0], TSNR[1])
+    inoise = getnoisewin(t, t1, TSNR[0], TSNR[1])
 
     # get signal window
-        isignal = getsignalwin(t, t0, TSNR[2])
+    isignal = getsignalwin(t, t1, TSNR[2])
     if np.size(inoise) < 1:
         print 'getSNR: Empty array inoise, check noise window!'
         return
@@ -354,18 +340,9 @@ def getSNR(st, TSNR, t0):
     x = x - np.mean(x[inoise])
 
     # calculate ratios
-        new_noiselevel = np.sqrt(np.mean(np.square(x[inoise])))
+    noiselevel = np.sqrt(np.mean(np.square(x[inoise])))
     signallevel = np.sqrt(np.mean(np.square(x[isignal])))
-        newSNR = signallevel / new_noiselevel
+    SNR = signallevel / noiselevel
-        if not SNR or newSNR > SNR:
-            SNR = newSNR
-            noiselevel = new_noiselevel
-
-    if not SNR or not noiselevel:
-        raise ValueError('signal to noise ratio could not be calculated:\n'
-                         'noiselevel: {0}\n'
-                         'SNR: {1}'.format(noiselevel, SNR))
-
     SNRdB = 10 * np.log10(SNR)
 
     return SNR, SNRdB, noiselevel
@@ -496,6 +473,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
            Spicks.append(UTCSpick.timestamp)
            SPtimes.append(spt)
 
+
     if len(SPtimes) >= 3:
     	# calculate slope
     	p1 = np.polyfit(Ppicks, SPtimes, 1)
@@ -527,8 +505,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
                     checkedPpicks.append(checkedPpick.timestamp)
                     checkedSpick = UTCDateTime(pickdic[key]['S']['mpp'])
                     checkedSpicks.append(checkedSpick.timestamp)
-                    checkedSPtime = pickdic[key]['S']['mpp'] - \
+                    checkedSPtime = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp']
-                                    pickdic[key]['P']['mpp']
                     checkedSPtimes.append(checkedSPtime)
 
                 pickdic[key]['S']['marked'] = marker
@@ -563,8 +540,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
                 plt.title('Wadati-Diagram, %d S-P Times, Vp/Vs(raw)=%5.2f,' \
                           'Vp/Vs(checked)=%5.2f' % (len(SPtimes), vpvsr, cvpvsr))
                 plt.legend([f1, f2, f3, f4], ['Skipped S-Picks', 'Wadati 1', \
-                                          'Reliable S-Picks', 'Wadati 2'],
+                           'Reliable S-Picks', 'Wadati 2'], loc='best')
-                       loc='best')
         else:
         	plt.title('Wadati-Diagram, %d S-P Times' % len(SPtimes))
 
@@ -649,8 +625,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
                         [minsiglevel, minsiglevel], 'g')
         p5, = plt.plot([pick, pick], [min(x), max(x)], 'b', linewidth=2)
         plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \
-                                          'Envelope Signal Window',
+                    'Envelope Signal Window', 'Minimum Signal Level', \
-                                          'Minimum Signal Level', \
                     'Onset'], loc='best')
         plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
         plt.ylabel('Counts')
@@ -713,8 +688,7 @@ def checkPonsets(pickdic, dttolerance, iplot):
     badstations = np.array(stations)[ibad]
 
     print 'checkPonset: Skipped %d P onsets out of %d' % (len(badstations) \
-                                                          + len(badjkstations),
+            + len(badjkstations), len(stations))
-                                                          len(stations))
 
     goodmarker = 'goodPonsetcheck'
     badmarker = 'badPonsetcheck'
@@ -819,7 +793,108 @@ def jackknife(X, phi, h):
     return PHI_jack, PHI_pseudo, PHI_sub
 
 
+def checkZ4S(X, pick, zfac, checkwin, iplot):
+    '''
+    Function to compare energy content of vertical trace with 
+    energy content of horizontal traces to detect spuriously 
+    picked S onsets instead of P onsets. Usually, P coda shows
+    larger longitudal energy on vertical trace than on horizontal
+    traces, where the transversal energy is larger within S coda.
+    Be careful: there are special circumstances, where this is not 
+    the case!
+
+    : param: X, fitered(!) time series, three traces
+    : type:  `~obspy.core.stream.Stream` 
+
+    : param: pick, initial (AIC) P onset time
+    : type:  float
+    
+    : param: zfac, factor for threshold determination, 
+             vertical energy must exceed coda level times zfac
+             to declare a pick as P onset
+    : type:  float
+
+    : param: checkwin, window length [s] for calculating P-coda
+             energy content
+    : type:  float
+
+    : param: iplot, if iplot > 1, energy content and threshold
+             are shown
+    : type:  int
+    '''
+
+    assert isinstance(X, Stream), "%s is not a stream object" % str(X)
+
+    print 'Check for spuriously picked S onset instead of P onset ...'
+
+    returnflag = 0
+
+    # split components
+    zdat = X.select(component="Z")
+    edat = X.select(component="E")
+    if len(edat) == 0:  # check for other components
+        edat = X.select(component="2")
+    ndat = X.select(component="N")
+    if len(ndat) == 0:  # check for other components
+        ndat = X.select(component="1")
+    
+
+    z = zdat[0].data
+    tz = np.arange(0, zdat[0].stats.npts / zdat[0].stats.sampling_rate,
+                  zdat[0].stats.delta)
+
+    # calculate RMS trace from vertical component
+    absz = np.sqrt(np.power(z, 2))
+    # calculate RMS trace from both horizontal traces
+    # make sure, both traces have equal lengths
+    lene = len(edat[0].data)
+    lenn = len(ndat[0].data)
+    minlen = min([lene, lenn])
+    absen = np.sqrt(np.power(edat[0].data[0:minlen - 1], 2) \
+                    + np.power(ndat[0].data[0:minlen - 1], 2))
+
+    # get signal window
+    isignal = getsignalwin(tz, pick, checkwin)
+
+    # calculate energy levels
+    zcodalevel = max(absz[isignal])
+    encodalevel = max(absen[isignal])
+    
+    # calculate threshold
+    minsiglevel = encodalevel * zfac
+
+    # vertical P-coda level must exceed horizontal P-coda level
+    # zfac times encodalevel
+    if zcodalevel < minsiglevel:
+    	print 'checkZ4S: Maybe S onset? Skip this P pick!'
+    else:
+        print 'checkZ4S: P onset passed checkZ4S test!'
+        returnflag = 1      
+
+    if iplot > 1:
+    	te = np.arange(0, edat[0].stats.npts / edat[0].stats.sampling_rate,
+                  edat[0].stats.delta)
+    	tn = np.arange(0, ndat[0].stats.npts / ndat[0].stats.sampling_rate,
+                  ndat[0].stats.delta)
+    	plt.plot(tz, z / max(z), 'k')
+        plt.plot(tz[isignal], z[isignal] / max(z), 'r')
+        plt.plot(te, edat[0].data / max(edat[0].data) + 1, 'k')
+        plt.plot(te[isignal], edat[0].data[isignal] / max(edat[0].data) + 1, 'r')
+        plt.plot(tn, ndat[0].data / max(ndat[0].data) + 2, 'k')
+        plt.plot(tn[isignal], ndat[0].data[isignal] / max(ndat[0].data) + 2, 'r')
+        plt.plot([tz[isignal[0]], tz[isignal[len(isignal) - 1]]], \
+                  [minsiglevel / max(z), minsiglevel / max(z)], 'g', \
+                   linewidth=2)
+        plt.xlabel('Time [s] since %s' % zdat[0].stats.starttime)
+        plt.ylabel('Normalized Counts')
+        plt.yticks([0, 1, 2], [zdat[0].stats.channel, edat[0].stats.channel, \
+                                        ndat[0].stats.channel])
+        plt.title('CheckZ4S, Station %s' % zdat[0].stats.station)
+        plt.show()
+        raw_input()
+
+    return returnflag
+
 if __name__ == '__main__':
     import doctest
-
     doctest.testmod()