Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gui/snuffler/snufflings/ampspec.py: 16%
101 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-23 12:35 +0000
1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from ..snuffling import Snuffling, Param, Switch
7import numpy as num
8from matplotlib import cm
11def window(freqs, fc, b):
12 w = num.zeros(len(freqs))
13 if fc == 0.:
14 w[freqs == 0.0] = 1.
15 return w
16 mask = num.logical_and(freqs != 0.0, freqs != fc)
17 T = num.log10(freqs[mask]/fc)*b
18 w[mask] = (num.sin(T)/T)**4
19 w[freqs == fc] = 1.
20 w /= num.sum(w)
21 return w
24class AmpSpec(Snuffling):
26 def help(self):
27 return '''
28<html>
29<head>
30<style type="text/css">
31 body { margin-left:10px };
32</style>
33<body>
34<h1 align="center">Plot Amplitude Spectrum</h1>
35<p>
36When smoothing is activated, a smoothing algorithm is applied as
37proposed
38by Konno and Ohmachi, (1998). </p>
39<p style='font-family:courier'>
40 Konno, K. and Omachi, T. (1998). Ground-motion characteristics
41 estimated from spectral ratio between horizontal and vertical
42 components of microtremor, Bull. Seism. Soc. Am., 88, 228-241
43</p>
44</body>
45</html>
46 '''
48 def setup(self):
49 '''Customization of the snuffling.'''
51 self.set_name('Plot Amplitude Spectrum')
52 self.add_parameter(Switch('Smoothing', 'want_smoothing', False))
53 self.add_parameter(Param('Smoothing band width', 'b', 40., 1., 100.))
54 self.set_live_update(False)
55 self._wins = {}
57 def call(self):
58 '''Main work routine of the snuffling.'''
60 all = []
61 for traces in self.chopper_selected_traces(fallback=True):
62 for trace in traces:
63 all.append(trace)
65 colors = iter(cm.Accent(num.linspace(0., 1., len(all))))
66 if self.want_smoothing:
67 alpha = 0.2
68 additional = 'and smoothing'
69 else:
70 alpha = 1.
71 additional = ''
73 flimits = []
74 alimits = []
75 pblabel = 'Calculating amplitude spectra %s' % additional
76 pb = self.get_viewer().parent().get_progressbars()
77 pb.set_status(pblabel, 0)
78 num_traces = len(all)
79 maxval = float(num_traces)
80 plot_data = []
81 plot_data_supplement = []
82 for i_tr, tr in enumerate(all):
83 val = i_tr/maxval*100.
84 pb.set_status(pblabel, val)
86 tr.ydata = tr.ydata.astype(float)
87 tr.ydata -= tr.ydata.mean()
88 f, a = tr.spectrum()
89 flimits.append(f[1])
90 flimits.append(f[-1])
91 absa = num.abs(a)
92 labsa = num.log(absa)
93 stdabsa = num.std(labsa)
94 meanabsa = num.mean(labsa)
95 lamin, lamax = meanabsa - 5*stdabsa, meanabsa + 5*stdabsa
96 alimits.append(num.exp(lamin))
97 alimits.append(num.exp(lamax))
98 c = next(colors)
99 plot_data.append((f, num.abs(a)))
100 plot_data_supplement.append((c, alpha, '.'.join(tr.nslc_id)))
101 if self.want_smoothing:
102 smoothed = self.konnoohmachi(num.abs(a), f, self.b)
103 plot_data.append((f, num.abs(smoothed)))
104 plot_data_supplement.append((c, 1., '.'.join(tr.nslc_id)))
105 self.get_viewer().update()
107 pb.set_status(pblabel, 100.)
109 fig = self.figure(name='Amplitude Spectra')
110 p = fig.add_subplot(111)
111 args = ('c', 'alpha', 'label')
112 for d, s in zip(plot_data, plot_data_supplement):
113 p.plot(*d, **dict(zip(args, s)))
115 p.set_xscale('log')
116 p.set_yscale('log')
118 amin, amax = num.min(alimits), num.max(alimits)
119 p.set_ylim(amin, amax)
121 fmin, fmax = num.min(flimits), num.max(flimits)
122 p.set_xlim(fmin, fmax)
124 p.set_xlabel('Frequency [Hz]')
125 p.set_ylabel('Counts')
127 handles, labels = p.get_legend_handles_labels()
128 leg_dict = dict(zip(labels, handles))
129 if num_traces > 1:
130 p.legend(list(leg_dict.values()), list(leg_dict.keys()),
131 loc=2,
132 borderaxespad=0.,
133 bbox_to_anchor=((1.05, 1.)))
134 fig.subplots_adjust(right=0.8,
135 left=0.1)
136 else:
137 p.set_title(list(leg_dict.keys())[0], fontsize=16)
138 fig.subplots_adjust(right=0.9,
139 left=0.1)
140 fig.canvas.draw()
142 def konnoohmachi(self, amps, freqs, b=20):
143 smooth = num.zeros(len(freqs), dtype=freqs.dtype)
144 amps = num.array(amps)
145 for i, fc in enumerate(freqs):
146 fkey = tuple((b, fc, freqs[0], freqs[1], freqs[-1]))
147 if fkey in self._wins.keys():
148 win = self._wins[fkey]
149 else:
150 win = window(freqs, fc, b)
151 self._wins[fkey] = win
152 smooth[i] = num.sum(win*amps)
154 return smooth
157def __snufflings__():
158 '''
159 Returns a list of snufflings to be exported by this module.
160 '''
162 return [AmpSpec()]