Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gui/snuffler/snufflings/stalta.py: 28%

103 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import numpy as num 

7 

8from ..snuffling import Param, Snuffling, Switch, Choice 

9from ..marker import Marker 

10from pyrocko import trace 

11 

12h = 3600. 

13 

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)]) 

16 

17 

18class DetectorSTALTA(Snuffling): 

19 

20 def help(self): 

21 return ''' 

22<html> 

23<head> 

24<style type="text/css"> 

25 body { margin-left:10px }; 

26</style> 

27</head> 

28<body> 

29<h1 align="center">STA/LTA</h1> 

30<p> 

31Detect onsets automatically using the Short-Time-Average/Long-Time-Average 

32ratio.<br/> 

33This snuffling uses the method: 

34<a href="http://emolch.github.io/pyrocko/v0.3/trace.html#pyrocko.trace.\ 

35Trace.sta_lta_centered" style="text-decoration:none"> 

36 <pre>pyrocko.trace.Trace.sta_lta_centered</pre> 

37</a></p> 

38<p> 

39<b>Parameters:</b><br /> 

40 <b>&middot; Highpass [Hz]</b> 

41 - Apply high pass filter before analysing.<br /> 

42 <b>&middot; Lowpass [Hz]</b> 

43 - Apply low pass filter before analysing.<br /> 

44 <b>&middot; Short Window [s]</b> 

45 - Window length of the short window.<br /> 

46 <b>&middot; Ratio</b> 

47 - Long window length is the short window length times the 

48 <b>Ratio</b>.<br /> 

49 <b>&middot; Level</b> 

50 - Define a trigger threshold. A marker is added where STA/LTA 

51 ratios exceed this threshold. <br /> 

52 <b>&middot; Processing Block length</b> 

53 - Subdivide dataset in blocks for analysis. <br /> 

54 <b>&middot; Show trigger level traces </b> 

55 - Add level traces showing the STA/LTA ration for each 

56 trace.<br /> 

57 <b>&middot; Apply to full dataset</b> 

58 - If marked entire loaded dataset will be analyzed. <br /> 

59 <b>&middot; Scaling/Normalization method</b> 

60 - Select how output of the STA/LTA should be scaled.</ br> 

61</p> 

62<p> 

63A helpfull description of how to tune the STA/LTA's parameters can be found 

64in the the following ebook chapter by Amadej Trnkoczy: <a 

65href="http://ebooks.gfz-potsdam.de/pubman/item/escidoc:4097:3/component/\ 

66escidoc:4098/IS_8.1_rev1.pdf">Understanding 

67and parameter setting of STA/LTA trigger algorithm</a> 

68</p> 

69</body> 

70</html> 

71''' 

72 

73 def setup(self): 

74 

75 self.set_name('STA LTA Detector') 

76 

77 self.add_parameter(Param( 

78 'Highpass [Hz]', 'highpass', None, 0.001, 1000., 

79 low_is_none=True)) 

80 

81 self.add_parameter(Param( 

82 'Lowpass [Hz]', 'lowpass', None, 0.001, 1000., 

83 high_is_none=True)) 

84 

85 self.add_parameter(Param( 

86 'Short Window [s]', 'swin', 30., 0.01, 2*h)) 

87 self.add_parameter(Param( 

88 'Ratio', 'ratio', 3., 1.1, 20.)) 

89 self.add_parameter(Param( 

90 'Level', 'level', 0.5, 0., 1.)) 

91 self.add_parameter(Switch( 

92 'Show trigger level traces', 'show_level_traces', False)) 

93 self.add_parameter(Choice( 

94 'Variant', 'variant', 'centered', ['centered', 'right'])) 

95 self.add_parameter(Choice( 

96 'Scaling/Normalization method', 'scalingmethod', '[0-1]', 

97 scalingmethods)) 

98 self.add_parameter( 

99 Switch('Detect on sum trace', 'apply_to_sum', False)) 

100 

101 self.add_trigger( 

102 'Copy passband from Main', self.copy_passband) 

103 

104 self.set_live_update(False) 

105 

106 def panel_visibility_changed(self, visible): 

107 viewer = self.get_viewer() 

108 if visible: 

109 viewer.pile_has_changed_signal.connect(self.adjust_controls) 

110 self.adjust_controls() 

111 

112 else: 

113 viewer.pile_has_changed_signal.disconnect(self.adjust_controls) 

114 

115 def adjust_controls(self): 

116 viewer = self.get_viewer() 

117 dtmin, dtmax = viewer.content_deltat_range() 

118 maxfreq = 0.5/dtmin 

119 minfreq = (0.5/dtmax)*0.001 

120 self.set_parameter_range('lowpass', minfreq, maxfreq) 

121 self.set_parameter_range('highpass', minfreq, maxfreq) 

122 

123 def copy_passband(self): 

124 viewer = self.get_viewer() 

125 self.set_parameter('lowpass', viewer.lowpass) 

126 self.set_parameter('highpass', viewer.highpass) 

127 

128 def call(self): 

129 ''' 

130 Main work routine of the snuffling. 

131 ''' 

132 

133 self.cleanup() 

134 

135 swin, ratio = self.swin, self.ratio 

136 lwin = swin * ratio 

137 tpad = lwin 

138 

139 data_pile = self.get_pile() 

140 

141 viewer = self.get_viewer() 

142 deltat_min = viewer.content_deltat_range()[0] 

143 

144 tinc = max(lwin * 2., 500000. * deltat_min) 

145 

146 show_level_traces = self.show_level_traces 

147 

148 nsamples_added = [0] 

149 

150 def update_sample_count(traces): 

151 for tr in traces: 

152 nsamples_added[0] += tr.data_len() 

153 

154 markers = [] 

155 

156 for batch in self.chopper_selected_traces( 

157 tinc=tinc, tpad=tpad, 

158 want_incomplete=False, 

159 fallback=True, 

160 style='batch', 

161 mode='visible', 

162 progress='Calculating STA/LTA', 

163 responsive=True, 

164 marker_selector=lambda marker: marker.tmin != marker.tmax, 

165 trace_selector=lambda x: not (x.meta and x.meta.get( 

166 'tabu', False))): 

167 

168 sumtrace = None 

169 isum = 0 

170 for tr in batch.traces: 

171 if self.lowpass is not None: 

172 tr.lowpass(4, self.lowpass, nyquist_exception=True) 

173 

174 if self.highpass is not None: 

175 tr.highpass(4, self.highpass, nyquist_exception=True) 

176 

177 sta_lta = { 

178 'centered': tr.sta_lta_centered, 

179 'right': tr.sta_lta_right}[self.variant] 

180 

181 sta_lta( 

182 swin, lwin, 

183 scalingmethod=scalingmethod_map[self.scalingmethod]) 

184 

185 tr.chop(batch.tmin, batch.tmax) 

186 

187 if not self.apply_to_sum: 

188 markers.extend(trace_to_pmarkers(tr, self.level, swin)) 

189 

190 tr.set_codes(location='cg') 

191 tr.meta = {'tabu': True} 

192 

193 if sumtrace is None: 

194 ny = int((tr.tmax - tr.tmin) / data_pile.deltatmin) 

195 sumtrace = trace.Trace( 

196 deltat=data_pile.deltatmin, 

197 tmin=tr.tmin, 

198 ydata=num.zeros(ny)) 

199 

200 sumtrace.set_codes( 

201 network='', station='SUM', 

202 location='cg', channel='') 

203 

204 sumtrace.meta = {'tabu': True} 

205 

206 sumtrace.add(tr, left=None, right=None) 

207 isum += 1 

208 

209 if sumtrace is not None: 

210 sumtrace.ydata /= float(isum) 

211 if self.apply_to_sum: 

212 markers.extend( 

213 trace_to_pmarkers(sumtrace, self.level, swin, 

214 [('*', '*', '*', '*')])) 

215 

216 if show_level_traces: 

217 update_sample_count([sumtrace]) 

218 self.add_trace(sumtrace) 

219 

220 self.add_markers(markers) 

221 markers = [] 

222 

223 if show_level_traces: 

224 update_sample_count(batch.traces) 

225 self.add_traces(batch.traces) 

226 

227 if show_level_traces and nsamples_added[0] > 10000000: 

228 self.error( 

229 'Limit reached. Turning off further display of level ' 

230 'traces to prevent memory exhaustion.') 

231 

232 show_level_traces = False 

233 

234 

235def trace_to_pmarkers(tr, level, swin, nslc_ids=None): 

236 markers = [] 

237 tpeaks, apeaks = tr.peaks(level, swin) 

238 for t, a in zip(tpeaks, apeaks): 

239 ids = nslc_ids or [tr.nslc_id] 

240 mark = Marker(ids, t, t, ) 

241 print(mark, a) 

242 markers.append(mark) 

243 

244 return markers 

245 

246 

247def __snufflings__(): 

248 return [DetectorSTALTA()]