1import numpy as num
2from pyrocko import gf, trace
3from pyrocko.guts import Object, Float, StringChoice, List, String
4from pyrocko.gui import marker
5from grond.meta import store_t
8guts_prefix = 'grond'
11response_wwssn_lp = trace.PoleZeroResponse(
12 poles=[
13 -0.2513 + 0.3351J,
14 -0.2513 - 0.3351J,
15 -0.0628 + 0.0304J,
16 -0.0628 - 0.0304J],
17 zeros=[0., 0., 0.],
18 constant=383.6)
21response_wwssn_lp_2 = trace.PoleZeroResponse(
22 poles=[
23 -0.40180 + 0.08559J,
24 -0.40180 - 0.08559J,
25 -0.04841,
26 -0.08816],
27 zeros=[0., 0., 0.])
29response_wa = trace.PoleZeroResponse(
30 poles=[
31 -5.49779 - 5.60886J,
32 -5.49779 + 5.60886J],
33 zeros=[0., 0.])
36class NamedResponse(StringChoice):
37 choices = [
38 'wood-anderson',
39 'wwssn-lp']
41 map = {
42 'wood-anderson': response_wa,
43 'wwssn-lp': response_wwssn_lp}
46class BruneResponse(trace.FrequencyResponse):
48 duration = Float.T()
50 def evaluate(self, freqs):
51 return 1.0 / (1.0 + (freqs*self.duration)**2)
53# from pyrocko.plot import response as pr
54# pr.plot([response_wwssn_lp, response_wwssn_lp_2, response_wa])
57class FeatureMethod(StringChoice):
58 choices = [
59 'peak_component',
60 'peak_to_peak_component',
61 'peak_absolute_vector',
62 'spectral_average']
65class WaveformQuantity(StringChoice):
66 choices = [
67 'displacement',
68 'velocity',
69 'acceleration']
72class FeatureMeasurementFailed(Exception):
73 pass
76class FeatureMeasure(Object):
77 name = String.T()
78 timing_tmin = gf.Timing.T(default='vel:8')
79 timing_tmax = gf.Timing.T(default='vel:2')
80 fmin = Float.T(optional=True)
81 fmax = Float.T(optional=True)
82 response = trace.FrequencyResponse.T(optional=True)
83 named_response = NamedResponse.T(optional=True)
84 channels = List.T(String.T())
85 quantity = WaveformQuantity.T(default='displacement')
86 method = FeatureMethod.T(default='peak_component')
88 def get_nmodelling_targets(self):
89 return len(self.channels)
91 def get_modelling_targets(
92 self, codes, lat, lon, depth, store_id, backazimuth):
94 mtargets = []
95 for channel in self.channels:
96 target = gf.Target(
97 quantity='displacement',
98 codes=codes + (channel, ),
99 lat=lat,
100 lon=lon,
101 depth=depth,
102 store_id=store_id)
104 if channel == 'R':
105 target.azimuth = backazimuth - 180.
106 target.dip = 0.
107 elif channel == 'T':
108 target.azimuth = backazimuth - 90.
109 target.dip = 0.
110 elif channel == 'Z':
111 target.azimuth = 0.
112 target.dip = -90.
114 mtargets.append(target)
116 return mtargets
118 def evaluate(
119 self, engine, source, targets,
120 dataset=None,
121 trs=None,
122 extra_responses=[],
123 debug=False):
125 from ..waveform import target as base
127 trs_processed = []
128 trs_orig = []
129 for itarget, target in enumerate(targets):
130 if target.codes[-1] not in self.channels:
131 continue
133 store = engine.get_store(target.store_id)
135 tmin = source.time + store_t(
136 store, self.timing_tmin, source, target)
137 tmax = source.time + store_t(
138 store, self.timing_tmax, source, target)
140 if self.fmin is not None and self.fmax is not None:
141 freqlimits = [
142 self.fmin/2.,
143 self.fmin,
144 self.fmax,
145 self.fmax*2.]
146 tfade = 1./self.fmin
148 else:
149 freqlimits = None
150 tfade = 0.0
152 if dataset is not None:
153 bazi = base.backazimuth_for_waveform(
154 target.azimuth, target.codes)
156 tr = dataset.get_waveform(
157 target.codes,
158 tinc_cache=1.0/self.fmin,
159 quantity=self.quantity,
160 tmin=tmin,
161 tmax=tmax,
162 freqlimits=freqlimits,
163 tfade=tfade,
164 deltat=store.config.deltat,
165 cache=True,
166 backazimuth=bazi)
168 else:
169 tr = trs[itarget]
171 tr.extend(
172 tmin - tfade,
173 tmax + tfade,
174 fillmethod='repeat')
176 tr = tr.transfer(
177 freqlimits=freqlimits,
178 tfade=tfade)
180 trs_orig.append(tr)
182 tr = tr.copy()
184 responses = []
185 responses.extend(extra_responses)
187 ndiff = \
188 WaveformQuantity.choices.index(self.quantity) - \
189 WaveformQuantity.choices.index(target.quantity)
191 if ndiff > 0:
192 responses.append(trace.DifferentiationResponse(ndiff))
194 if ndiff < 0:
195 responses.append(trace.IntegrationResponse(-ndiff))
197 if self.response:
198 responses.append(self.response)
200 if self.named_response:
201 responses.append(
202 NamedResponse.map[self.named_response])
204 if responses:
205 trans = trace.MultiplyResponse(responses)
206 try:
207 tr = tr.transfer(transfer_function=trans)
209 except trace.TraceTooShort:
210 raise FeatureMeasurementFailed(
211 'transfer: trace too short')
213 if tmin is None or tmax is None:
214 raise FeatureMeasurementFailed(
215 'timing determination failed (phase unavailable?)')
217 tr.chop(tmin, tmax)
219 tr.set_location(tr.location + '-' + self.name + '-proc')
220 trs_processed.append(tr)
222 markers = []
223 marker_candidates = []
224 if self.method in ['peak_component', 'peak_to_peak_component']:
225 component_amp_maxs = []
226 for tr in trs_processed:
227 y = tr.get_ydata()
228 if self.method == 'peak_component':
229 yabs = num.abs(y)
230 i_at_amax = num.argmax(yabs)
231 amax = yabs[i_at_amax]
232 if debug:
233 t_at_amax = tr.tmin + i_at_amax * tr.deltat
234 mark = marker.Marker(
235 [tr.nslc_id],
236 t_at_amax,
237 t_at_amax,
238 0)
240 marker_candidates.append(mark)
242 component_amp_maxs.append(amax)
243 else:
244 i_at_amax = num.argmax(y)
245 i_at_amin = num.argmin(y)
246 amax = y[i_at_amax]
247 amin = y[i_at_amin]
248 if debug:
249 t_at_amax = tr.tmin + i_at_amax * tr.deltat
250 t_at_amin = tr.tmin + i_at_amin * tr.deltat
251 ts = sorted([t_at_amax, t_at_amin])
252 mark = marker.Marker(
253 [tr.nslc_id],
254 ts[0],
255 ts[1],
256 0)
258 marker_candidates.append(mark)
260 component_amp_maxs.append(amax - amin)
262 i_at_amax = num.argmax(component_amp_maxs)
263 if debug:
264 markers.append(marker_candidates[i_at_amax])
265 amp_max = component_amp_maxs[i_at_amax]
267 elif self.method == 'peak_absolute_vector':
268 trsum = None
269 for tr in trs_processed:
270 tr.set_ydata(tr.get_ydata()**2)
271 if trsum is None:
272 trsum = tr
273 else:
274 trsum.add(tr)
276 trsum.set_ydata(num.sqrt(tr.get_ydata))
277 trsum.set_codes(channel='SUM')
279 yabs = trsum.get_ydata()
281 i_at_amax = num.argmax(yabs)
282 amax = yabs[i_at_amax]
283 t_at_amax = tr.tmin + i_at_amax * tr.deltat
284 amp_max = amax
286 if debug:
287 markers.append(marker.Marker(
288 [trsum.nslc_id],
289 t_at_amax,
290 t_at_amax,
291 0))
293 trs_processed.append(trsum)
295 elif self.method == 'spectral_average':
296 component_amp_maxs = []
297 for tr in trs_processed:
298 freqs, values = tr.spectrum()
299 component_amp_maxs.append(num.mean(num.abs(values[
300 num.logical_and(self.fmin <= freqs, freqs <= self.fmax)])))
302 amp_max = num.mean(component_amp_maxs)
304 if debug:
305 trs_out = []
306 for tr in trs_orig:
307 tr_out = tr.copy()
308 tr_out.set_location(tr_out.location + '-' + self.name)
309 trs_out.append(tr_out)
311 return amp_max, (trs_out + trs_processed, markers)
313 return amp_max, None