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 for iwin in range(nwin):
416 tmin_win = tmin + iwin*tinc
417 tmax_win = tmin + (iwin+1)*tinc
419 if have_waveforms(tmin_win, tmax_win):
420 continue
422 if pbar is None:
423 pbar = util.progressbar('Generating waveforms', (nwin-iwin))
425 pbar.update(iwin)
427 trs = self.get_waveforms(engine, sources, tmin_win, tmax_win)
429 try:
430 wpaths = io.save(
431 trs, path_traces,
432 additional=dict(
433 wmin_year=tts(tmin_win, format='%Y'),
434 wmin_month=tts(tmin_win, format='%m'),
435 wmin_day=tts(tmin_win, format='%d'),
436 wmin=tts(tmin_win, format='%Y-%m-%d_%H-%M-%S'),
437 wmax_year=tts(tmax_win, format='%Y'),
438 wmax_month=tts(tmax_win, format='%m'),
439 wmax_day=tts(tmax_win, format='%d'),
440 wmax=tts(tmax_win, format='%Y-%m-%d_%H-%M-%S')))
442 for wpath in wpaths:
443 logger.debug('Generated file: %s' % wpath)
445 add_files(wpaths)
447 except FileSaveError as e:
448 raise ScenarioError(str(e))
450 if pbar is not None:
451 pbar.finish()
453 def ensure_responses(self, path):
454 from pyrocko.io import stationxml
456 path_responses = op.join(path, 'meta')
457 util.ensuredir(path_responses)
459 fn_stationxml = op.join(path_responses, 'waveform_response.xml')
460 i = 0
461 while op.exists(fn_stationxml):
462 fn_stationxml = op.join(
463 path_responses, 'waveform_response-%s.xml' % i)
464 i += 1
466 logger.debug('Writing waveform meta information to StationXML...')
468 stations = self.get_stations()
469 sxml = stationxml.FDSNStationXML.from_pyrocko_stations(stations)
471 sunit = {
472 'displacement': 'M',
473 'velocity': 'M/S',
474 'acceleration': 'M/S**2',
475 'counts': 'COUNTS'}[self.seismogram_quantity]
477 response = stationxml.Response(
478 instrument_sensitivity=stationxml.Sensitivity(
479 value=1.,
480 frequency=1.,
481 input_units=stationxml.Units(sunit),
482 output_units=stationxml.Units('COUNTS')),
483 stage_list=[])
485 for net, station, channel in sxml.iter_network_station_channels():
486 channel.response = response
488 sxml.dump_xml(filename=fn_stationxml)
490 def add_map_artists(self, engine, sources, automap):
491 automap.add_stations(self.get_stations())