1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import print_function, absolute_import
6import numpy as num
8from pyrocko.gui.snuffling import Param, Snuffling, Switch, Choice
9from pyrocko.gui.util import Marker
10from pyrocko import trace
12h = 3600.
14scalingmethods = ('[0-1]', '[-1/ratio,1]', '[-1/ratio,1] clipped to [0,1]')
15scalingmethod_map = dict([(m, i+1) for (i, m) in enumerate(scalingmethods)])
18class DetectorSTALTA(Snuffling):
20 '''
21 <html>
22 <head>
23 <style type="text/css">
24 body { margin-left:10px };
25 </style>
26 </head>
27 <body>
28 <h1 align="center">STA/LTA</h1>
29 <p>
30 Detect onsets automatically using the Short-Time-Average/Long-Time-Average
31 ratio.<br/>
32 This snuffling uses the method:
33 <a href="http://emolch.github.io/pyrocko/v0.3/trace.html#pyrocko.trace.\
34Trace.sta_lta_centered" style="text-decoration:none">
35 <pre>pyrocko.trace.Trace.sta_lta_centered</pre>
36 </a></p>
37 <p>
38 <b>Parameters:</b><br />
39 <b>· Highpass [Hz]</b>
40 - Apply high pass filter before analysing.<br />
41 <b>· Lowpass [Hz]</b>
42 - Apply low pass filter before analysing.<br />
43 <b>· Short Window [s]</b>
44 - Window length of the short window.<br />
45 <b>· Ratio</b>
46 - Long window length is the short window length times the
47 <b>Ratio</b>.<br />
48 <b>· Level</b>
49 - Define a trigger threshold. A marker is added where STA/LTA
50 ratios exceed this threshold. <br />
51 <b>· Processing Block length</b>
52 - Subdivide dataset in blocks for analysis. <br />
53 <b>· Show trigger level traces </b>
54 - Add level traces showing the STA/LTA ration for each
55 trace.<br />
56 <b>· Apply to full dataset</b>
57 - If marked entire loaded dataset will be analyzed. <br />
58 <b>· Scaling/Normalization method</b>
59 - Select how output of the STA/LTA should be scaled.</ br>
60 </p>
61 <p>
62 A helpfull description of how to tune the STA/LTA's parameters can be found
63 in the the following ebook chapter by Amadej Trnkoczy: <a
64 href="http://ebooks.gfz-potsdam.de/pubman/item/escidoc:4097:3/component/\
65escidoc:4098/IS_8.1_rev1.pdf">Understanding
66 and parameter setting of STA/LTA trigger algorithm</a>
67 </p>
68 </body>
69 </html>
70 '''
71 def setup(self):
73 self.set_name('STA LTA Detector')
75 self.add_parameter(Param(
76 'Highpass [Hz]', 'highpass', None, 0.001, 1000.,
77 low_is_none=True))
79 self.add_parameter(Param(
80 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000.,
81 high_is_none=True))
83 self.add_parameter(Param(
84 'Short Window [s]', 'swin', 30., 0.01, 2*h))
85 self.add_parameter(Param(
86 'Ratio', 'ratio', 3., 1.1, 20.))
87 self.add_parameter(Param(
88 'Level', 'level', 0.5, 0., 1.))
89 self.add_parameter(Switch(
90 'Show trigger level traces', 'show_level_traces', False))
91 self.add_parameter(Choice(
92 'Variant', 'variant', 'centered', ['centered', 'right']))
93 self.add_parameter(Choice(
94 'Scaling/Normalization method', 'scalingmethod', '[0-1]',
95 scalingmethods))
96 self.add_parameter(
97 Switch('Detect on sum trace', 'apply_to_sum', False))
99 self.add_trigger(
100 'Copy passband from Main', self.copy_passband)
102 self.set_live_update(False)
104 def copy_passband(self):
105 viewer = self.get_viewer()
106 self.set_parameter('lowpass', viewer.lowpass)
107 self.set_parameter('highpass', viewer.highpass)
109 def call(self):
110 '''
111 Main work routine of the snuffling.
112 '''
114 self.cleanup()
116 swin, ratio = self.swin, self.ratio
117 lwin = swin * ratio
118 tpad = lwin
120 data_pile = self.get_pile()
122 viewer = self.get_viewer()
123 deltat_min = viewer.content_deltat_range()[0]
125 tinc = max(lwin * 2., 500000. * deltat_min)
127 show_level_traces = self.show_level_traces
129 nsamples_added = [0]
131 def update_sample_count(traces):
132 for tr in traces:
133 nsamples_added[0] += tr.data_len()
135 markers = []
137 for batch in self.chopper_selected_traces(
138 tinc=tinc, tpad=tpad,
139 want_incomplete=False,
140 fallback=True,
141 style='batch',
142 mode='visible',
143 progress='Calculating STA/LTA',
144 responsive=True,
145 marker_selector=lambda marker: marker.tmin != marker.tmax,
146 trace_selector=lambda x: not (x.meta and x.meta.get(
147 'tabu', False))):
149 sumtrace = None
150 isum = 0
151 for tr in batch.traces:
152 if self.lowpass is not None:
153 tr.lowpass(4, self.lowpass, nyquist_exception=True)
155 if self.highpass is not None:
156 tr.highpass(4, self.highpass, nyquist_exception=True)
158 sta_lta = {
159 'centered': tr.sta_lta_centered,
160 'right': tr.sta_lta_right}[self.variant]
162 sta_lta(
163 swin, lwin,
164 scalingmethod=scalingmethod_map[self.scalingmethod])
166 tr.chop(batch.tmin, batch.tmax)
168 if not self.apply_to_sum:
169 markers.extend(trace_to_pmarkers(tr, self.level, swin))
171 tr.set_codes(location='cg')
172 tr.meta = {'tabu': True}
174 if sumtrace is None:
175 ny = int((tr.tmax - tr.tmin) / data_pile.deltatmin)
176 sumtrace = trace.Trace(
177 deltat=data_pile.deltatmin,
178 tmin=tr.tmin,
179 ydata=num.zeros(ny))
181 sumtrace.set_codes(
182 network='', station='SUM',
183 location='cg', channel='')
185 sumtrace.meta = {'tabu': True}
187 sumtrace.add(tr, left=None, right=None)
188 isum += 1
190 if sumtrace is not None:
191 sumtrace.ydata /= float(isum)
192 if self.apply_to_sum:
193 markers.extend(
194 trace_to_pmarkers(sumtrace, self.level, swin,
195 [('*', '*', '*', '*')]))
197 if show_level_traces:
198 update_sample_count([sumtrace])
199 self.add_trace(sumtrace)
201 self.add_markers(markers)
202 markers = []
204 if show_level_traces:
205 update_sample_count(batch.traces)
206 self.add_traces(batch.traces)
208 if show_level_traces and nsamples_added[0] > 10000000:
209 self.error(
210 'Limit reached. Turning off further display of level '
211 'traces to prevent memory exhaustion.')
213 show_level_traces = False
216def trace_to_pmarkers(tr, level, swin, nslc_ids=None):
217 markers = []
218 tpeaks, apeaks = tr.peaks(level, swin)
219 for t, a in zip(tpeaks, apeaks):
220 ids = nslc_ids or [tr.nslc_id]
221 mark = Marker(ids, t, t, )
222 print(mark, a)
223 markers.append(mark)
225 return markers
228def __snufflings__():
229 return [DetectorSTALTA()]