1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
8from ..snuffling import Param, Snuffling, Switch, Choice
9from ..marker 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 panel_visibility_changed(self, visible):
105 viewer = self.get_viewer()
106 if visible:
107 viewer.pile_has_changed_signal.connect(self.adjust_controls)
108 self.adjust_controls()
110 else:
111 viewer.pile_has_changed_signal.disconnect(self.adjust_controls)
113 def adjust_controls(self):
114 viewer = self.get_viewer()
115 dtmin, dtmax = viewer.content_deltat_range()
116 maxfreq = 0.5/dtmin
117 minfreq = (0.5/dtmax)*0.001
118 self.set_parameter_range('lowpass', minfreq, maxfreq)
119 self.set_parameter_range('highpass', minfreq, maxfreq)
121 def copy_passband(self):
122 viewer = self.get_viewer()
123 self.set_parameter('lowpass', viewer.lowpass)
124 self.set_parameter('highpass', viewer.highpass)
126 def call(self):
127 '''
128 Main work routine of the snuffling.
129 '''
131 self.cleanup()
133 swin, ratio = self.swin, self.ratio
134 lwin = swin * ratio
135 tpad = lwin
137 data_pile = self.get_pile()
139 viewer = self.get_viewer()
140 deltat_min = viewer.content_deltat_range()[0]
142 tinc = max(lwin * 2., 500000. * deltat_min)
144 show_level_traces = self.show_level_traces
146 nsamples_added = [0]
148 def update_sample_count(traces):
149 for tr in traces:
150 nsamples_added[0] += tr.data_len()
152 markers = []
154 for batch in self.chopper_selected_traces(
155 tinc=tinc, tpad=tpad,
156 want_incomplete=False,
157 fallback=True,
158 style='batch',
159 mode='visible',
160 progress='Calculating STA/LTA',
161 responsive=True,
162 marker_selector=lambda marker: marker.tmin != marker.tmax,
163 trace_selector=lambda x: not (x.meta and x.meta.get(
164 'tabu', False))):
166 sumtrace = None
167 isum = 0
168 for tr in batch.traces:
169 if self.lowpass is not None:
170 tr.lowpass(4, self.lowpass, nyquist_exception=True)
172 if self.highpass is not None:
173 tr.highpass(4, self.highpass, nyquist_exception=True)
175 sta_lta = {
176 'centered': tr.sta_lta_centered,
177 'right': tr.sta_lta_right}[self.variant]
179 sta_lta(
180 swin, lwin,
181 scalingmethod=scalingmethod_map[self.scalingmethod])
183 tr.chop(batch.tmin, batch.tmax)
185 if not self.apply_to_sum:
186 markers.extend(trace_to_pmarkers(tr, self.level, swin))
188 tr.set_codes(location='cg')
189 tr.meta = {'tabu': True}
191 if sumtrace is None:
192 ny = int((tr.tmax - tr.tmin) / data_pile.deltatmin)
193 sumtrace = trace.Trace(
194 deltat=data_pile.deltatmin,
195 tmin=tr.tmin,
196 ydata=num.zeros(ny))
198 sumtrace.set_codes(
199 network='', station='SUM',
200 location='cg', channel='')
202 sumtrace.meta = {'tabu': True}
204 sumtrace.add(tr, left=None, right=None)
205 isum += 1
207 if sumtrace is not None:
208 sumtrace.ydata /= float(isum)
209 if self.apply_to_sum:
210 markers.extend(
211 trace_to_pmarkers(sumtrace, self.level, swin,
212 [('*', '*', '*', '*')]))
214 if show_level_traces:
215 update_sample_count([sumtrace])
216 self.add_trace(sumtrace)
218 self.add_markers(markers)
219 markers = []
221 if show_level_traces:
222 update_sample_count(batch.traces)
223 self.add_traces(batch.traces)
225 if show_level_traces and nsamples_added[0] > 10000000:
226 self.error(
227 'Limit reached. Turning off further display of level '
228 'traces to prevent memory exhaustion.')
230 show_level_traces = False
233def trace_to_pmarkers(tr, level, swin, nslc_ids=None):
234 markers = []
235 tpeaks, apeaks = tr.peaks(level, swin)
236 for t, a in zip(tpeaks, apeaks):
237 ids = nslc_ids or [tr.nslc_id]
238 mark = Marker(ids, t, t, )
239 print(mark, a)
240 markers.append(mark)
242 return markers
245def __snufflings__():
246 return [DetectorSTALTA()]