Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/waveform_phase_ratio/measure.py: 22%

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import numpy as num 

5from pyrocko import gf, trace 

6from pyrocko.guts import Object, Float, StringChoice, List, String 

7from pyrocko.gui import marker 

8from grond.meta import store_t 

9 

10 

11guts_prefix = 'grond' 

12 

13 

14response_wwssn_lp = trace.PoleZeroResponse( 

15 poles=[ 

16 -0.2513 + 0.3351J, 

17 -0.2513 - 0.3351J, 

18 -0.0628 + 0.0304J, 

19 -0.0628 - 0.0304J], 

20 zeros=[0., 0., 0.], 

21 constant=383.6) 

22 

23 

24response_wwssn_lp_2 = trace.PoleZeroResponse( 

25 poles=[ 

26 -0.40180 + 0.08559J, 

27 -0.40180 - 0.08559J, 

28 -0.04841, 

29 -0.08816], 

30 zeros=[0., 0., 0.]) 

31 

32response_wa = trace.PoleZeroResponse( 

33 poles=[ 

34 -5.49779 - 5.60886J, 

35 -5.49779 + 5.60886J], 

36 zeros=[0., 0.]) 

37 

38 

39class NamedResponse(StringChoice): 

40 choices = [ 

41 'wood-anderson', 

42 'wwssn-lp'] 

43 

44 map = { 

45 'wood-anderson': response_wa, 

46 'wwssn-lp': response_wwssn_lp} 

47 

48 

49class BruneResponse(trace.FrequencyResponse): 

50 

51 duration = Float.T() 

52 

53 def evaluate(self, freqs): 

54 return 1.0 / (1.0 + (freqs*self.duration)**2) 

55 

56# from pyrocko.plot import response as pr 

57# pr.plot([response_wwssn_lp, response_wwssn_lp_2, response_wa]) 

58 

59 

60class FeatureMethod(StringChoice): 

61 choices = [ 

62 'peak_component', 

63 'peak_to_peak_component', 

64 'peak_absolute_vector', 

65 'spectral_average'] 

66 

67 

68class WaveformQuantity(StringChoice): 

69 choices = [ 

70 'displacement', 

71 'velocity', 

72 'acceleration'] 

73 

74 

75class FeatureMeasurementFailed(Exception): 

76 pass 

77 

78 

79class FeatureMeasure(Object): 

80 name = String.T() 

81 timing_tmin = gf.Timing.T(default='vel:8') 

82 timing_tmax = gf.Timing.T(default='vel:2') 

83 fmin = Float.T(optional=True) 

84 fmax = Float.T(optional=True) 

85 response = trace.FrequencyResponse.T(optional=True) 

86 named_response = NamedResponse.T(optional=True) 

87 channels = List.T(String.T()) 

88 quantity = WaveformQuantity.T(default='displacement') 

89 method = FeatureMethod.T(default='peak_component') 

90 

91 def get_nmodelling_targets(self): 

92 return len(self.channels) 

93 

94 def get_modelling_targets( 

95 self, codes, lat, lon, depth, store_id, backazimuth): 

96 

97 mtargets = [] 

98 for channel in self.channels: 

99 target = gf.Target( 

100 quantity='displacement', 

101 codes=codes + (channel, ), 

102 lat=lat, 

103 lon=lon, 

104 depth=depth, 

105 store_id=store_id) 

106 

107 if channel == 'R': 

108 target.azimuth = backazimuth - 180. 

109 target.dip = 0. 

110 elif channel == 'T': 

111 target.azimuth = backazimuth - 90. 

112 target.dip = 0. 

113 elif channel == 'Z': 

114 target.azimuth = 0. 

115 target.dip = -90. 

116 

117 mtargets.append(target) 

118 

119 return mtargets 

120 

121 def evaluate( 

122 self, engine, source, targets, 

123 dataset=None, 

124 trs=None, 

125 extra_responses=[], 

126 debug=False): 

127 

128 from ..waveform import target as base 

129 

130 trs_processed = [] 

131 trs_orig = [] 

132 for itarget, target in enumerate(targets): 

133 if target.codes[-1] not in self.channels: 

134 continue 

135 

136 store = engine.get_store(target.store_id) 

137 

138 tmin = source.time + store_t( 

139 store, self.timing_tmin, source, target) 

140 tmax = source.time + store_t( 

141 store, self.timing_tmax, source, target) 

142 

143 if self.fmin is not None and self.fmax is not None: 

144 freqlimits = [ 

145 self.fmin/2., 

146 self.fmin, 

147 self.fmax, 

148 self.fmax*2.] 

149 tfade = 1./self.fmin 

150 

151 else: 

152 freqlimits = None 

153 tfade = 0.0 

154 

155 if dataset is not None: 

156 bazi = base.backazimuth_for_waveform( 

157 target.azimuth, target.codes) 

158 

159 tr = dataset.get_waveform( 

160 target.codes, 

161 tinc_cache=1.0/self.fmin, 

162 quantity=self.quantity, 

163 tmin=tmin, 

164 tmax=tmax, 

165 freqlimits=freqlimits, 

166 tfade=tfade, 

167 deltat=store.config.deltat, 

168 cache=True, 

169 backazimuth=bazi) 

170 

171 else: 

172 tr = trs[itarget] 

173 

174 tr.extend( 

175 tmin - tfade, 

176 tmax + tfade, 

177 fillmethod='repeat') 

178 

179 tr = tr.transfer( 

180 freqlimits=freqlimits, 

181 tfade=tfade) 

182 

183 trs_orig.append(tr) 

184 

185 tr = tr.copy() 

186 

187 responses = [] 

188 responses.extend(extra_responses) 

189 

190 ndiff = \ 

191 WaveformQuantity.choices.index(self.quantity) - \ 

192 WaveformQuantity.choices.index(target.quantity) 

193 

194 if ndiff > 0: 

195 responses.append(trace.DifferentiationResponse(ndiff)) 

196 

197 if ndiff < 0: 

198 responses.append(trace.IntegrationResponse(-ndiff)) 

199 

200 if self.response: 

201 responses.append(self.response) 

202 

203 if self.named_response: 

204 responses.append( 

205 NamedResponse.map[self.named_response]) 

206 

207 if responses: 

208 trans = trace.MultiplyResponse(responses) 

209 try: 

210 tr = tr.transfer(transfer_function=trans) 

211 

212 except trace.TraceTooShort: 

213 raise FeatureMeasurementFailed( 

214 'transfer: trace too short') 

215 

216 if tmin is None or tmax is None: 

217 raise FeatureMeasurementFailed( 

218 'timing determination failed (phase unavailable?)') 

219 

220 tr.chop(tmin, tmax) 

221 

222 tr.set_location(tr.location + '-' + self.name + '-proc') 

223 trs_processed.append(tr) 

224 

225 markers = [] 

226 marker_candidates = [] 

227 if self.method in ['peak_component', 'peak_to_peak_component']: 

228 component_amp_maxs = [] 

229 for tr in trs_processed: 

230 y = tr.get_ydata() 

231 if self.method == 'peak_component': 

232 yabs = num.abs(y) 

233 i_at_amax = num.argmax(yabs) 

234 amax = yabs[i_at_amax] 

235 if debug: 

236 t_at_amax = tr.tmin + i_at_amax * tr.deltat 

237 mark = marker.Marker( 

238 [tr.nslc_id], 

239 t_at_amax, 

240 t_at_amax, 

241 0) 

242 

243 marker_candidates.append(mark) 

244 

245 component_amp_maxs.append(amax) 

246 else: 

247 i_at_amax = num.argmax(y) 

248 i_at_amin = num.argmin(y) 

249 amax = y[i_at_amax] 

250 amin = y[i_at_amin] 

251 if debug: 

252 t_at_amax = tr.tmin + i_at_amax * tr.deltat 

253 t_at_amin = tr.tmin + i_at_amin * tr.deltat 

254 ts = sorted([t_at_amax, t_at_amin]) 

255 mark = marker.Marker( 

256 [tr.nslc_id], 

257 ts[0], 

258 ts[1], 

259 0) 

260 

261 marker_candidates.append(mark) 

262 

263 component_amp_maxs.append(amax - amin) 

264 

265 i_at_amax = num.argmax(component_amp_maxs) 

266 if debug: 

267 markers.append(marker_candidates[i_at_amax]) 

268 amp_max = component_amp_maxs[i_at_amax] 

269 

270 elif self.method == 'peak_absolute_vector': 

271 trsum = None 

272 for tr in trs_processed: 

273 tr.set_ydata(tr.get_ydata()**2) 

274 if trsum is None: 

275 trsum = tr 

276 else: 

277 trsum.add(tr) 

278 

279 trsum.set_ydata(num.sqrt(tr.get_ydata)) 

280 trsum.set_codes(channel='SUM') 

281 

282 yabs = trsum.get_ydata() 

283 

284 i_at_amax = num.argmax(yabs) 

285 amax = yabs[i_at_amax] 

286 t_at_amax = tr.tmin + i_at_amax * tr.deltat 

287 amp_max = amax 

288 

289 if debug: 

290 markers.append(marker.Marker( 

291 [trsum.nslc_id], 

292 t_at_amax, 

293 t_at_amax, 

294 0)) 

295 

296 trs_processed.append(trsum) 

297 

298 elif self.method == 'spectral_average': 

299 component_amp_maxs = [] 

300 for tr in trs_processed: 

301 freqs, values = tr.spectrum() 

302 component_amp_maxs.append(num.mean(num.abs(values[ 

303 num.logical_and(self.fmin <= freqs, freqs <= self.fmax)]))) 

304 

305 amp_max = num.mean(component_amp_maxs) 

306 

307 if debug: 

308 trs_out = [] 

309 for tr in trs_orig: 

310 tr_out = tr.copy() 

311 tr_out.set_location(tr_out.location + '-' + self.name) 

312 trs_out.append(tr_out) 

313 

314 return amp_max, (trs_out + trs_processed, markers) 

315 

316 return amp_max, None