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.response import DifferentiationResponse
19from pyrocko.io.io_common import FileSaveError
20from pyrocko import pile
22from ..station import StationGenerator, RandomStationGenerator
23from .base import TargetGenerator, NoiseGenerator
24from ..error import ScenarioError
27DEFAULT_STORE_ID = 'global_2s'
29logger = logging.getLogger('pyrocko.scenario.targets.waveform')
30guts_prefix = 'pf.scenario'
33class WaveformNoiseGenerator(NoiseGenerator):
35 def get_time_increment(self, deltat):
36 return deltat * 1024
38 def get_intersecting_snippets(self, deltat, codes, tmin, tmax):
39 raise NotImplementedError()
41 def add_noise(self, tr):
42 for ntr in self.get_intersecting_snippets(
43 tr.deltat, tr.nslc_id, tr.tmin, tr.tmax):
44 tr.add(ntr)
47class WhiteNoiseGenerator(WaveformNoiseGenerator):
49 scale = Float.T(default=1e-6)
51 def get_seed_offset2(self, deltat, iw, codes):
52 m = hashlib.sha1(('%e %i %s.%s.%s.%s' % ((deltat, iw) + codes))
53 .encode('utf8'))
54 return int(m.hexdigest(), base=16) % 10000000
56 def get_intersecting_snippets(self, deltat, codes, tmin, tmax):
57 tinc = self.get_time_increment(deltat)
58 iwmin = int(math.floor(tmin / tinc))
59 iwmax = int(math.floor(tmax / tinc))
61 trs = []
62 for iw in range(iwmin, iwmax+1):
63 seed_offset = self.get_seed_offset2(deltat, iw, codes)
64 rstate = self.get_rstate(seed_offset)
66 n = int(round(tinc // deltat))
68 trs.append(trace.Trace(
69 codes[0], codes[1], codes[2], codes[3],
70 deltat=deltat,
71 tmin=iw*tinc,
72 ydata=rstate.normal(loc=0, scale=self.scale, size=n)))
74 return trs
77class WaveformGenerator(TargetGenerator):
79 station_generators = List.T(
80 StationGenerator.T(),
81 default=[RandomStationGenerator.D()],
82 help='List of StationGenerators.')
84 noise_generator = WaveformNoiseGenerator.T(
85 default=WhiteNoiseGenerator.D(),
86 help='Add Synthetic noise on the waveforms.')
88 store_id = gf.StringID.T(
89 default=DEFAULT_STORE_ID,
90 help='The GF store to use for forward-calculations.')
92 seismogram_quantity = StringChoice.T(
93 choices=['displacement', 'velocity', 'acceleration', 'counts'],
94 default='displacement')
96 vmin_cut = Float.T(
97 default=2000.,
98 help='Minimum velocity to seismic velocity to consider in the model.')
99 vmax_cut = Float.T(
100 default=8000.,
101 help='Maximum velocity to seismic velocity to consider in the model.')
103 fmin = Float.T(
104 default=0.01,
105 help='Minimum frequency/wavelength to resolve in the'
106 ' synthetic waveforms.')
108 tabulated_phases = List.T(
109 gf.meta.TPDef.T(), optional=True,
110 help='Define seismic phases to be calculated.')
112 tabulated_phases_from_store = Bool.T(
113 default=False,
114 help='Calculate seismic phase arrivals for all travel-time tables '
115 'defined in GF store.')
117 tabulated_phases_noise_scale = Float.T(
118 default=0.0,
119 help='Standard deviation of normally distributed noise added to '
120 'calculated phase arrivals.')
122 taper = trace.Taper.T(
123 optional=True,
124 help='Time domain taper applied to synthetic waveforms.')
126 compensate_synthetic_offsets = Bool.T(
127 default=False,
128 help='Center synthetic trace amplitudes using mean of waveform tips.')
130 tinc = Float.T(
131 optional=True,
132 help='Time increment of waveforms.')
134 continuous = Bool.T(
135 default=True,
136 help='Only produce traces that intersect with events.')
138 def __init__(self, *args, **kwargs):
139 super(WaveformGenerator, self).__init__(*args, **kwargs)
140 self._targets = []
141 self._piles = {}
143 def _get_pile(self, path):
144 apath = op.abspath(path)
145 assert op.isdir(apath)
147 if apath not in self._piles:
148 fns = util.select_files(
149 [apath], show_progress=False)
151 p = pile.Pile()
152 if fns:
153 p.load_files(fns, fileformat='mseed', show_progress=False)
155 self._piles[apath] = p
157 return self._piles[apath]
159 def get_stations(self):
160 stations = []
161 for station_generator in self.station_generators:
162 stations.extend(station_generator.get_stations())
163 return stations
165 def get_distance_range(self, sources):
166 distances = num.array(
167 [sg.get_distance_range(sources)
168 for sg in self.station_generators])
169 return (distances[:, 0].min(), distances[:, 1].max())
171 def get_targets(self):
172 if self._targets:
173 return self._targets
175 for station in self.get_stations():
176 channel_data = []
177 channels = station.get_channels()
178 if channels:
179 for channel in channels:
180 channel_data.append([
181 channel.name,
182 channel.azimuth,
183 channel.dip])
185 else:
186 for c_name in ['BHZ', 'BHE', 'BHN']:
187 channel_data.append([
188 c_name,
189 model.guess_azimuth_from_name(c_name),
190 model.guess_dip_from_name(c_name)])
192 for c_name, c_azi, c_dip in channel_data:
194 target = gf.Target(
195 codes=(
196 station.network,
197 station.station,
198 station.location,
199 c_name),
200 quantity='displacement',
201 lat=station.lat,
202 lon=station.lon,
203 north_shift=station.north_shift,
204 east_shift=station.east_shift,
205 depth=station.depth,
206 store_id=self.store_id,
207 optimization='enable',
208 interpolation='nearest_neighbor',
209 azimuth=c_azi,
210 dip=c_dip)
212 self._targets.append(target)
214 return self._targets
216 def get_time_range(self, sources):
217 dmin, dmax = self.get_distance_range(sources)
219 times = num.array([source.time for source in sources],
220 dtype=util.get_time_dtype())
222 tmin_events = num.min(times)
223 tmax_events = num.max(times)
225 tmin = tmin_events + dmin / self.vmax_cut - 10.0 / self.fmin
226 tmax = tmax_events + dmax / self.vmin_cut + 10.0 / self.fmin
228 return tmin, tmax
230 def get_codes_to_deltat(self, engine, sources):
231 deltats = {}
233 targets = self.get_targets()
234 for source in sources:
235 for target in targets:
236 deltats[target.codes] = engine.get_store(
237 target.store_id).config.deltat
239 return deltats
241 def get_useful_time_increment(self, engine, sources):
242 _, dmax = self.get_distance_range(sources)
243 tinc = dmax / self.vmin_cut + 2.0 / self.fmin
245 deltats = set(self.get_codes_to_deltat(engine, sources).values())
246 deltat = reduce(util.lcm, deltats)
247 tinc = int(round(tinc / deltat)) * deltat
248 return tinc
250 def get_relevant_sources(self, sources, tmin, tmax):
251 dmin, dmax = self.get_distance_range(sources)
252 trange = tmax - tmin
253 tmax_pad = trange + tmax + dmin / self.vmax_cut
254 tmin_pad = tmin - (dmax / self.vmin_cut + trange)
256 return [s for s in sources if s.time < tmax_pad and s.time > tmin_pad]
258 def get_waveforms(self, engine, sources, tmin, tmax):
260 sources_relevant = self.get_relevant_sources(sources, tmin, tmax)
261 if not (self.continuous or sources_relevant):
262 return []
264 trs = {}
265 tts = util.time_to_str
267 for nslc, deltat in self.get_codes_to_deltat(engine, sources).items():
268 tr_tmin = round(tmin / deltat) * deltat
269 tr_tmax = (round(tmax / deltat)-1) * deltat
270 nsamples = int(round((tr_tmax - tr_tmin) / deltat)) + 1
272 tr = trace.Trace(
273 *nslc,
274 tmin=tr_tmin,
275 ydata=num.zeros(nsamples),
276 deltat=deltat)
278 self.noise_generator.add_noise(tr)
280 trs[nslc] = tr
282 logger.debug('Forward modelling waveforms between %s - %s...'
283 % (tts(tmin, format='%Y-%m-%d_%H-%M-%S'),
284 tts(tmax, format='%Y-%m-%d_%H-%M-%S')))
286 if not sources_relevant:
287 return list(trs.values())
289 targets = self.get_targets()
290 response = engine.process(sources_relevant, targets)
291 for source, target, res in response.iter_results(
292 get='results'):
294 if isinstance(res, gf.SeismosizerError):
295 logger.warning(
296 'Out of bounds! \nTarget: %s\nSource: %s\n' % (
297 '.'.join(target.codes)), source)
298 continue
300 tr = res.trace.pyrocko_trace()
302 candidate = trs[target.codes]
303 if not candidate.overlaps(tr.tmin, tr.tmax):
304 continue
306 if self.compensate_synthetic_offsets:
307 tr.ydata -= (num.mean(tr.ydata[-3:-1]) +
308 num.mean(tr.ydata[1:3])) / 2.
310 if self.taper:
311 tr.taper(self.taper)
313 resp = self.get_transfer_function(target.codes)
314 if resp:
315 tr = tr.transfer(transfer_function=resp)
317 candidate.add(tr)
318 trs[target.codes] = candidate
320 return list(trs.values())
322 def get_onsets(self, engine, sources, *args, **kwargs):
324 targets = {t.codes[:3]: t for t in self.get_targets()}
326 markers = []
327 for source in sources:
328 ev = source.pyrocko_event()
329 markers.append(EventMarker(ev))
330 for nsl, target in targets.items():
331 store = engine.get_store(target.store_id)
332 if self.tabulated_phases:
333 tabulated_phases = self.tabulated_phases
335 elif self.tabulated_phases_from_store:
336 tabulated_phases = store.config.tabulated_phases
337 else:
338 tabulated_phases = []
340 for phase in tabulated_phases:
341 t = store.t(phase.id, source, target)
342 if not t:
343 continue
345 noise_scale = self.tabulated_phases_noise_scale
346 if noise_scale != 0.0:
347 t += num.random.normal(scale=noise_scale)
349 t += source.time
350 markers.append(
351 PhaseMarker(
352 phasename=phase.id,
353 tmin=t,
354 tmax=t,
355 event=source.pyrocko_event(),
356 nslc_ids=(nsl+('*',),)
357 )
358 )
359 return markers
361 def get_transfer_function(self, codes):
362 if self.seismogram_quantity == 'displacement':
363 return None
364 elif self.seismogram_quantity == 'velocity':
365 return DifferentiationResponse(1)
366 elif self.seismogram_quantity == 'acceleration':
367 return DifferentiationResponse(2)
368 elif self.seismogram_quantity == 'counts':
369 raise NotImplementedError()
371 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
372 self.ensure_waveforms(engine, sources, path, tmin, tmax)
373 self.ensure_responses(path)
375 def ensure_waveforms(self, engine, sources, path, tmin=None, tmax=None):
377 path_waveforms = op.join(path, 'waveforms')
378 util.ensuredir(path_waveforms)
380 p = self._get_pile(path_waveforms)
382 nslc_ids = set(target.codes for target in self.get_targets())
384 def have_waveforms(tmin, tmax):
385 trs_have = p.all(
386 tmin=tmin, tmax=tmax,
387 load_data=False, degap=False,
388 trace_selector=lambda tr: tr.nslc_id in nslc_ids)
390 return any(tr.data_len() > 0 for tr in trs_have)
392 def add_files(paths):
393 p.load_files(paths, fileformat='mseed', show_progress=False)
395 path_traces = op.join(
396 path_waveforms,
397 '%(wmin_year)s',
398 '%(wmin_month)s',
399 '%(wmin_day)s',
400 'waveform_%(network)s_%(station)s_' +
401 '%(location)s_%(channel)s_%(tmin)s_%(tmax)s.mseed')
403 tmin_all, tmax_all = self.get_time_range(sources)
404 tmin = tmin if tmin is not None else tmin_all
405 tmax = tmax if tmax is not None else tmax_all
406 tts = util.time_to_str
408 tinc = self.tinc or self.get_useful_time_increment(engine, sources)
409 tmin = num.floor(tmin / tinc) * tinc
410 tmax = num.ceil(tmax / tinc) * tinc
412 nwin = int(round((tmax - tmin) / tinc))
414 pbar = None
415 try:
416 for iwin in range(nwin):
417 tmin_win = tmin + iwin*tinc
418 tmax_win = tmin + (iwin+1)*tinc
420 if have_waveforms(tmin_win, tmax_win):
421 continue
423 if pbar is None:
424 pbar = util.progressbar(
425 'Generating waveforms', (nwin-iwin))
427 pbar.update(iwin)
429 trs = self.get_waveforms(engine, sources, tmin_win, tmax_win)
431 try:
432 wpaths = io.save(
433 trs, path_traces,
434 additional=dict(
435 wmin_year=tts(tmin_win, format='%Y'),
436 wmin_month=tts(tmin_win, format='%m'),
437 wmin_day=tts(tmin_win, format='%d'),
438 wmin=tts(tmin_win, format='%Y-%m-%d_%H-%M-%S'),
439 wmax_year=tts(tmax_win, format='%Y'),
440 wmax_month=tts(tmax_win, format='%m'),
441 wmax_day=tts(tmax_win, format='%d'),
442 wmax=tts(tmax_win, format='%Y-%m-%d_%H-%M-%S')))
444 for wpath in wpaths:
445 logger.debug('Generated file: %s' % wpath)
447 add_files(wpaths)
449 except FileSaveError as e:
450 raise ScenarioError(str(e))
452 finally:
453 if pbar is not None:
454 pbar.finish()
456 def ensure_responses(self, path):
457 from pyrocko.io import stationxml
459 path_responses = op.join(path, 'meta')
460 util.ensuredir(path_responses)
462 fn_stationxml = op.join(path_responses, 'waveform_response.xml')
463 i = 0
464 while op.exists(fn_stationxml):
465 fn_stationxml = op.join(
466 path_responses, 'waveform_response-%s.xml' % i)
467 i += 1
469 logger.debug('Writing waveform meta information to StationXML...')
471 stations = self.get_stations()
472 sxml = stationxml.FDSNStationXML.from_pyrocko_stations(stations)
474 sunit = {
475 'displacement': 'M',
476 'velocity': 'M/S',
477 'acceleration': 'M/S**2',
478 'counts': 'COUNTS'}[self.seismogram_quantity]
480 response = stationxml.Response(
481 instrument_sensitivity=stationxml.Sensitivity(
482 value=1.,
483 frequency=1.,
484 input_units=stationxml.Units(sunit),
485 output_units=stationxml.Units('COUNTS')),
486 stage_list=[])
488 for net, station, channel in sxml.iter_network_station_channels():
489 channel.response = response
491 sxml.dump_xml(filename=fn_stationxml)
493 def add_map_artists(self, engine, sources, automap):
494 automap.add_stations(self.get_stations())