1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import time
7import logging
8from collections import defaultdict
10import numpy as num
12from ..snuffling import Snuffling, Param, Switch
13from pyrocko.trace import Trace, NoData
14from pyrocko import util
16logger = logging.getLogger('pyrocko.gui.snuffling.rms')
19class RootMeanSquareSnuffling(Snuffling):
21 '''
22 Create traces with blockwise root mean square values.
23 '''
25 def setup(self):
26 '''
27 Customization of the snuffling.
28 '''
30 self.set_name('RMS')
32 self.add_parameter(Param(
33 'Highpass [Hz]', 'highpass', None, 0.001, 1000.,
34 low_is_none=True))
36 self.add_parameter(Param(
37 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000.,
38 high_is_none=True))
40 self.add_parameter(Param(
41 'Block Length [s]', 'block_length', 100., 0.1, 3600.))
43 self.add_parameter(Switch(
44 'Group channels', 'group_channels', True))
46 self.add_parameter(Switch(
47 'Log RMS', 'log', False))
49 self.add_trigger(
50 'Copy passband from Main', self.copy_passband)
52 self.set_live_update(False)
54 def copy_passband(self):
55 viewer = self.get_viewer()
56 self.set_parameter('lowpass', viewer.lowpass)
57 self.set_parameter('highpass', viewer.highpass)
59 def call(self):
60 '''
61 Main work routine of the snuffling.
62 '''
64 self.cleanup()
66 tinc = self.block_length
68 tpad = 0.0
69 for freq in (self.highpass, self.lowpass):
70 if freq is not None:
71 tpad = max(tpad, 1./freq)
73 targets = {}
74 tlast = time.time()
75 for batch in self.chopper_selected_traces(
76 tinc=tinc,
77 tpad=tpad,
78 want_incomplete=False,
79 style='batch',
80 mode='visible',
81 progress='Calculating RMS',
82 responsive=True,
83 fallback=True):
85 tcur = batch.tmin + 0.5 * tinc
87 if self.group_channels:
88 def grouping(nslc):
89 return nslc[:3] + (nslc[3][:-1],)
91 else:
92 def grouping(nslc):
93 return nslc
95 rms = defaultdict(list)
96 for tr in batch.traces:
98 # don't work traces produced by this (and other) snuffling
99 if tr.meta and 'tabu' in tr.meta and tr.meta['tabu']:
100 continue
102 if self.lowpass is not None:
103 tr.lowpass(4, self.lowpass, nyquist_exception=True)
105 if self.highpass is not None:
106 tr.highpass(4, self.highpass, nyquist_exception=True)
108 try:
109 tr.chop(batch.tmin, batch.tmax)
110 val = 0.0
111 if tr.ydata.size != 0:
112 y = tr.ydata.astype(float)
113 if num.any(num.isnan(y)):
114 logger.error(
115 'NaN value in trace %s (%s - %s)' % (
116 '.'.join(tr.nslc_id),
117 util.time_to_str(tr.tmin),
118 util.time_to_str(tr.tmax)))
120 val = num.sqrt(num.sum(y**2)/tr.ydata.size)
122 rms[grouping(tr.nslc_id)].append(val)
124 except NoData:
125 continue
127 tnow = time.time()
129 insert_now = False
130 if tnow > tlast + 0.8:
131 insert_now = True
132 tlast = tnow
134 to_add = []
135 for key, values in rms.items():
136 target = targets.get(key, None)
138 if not target \
139 or abs((target.tmax + tinc) - tcur) > 0.01 * tinc \
140 or insert_now:
142 if target:
143 to_add.append(target)
145 target = targets[key] = Trace(
146 network=key[0],
147 station=key[1],
148 location=key[2],
149 channel=key[3] + '-RMS',
150 tmin=tcur,
151 deltat=tinc,
152 ydata=num.zeros(0, dtype=float),
153 meta={'tabu': True})
155 value = num.atleast_1d(
156 num.sqrt(num.sum(num.array(values, dtype=float)**2)))
158 if self.log and value != 0.0:
159 value = num.log(value)
161 targets[key].append(value)
163 if to_add:
164 self.add_traces(to_add)
166 # add the newly created traces to the viewer
167 if targets.values():
168 self.add_traces(list(targets.values()))
171def __snufflings__():
172 '''
173 Returns a list of snufflings to be exported by this module.
174 '''
176 return [RootMeanSquareSnuffling()]