ObsPy Trigger/Picker Tutorial

This is a small tutorial for the UNESCO short course on triggering. Test data used in this tutorial can be downloaded here:  trigger_data.zip. The triggers are implemented as described in Withers et al. (1998). There are two scripts in the repository which we used for detection in the Bavarian network. They might be useful as reference code examples or serve as cookbook for similar applications.
Please also note the convenience method of ObsPy's  Stream and  Trace objects for triggering.

• branches/sandbox/stalta/stalta4baynet.py
• branches/sandbox/stalta/coincidence4baynet.py

Reading files into Python variables

The data files are read into an ObsPy Trace object using read().

>>> from obspy.core import read
>>> st = read("data/ev0_6.a01.gse2")
>>> st = st.select(component="Z")
>>> tr = st[0]

For more details on Trace and Stream objects, have a look at the ObsPy Tutorial.

The data format is automatically detected. Important in this tutorial are the Trace attributes:

• tr.data: contains the data as numpy.ndarray
• tr.stats: contains a dictionary-like class of header entries
• tr.stats.sampling_rate: the sampling rate
• tr.stats.npts: sample count of data

As an example, the header of the data file is printed and the data are plotted like this (for simpler preview plots use tr.plot()):

>>> import matplotlib.pylab as plt
>>> plt.ion() # turn on interactive mode
>>> plt.plot(tr.data,'k')
>>> print tr.stats
         network: 
         station: EV0_6
        location: 
         channel: EHZ
       starttime: 1970-01-01T01:00:00.000000Z
         endtime: 1970-01-01T01:00:59.995000Z
   sampling_rate: 200.0
           delta: 0.005
            npts: 12000
           calib: 1.0
         _format: GSE2
            gse2: AttribDict({'instype': '      ', 'datatype': 'CM6', 'hang': 0.0, 'auxid': '    ', 'vang': -1.0, 'calper': 1.0})

Available trigger methods

After loading the data, we are able to pass the waveform data to the following trigger routines defined in obspy.signal.trigger:

def recStalta(a, nsta, nlta):
def recStaltaPy(a, nsta, nlta):
def carlStaTrig(a, Nsta, Nlta, ratio, quiet):
def classicStaLta(a, Nsta, Nlta):
def delayedStaLta(a, Nsta, Nlta):
def zdetect(a, Nsta):
def pkBaer(reltrc,samp_int,tdownmax,tupevent,thr1,thr2,preset_len,p_dur):
def arPick(a,b,c,samp_rate,f1,f2,lta_p,sta_p,lta_s,sta_s,m_p,m_s,l_p,l_s, s_pick=True):

Help for each function is available  HTML formatted or in the usual Python manner:

from obspy.signal.trigger import *
help(classicStaLta)

The triggering itself mainly consists of the following two steps:

• Calculating the characteristic function
• Setting picks based on values of the characteristic function

Trigger Examples

For all the examples, the commands to read in the data and to load the modules are the following:

import matplotlib.pyplot as plt
from obspy.core import read
from obspy.signal.trigger import *
from obspy.imaging.waveform import plot_trigger
import numpy as np

trace = read("data/ev0_6.a01.gse2")[0]
df = trace.stats.sampling_rate

(plot_trigger is also part of the already mentioned zip file  trigger_data.zip and can be imported from the included file.)

Classic Sta Lta

cft = classicStaLta(trace.data, 5.*df, 10.*df)
plot_trigger(trace, cft, 1.5, 0.5)

Z-Detect

cft = zdetect(trace.data, 10.*df)
plot_trigger(trace, cft, -0.4, -0.3)

Recursive Sta Lta

cft = recStaltaPy(trace.data, 5.*df, 10.*df)
plot_trigger(trace, cft, 1.2, 0.5)

Carl-Sta-Trig

cft = carlStaTrig(trace.data, 5.*df, 10.*df, 0.8, 0.8)
plot_trigger(trace, cft, 20.0, -20.0)

Delayed Sta Lta

cft = delayedStaLta(trace.data, 5.*df, 10.*df)
plot_trigger(trace, cft, 5, 10)

Baer Picker

For pkBaer, input is in seconds, output is in samples

from obspy.signal.trigger import pkBaer
ntdownmax, ntupevent, thr1, thr2, npreset_len, np_dur = (20, 60, 7.0, 12.0, 100, 100)
ppick, phase_info = pkBaer(tr.data.astype('float32'), df, ntdownmax, ntupevent, thr1, thr2, npreset_len, np_dur)
print ppick/df, phase_info

This yields the output 34.465000000000003 EPU3, which means that a P pick was set at 34.46s with Phase information EPU3.

AR Picker

For arPick, input and output are in seconds.

from obspy.core import read
from obspy.signal.trigger import arPick
tr1 = read('data/loc_RJOB20050801145719850.z.gse2')[0]
tr2 = read('data/loc_RJOB20050801145719850.n.gse2')[0]
tr3 = read('data/loc_RJOB20050801145719850.e.gse2')[0]
a, b, c = tr1.data.astype('float32'), tr2.data.astype('float32'), tr3.data.astype('float32')
df = tr1.stats.sampling_rate
f1, f2, lta_p, sta_p, lta_s, sta_s, m_p, m_s, l_p, l_s = 1.0, 20.0, 1.0, 0.1, 4.0, 1.0, 2, 8, 0.1, 0.2
ppick, spick = arPick(a, b, c, df, f1, f2, lta_p, sta_p, lta_s, sta_s, m_p, m_s, l_p, l_s)
print ppick, spick

This gives the output 30.6350002289 31.2800006866, meaning that a P pick at 30.63s and an S pick at 31.28s were identified.

Advanced Example

A more complicated example, where the data are retrieved via arclink and results are plotted step by step, is shown here:

from obspy.core import UTCDateTime
from obspy.arclink import Client
from obspy.signal.trigger import recStaltaPy, triggerOnset
import matplotlib.pyplot as plt
import numpy as np

# Retrieve waveforms via ArcLink
client = Client(host="webdc.eu", port=18001)
t = UTCDateTime("2009-08-24 00:19:45")
st = client.getWaveform('BW', 'RTSH', '', 'EHZ', t, t+50)

# For convenience
tr = st[0] #only one trace in mseed volume
df = tr.stats.sampling_rate

# Characteristic function and trigger onsets
cft = recStaltaPy(tr.data,2.5*df,10.*df)
on_of = triggerOnset(cft,3.5,0.5)

# Plotting the results
plt.ion()
ax = plt.subplot(211)
plt.plot(tr.data,'k')
ymin, ymax = ax.get_ylim()
plt.vlines(on_of[:,0],ymin, ymax, color='r',linewidth=2)
plt.vlines(on_of[:,1],ymin, ymax, color='b',linewidth=2)
plt.subplot(212, sharex=ax)
plt.plot(cft, 'k')
plt.hlines([3.5, 0.5], 0, len(cft), color=['r','b'], linestyle='--')