1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import hashlib
8import math
9import logging
10import numpy as num
12from os import path as op
13from functools import reduce
15from pyrocko.guts import StringChoice, Float, List, Bool
16from pyrocko.gui.marker import PhaseMarker, EventMarker
17from pyrocko import gf, model, util, trace, io
18from pyrocko.io.io_common import FileSaveError
19from pyrocko import pile
21from ..station import StationGenerator, RandomStationGenerator
22from .base import TargetGenerator, NoiseGenerator
23from ..error import ScenarioError
26DEFAULT_STORE_ID = 'global_2s'
28logger = logging.getLogger('pyrocko.scenario.targets.waveform')
29guts_prefix = 'pf.scenario'
32class WaveformNoiseGenerator(NoiseGenerator):
34 def get_time_increment(self, deltat):
35 return deltat * 1024
37 def get_intersecting_snippets(self, deltat, codes, tmin, tmax):
38 raise NotImplementedError()
40 def add_noise(self, tr):
41 for ntr in self.get_intersecting_snippets(
42 tr.deltat, tr.nslc_id, tr.tmin, tr.tmax):
43 tr.add(ntr)
46class WhiteNoiseGenerator(WaveformNoiseGenerator):
48 scale = Float.T(default=1e-6)
50 def get_seed_offset2(self, deltat, iw, codes):
51 m = hashlib.sha1(('%e %i %s.%s.%s.%s' % ((deltat, iw) + codes))
52 .encode('utf8'))
53 return int(m.hexdigest(), base=16) % 10000000
55 def get_intersecting_snippets(self, deltat, codes, tmin, tmax):
56 tinc = self.get_time_increment(deltat)
57 iwmin = int(math.floor(tmin / tinc))
58 iwmax = int(math.floor(tmax / tinc))
60 trs = []
61 for iw in range(iwmin, iwmax+1):
62 seed_offset = self.get_seed_offset2(deltat, iw, codes)
63 rstate = self.get_rstate(seed_offset)
65 n = int(round(tinc // deltat))
67 trs.append(trace.Trace(
68 codes[0], codes[1], codes[2], codes[3],
69 deltat=deltat,
70 tmin=iw*tinc,
71 ydata=rstate.normal(loc=0, scale=self.scale, size=n)))
73 return trs
76class WaveformGenerator(TargetGenerator):
78 station_generator = StationGenerator.T(
79 default=RandomStationGenerator.D(),
80 help='The StationGenerator for creating the stations.')
82 noise_generator = WaveformNoiseGenerator.T(
83 default=WhiteNoiseGenerator.D(),
84 help='Add Synthetic noise on the waveforms.')
86 store_id = gf.StringID.T(
87 default=DEFAULT_STORE_ID,
88 help='The GF store to use for forward-calculations.')
90 seismogram_quantity = StringChoice.T(
91 choices=['displacement', 'velocity', 'acceleration', 'counts'],
92 default='displacement')
94 vmin_cut = Float.T(
95 default=2000.,
96 help='Minimum velocity to seismic velicty to consider in the model.')
97 vmax_cut = Float.T(
98 default=8000.,
99 help='Maximum velocity to seismic velicty to consider in the model.')
101 fmin = Float.T(
102 default=0.01,
103 help='Minimum frequency/wavelength to resolve in the'
104 ' synthetic waveforms.')
106 tabulated_phases = List.T(
107 gf.meta.TPDef.T(), optional=True,
108 help='Define seismic phases to be calculated.')
110 tabulated_phases_from_store = Bool.T(
111 default=False,
112 help='Calculate seismic phase arrivals for all travel-time tables '
113 'defined in GF store.')
115 tabulated_phases_noise_scale = Float.T(
116 default=0.0,
117 help='Standard deviation of normally distributed noise added to '
118 'calculated phase arrivals.')
120 taper = trace.Taper.T(
121 optional=True,
122 help='Time domain taper applied to synthetic waveforms.')
124 compensate_synthetic_offsets = Bool.T(
125 default=False,
126 help='Center synthetic trace amplitudes using mean of waveform tips.')
128 tinc = Float.T(
129 optional=True,
130 help='Time increment of waveforms.')
132 continuous = Bool.T(
133 default=True,
134 help='Only produce traces that intersect with events.')
136 def __init__(self, *args, **kwargs):
137 super(WaveformGenerator, self).__init__(*args, **kwargs)
138 self._targets = []
139 self._piles = {}
141 def _get_pile(self, path):
142 apath = op.abspath(path)
143 assert op.isdir(apath)
145 if apath not in self._piles:
146 fns = util.select_files(
147 [apath], show_progress=False)
149 p = pile.Pile()
150 if fns:
151 p.load_files(fns, fileformat='mseed', show_progress=False)
153 self._piles[apath] = p
155 return self._piles[apath]
157 def get_stations(self):
158 return self.station_generator.get_stations()
160 def get_targets(self):
161 if self._targets:
162 return self._targets
164 for station in self.get_stations():
165 channel_data = []
166 channels = station.get_channels()
167 if channels:
168 for channel in channels:
169 channel_data.append([
170 channel.name,
171 channel.azimuth,
172 channel.dip])
174 else:
175 for c_name in ['BHZ', 'BHE', 'BHN']:
176 channel_data.append([
177 c_name,
178 model.guess_azimuth_from_name(c_name),
179 model.guess_dip_from_name(c_name)])
181 for c_name, c_azi, c_dip in channel_data:
183 target = gf.Target(
184 codes=(
185 station.network,
186 station.station,
187 station.location,
188 c_name),
189 quantity='displacement',
190 lat=station.lat,
191 lon=station.lon,
192 north_shift=station.north_shift,
193 east_shift=station.east_shift,
194 depth=station.depth,
195 store_id=self.store_id,
196 optimization='enable',
197 interpolation='nearest_neighbor',
198 azimuth=c_azi,
199 dip=c_dip)
201 self._targets.append(target)
203 return self._targets
205 def get_time_range(self, sources):
206 dmin, dmax = self.station_generator.get_distance_range(sources)
208 times = num.array([source.time for source in sources],
209 dtype=util.get_time_dtype())
211 tmin_events = num.min(times)
212 tmax_events = num.max(times)
214 tmin = tmin_events + dmin / self.vmax_cut - 10.0 / self.fmin
215 tmax = tmax_events + dmax / self.vmin_cut + 10.0 / self.fmin
217 return tmin, tmax
219 def get_codes_to_deltat(self, engine, sources):
220 deltats = {}
222 targets = self.get_targets()
223 for source in sources:
224 for target in targets:
225 deltats[target.codes] = engine.get_store(
226 target.store_id).config.deltat
228 return deltats
230 def get_useful_time_increment(self, engine, sources):
231 _, dmax = self.station_generator.get_distance_range(sources)
232 tinc = dmax / self.vmin_cut + 2.0 / self.fmin
234 deltats = set(self.get_codes_to_deltat(engine, sources).values())
235 deltat = reduce(util.lcm, deltats)
236 tinc = int(round(tinc / deltat)) * deltat
237 return tinc
239 def get_relevant_sources(self, sources, tmin, tmax):
240 dmin, dmax = self.station_generator.get_distance_range(sources)
241 trange = tmax - tmin
242 tmax_pad = trange + tmax + dmin / self.vmax_cut
243 tmin_pad = tmin - (dmax / self.vmin_cut + trange)
245 return [s for s in sources if s.time < tmax_pad and s.time > tmin_pad]
247 def get_waveforms(self, engine, sources, tmin, tmax):
249 sources_relevant = self.get_relevant_sources(sources, tmin, tmax)
250 if not (self.continuous or sources_relevant):
251 return []
253 trs = {}
254 tts = util.time_to_str
256 for nslc, deltat in self.get_codes_to_deltat(engine, sources).items():
257 tr_tmin = round(tmin / deltat) * deltat
258 tr_tmax = (round(tmax / deltat)-1) * deltat
259 nsamples = int(round((tr_tmax - tr_tmin) / deltat)) + 1
261 tr = trace.Trace(
262 *nslc,
263 tmin=tr_tmin,
264 ydata=num.zeros(nsamples),
265 deltat=deltat)
267 self.noise_generator.add_noise(tr)
269 trs[nslc] = tr
271 logger.debug('Forward modelling waveforms between %s - %s...'
272 % (tts(tmin, format='%Y-%m-%d_%H-%M-%S'),
273 tts(tmax, format='%Y-%m-%d_%H-%M-%S')))
275 if not sources_relevant:
276 return list(trs.values())
278 targets = self.get_targets()
279 response = engine.process(sources_relevant, targets)
280 for source, target, res in response.iter_results(
281 get='results'):
283 if isinstance(res, gf.SeismosizerError):
284 logger.warning(
285 'Out of bounds! \nTarget: %s\nSource: %s\n' % (
286 '.'.join(target.codes)), source)
287 continue
289 tr = res.trace.pyrocko_trace()
291 candidate = trs[target.codes]
292 if not candidate.overlaps(tr.tmin, tr.tmax):
293 continue
295 if self.compensate_synthetic_offsets:
296 tr.ydata -= (num.mean(tr.ydata[-3:-1]) +
297 num.mean(tr.ydata[1:3])) / 2.
299 if self.taper:
300 tr.taper(self.taper)
302 resp = self.get_transfer_function(target.codes)
303 if resp:
304 tr = tr.transfer(transfer_function=resp)
306 candidate.add(tr)
307 trs[target.codes] = candidate
309 return list(trs.values())
311 def get_onsets(self, engine, sources, *args, **kwargs):
313 targets = {t.codes[:3]: t for t in self.get_targets()}
315 markers = []
316 for source in sources:
317 ev = source.pyrocko_event()
318 markers.append(EventMarker(ev))
319 for nsl, target in targets.items():
320 store = engine.get_store(target.store_id)
321 if self.tabulated_phases:
322 tabulated_phases = self.tabulated_phases
324 elif self.tabulated_phases_from_store:
325 tabulated_phases = store.config.tabulated_phases
326 else:
327 tabulated_phases = []
329 for phase in tabulated_phases:
330 t = store.t(phase.id, source, target)
331 if not t:
332 continue
334 noise_scale = self.tabulated_phases_noise_scale
335 if noise_scale != 0.0:
336 t += num.random.normal(scale=noise_scale)
338 t += source.time
339 markers.append(
340 PhaseMarker(
341 phasename=phase.id,
342 tmin=t,
343 tmax=t,
344 event=source.pyrocko_event(),
345 nslc_ids=(nsl+('*',),)
346 )
347 )
348 return markers
350 def get_transfer_function(self, codes):
351 if self.seismogram_quantity == 'displacement':
352 return None
353 elif self.seismogram_quantity == 'velocity':
354 return trace.DifferentiationResponse(1)
355 elif self.seismogram_quantity == 'acceleration':
356 return trace.DifferentiationResponse(2)
357 elif self.seismogram_quantity == 'counts':
358 raise NotImplementedError()
360 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
361 self.ensure_waveforms(engine, sources, path, tmin, tmax)
362 self.ensure_responses(path)
364 def ensure_waveforms(self, engine, sources, path, tmin=None, tmax=None):
366 path_waveforms = op.join(path, 'waveforms')
367 util.ensuredir(path_waveforms)
369 p = self._get_pile(path_waveforms)
371 nslc_ids = set(target.codes for target in self.get_targets())
373 def have_waveforms(tmin, tmax):
374 trs_have = p.all(
375 tmin=tmin, tmax=tmax,
376 load_data=False, degap=False,
377 trace_selector=lambda tr: tr.nslc_id in nslc_ids)
379 return any(tr.data_len() > 0 for tr in trs_have)
381 def add_files(paths):
382 p.load_files(paths, fileformat='mseed', show_progress=False)
384 path_traces = op.join(
385 path_waveforms,
386 '%(wmin_year)s',
387 '%(wmin_month)s',
388 '%(wmin_day)s',
389 'waveform_%(network)s_%(station)s_' +
390 '%(location)s_%(channel)s_%(tmin)s_%(tmax)s.mseed')
392 tmin_all, tmax_all = self.get_time_range(sources)
393 tmin = tmin if tmin is not None else tmin_all
394 tmax = tmax if tmax is not None else tmax_all
395 tts = util.time_to_str
397 tinc = self.tinc or self.get_useful_time_increment(engine, sources)
398 tmin = num.floor(tmin / tinc) * tinc
399 tmax = num.ceil(tmax / tinc) * tinc
401 nwin = int(round((tmax - tmin) / tinc))
403 pbar = None
404 for iwin in range(nwin):
405 tmin_win = tmin + iwin*tinc
406 tmax_win = tmin + (iwin+1)*tinc
408 if have_waveforms(tmin_win, tmax_win):
409 continue
411 if pbar is None:
412 pbar = util.progressbar('Generating waveforms', (nwin-iwin))
414 pbar.update(iwin)
416 trs = self.get_waveforms(engine, sources, tmin_win, tmax_win)
418 try:
419 wpaths = io.save(
420 trs, path_traces,
421 additional=dict(
422 wmin_year=tts(tmin_win, format='%Y'),
423 wmin_month=tts(tmin_win, format='%m'),
424 wmin_day=tts(tmin_win, format='%d'),
425 wmin=tts(tmin_win, format='%Y-%m-%d_%H-%M-%S'),
426 wmax_year=tts(tmax_win, format='%Y'),
427 wmax_month=tts(tmax_win, format='%m'),
428 wmax_day=tts(tmax_win, format='%d'),
429 wmax=tts(tmax_win, format='%Y-%m-%d_%H-%M-%S')))
431 for wpath in wpaths:
432 logger.debug('Generated file: %s' % wpath)
434 add_files(wpaths)
436 except FileSaveError as e:
437 raise ScenarioError(str(e))
439 if pbar is not None:
440 pbar.finish()
442 def ensure_responses(self, path):
443 from pyrocko.io import stationxml
445 path_responses = op.join(path, 'meta')
446 util.ensuredir(path_responses)
448 fn_stationxml = op.join(path_responses, 'stations.xml')
449 if op.exists(fn_stationxml):
450 return
452 logger.debug('Writing waveform meta information to StationXML...')
454 stations = self.station_generator.get_stations()
455 sxml = stationxml.FDSNStationXML.from_pyrocko_stations(stations)
457 sunit = {
458 'displacement': 'M',
459 'velocity': 'M/S',
460 'acceleration': 'M/S**2',
461 'counts': 'COUNTS'}[self.seismogram_quantity]
463 response = stationxml.Response(
464 instrument_sensitivity=stationxml.Sensitivity(
465 value=1.,
466 frequency=1.,
467 input_units=stationxml.Units(sunit),
468 output_units=stationxml.Units('COUNTS')),
469 stage_list=[])
471 for net, station, channel in sxml.iter_network_station_channels():
472 channel.response = response
474 sxml.dump_xml(filename=fn_stationxml)
476 def add_map_artists(self, engine, sources, automap):
477 automap.add_stations(self.get_stations())