Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/targets/waveform.py: 92%
239 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-01-03 09:20 +0000
1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Synthetic seismic waveform generator.
8'''
10import hashlib
11import math
12import logging
13import numpy as num
15from os import path as op
16from functools import reduce
18from pyrocko.guts import Float, List, Bool, Dict, String
19from pyrocko.gui.snuffler.marker import PhaseMarker, EventMarker
20from pyrocko import gf, model, util, trace, io
21from pyrocko.io.io_common import FileSaveError
22from pyrocko import pile
24from ..station import StationGenerator, RandomStationGenerator
25from .base import TargetGenerator, NoiseGenerator
26from ..error import ScenarioError
29DEFAULT_STORE_ID = 'global_2s'
31logger = logging.getLogger('pyrocko.scenario.targets.waveform')
32guts_prefix = 'pf.scenario'
35class WaveformNoiseGenerator(NoiseGenerator):
37 def get_time_increment(self, deltat):
38 return deltat * 1024
40 def get_intersecting_snippets(self, deltat, codes, tmin, tmax):
41 raise NotImplementedError()
43 def add_noise(self, tr):
44 for ntr in self.get_intersecting_snippets(
45 tr.deltat, tr.nslc_id, tr.tmin, tr.tmax):
46 tr.add(ntr)
49class WhiteNoiseGenerator(WaveformNoiseGenerator):
51 scale = Float.T(default=1e-6)
53 def get_seed_offset2(self, deltat, iw, codes):
54 m = hashlib.sha1(('%e %i %s.%s.%s.%s' % ((deltat, iw) + codes))
55 .encode('utf8'))
56 return int(m.hexdigest(), base=16) % 10000000
58 def get_intersecting_snippets(self, deltat, codes, tmin, tmax):
59 tinc = self.get_time_increment(deltat)
60 iwmin = int(math.floor(tmin / tinc))
61 iwmax = int(math.floor(tmax / tinc))
63 trs = []
64 for iw in range(iwmin, iwmax+1):
65 seed_offset = self.get_seed_offset2(deltat, iw, codes)
66 rstate = self.get_rstate(seed_offset)
68 n = int(round(tinc // deltat))
70 trs.append(trace.Trace(
71 codes[0], codes[1], codes[2], codes[3],
72 deltat=deltat,
73 tmin=iw*tinc,
74 ydata=rstate.normal(loc=0, scale=self.scale, size=n)))
76 return trs
79class WaveformGenerator(TargetGenerator):
81 station_generators = List.T(
82 StationGenerator.T(),
83 default=[RandomStationGenerator.D()],
84 help='List of StationGenerators.')
86 noise_generator = WaveformNoiseGenerator.T(
87 default=WhiteNoiseGenerator.D(),
88 help='Add Synthetic noise on the waveforms.')
90 store_id = gf.StringID.T(
91 default=DEFAULT_STORE_ID,
92 help='The GF store to use for forward-calculations.')
94 seismogram_quantity = gf.QuantityType.T(
95 default='displacement')
97 nsl_to_store_id = Dict.T(
98 String.T(), gf.StringID.T(),
99 help='Selectively use different GF stores for different stations.')
101 vmin_cut = Float.T(
102 default=2000.,
103 help='Minimum velocity to seismic velocity to consider in the model.')
104 vmax_cut = Float.T(
105 default=8000.,
106 help='Maximum velocity to seismic velocity to consider in the model.')
108 fmin = Float.T(
109 default=0.01,
110 help='Minimum frequency/wavelength to resolve in the'
111 ' synthetic waveforms.')
113 tabulated_phases = List.T(
114 gf.meta.TPDef.T(), optional=True,
115 help='Define seismic phases to be calculated.')
117 tabulated_phases_from_store = Bool.T(
118 default=False,
119 help='Calculate seismic phase arrivals for all travel-time tables '
120 'defined in GF store.')
122 tabulated_phases_noise_scale = Float.T(
123 default=0.0,
124 help='Standard deviation of normally distributed noise added to '
125 'calculated phase arrivals.')
127 taper = trace.Taper.T(
128 optional=True,
129 help='Time domain taper applied to synthetic waveforms.')
131 compensate_synthetic_offsets = Bool.T(
132 default=False,
133 help='Center synthetic trace amplitudes using mean of waveform tips.')
135 tinc = Float.T(
136 optional=True,
137 help='Time increment of waveforms.')
139 continuous = Bool.T(
140 default=True,
141 help='Only produce traces that intersect with events.')
143 def __init__(self, *args, **kwargs):
144 super(WaveformGenerator, self).__init__(*args, **kwargs)
145 self._targets = []
146 self._piles = {}
148 def _get_pile(self, path):
149 apath = op.abspath(path)
150 assert op.isdir(apath)
152 if apath not in self._piles:
153 fns = util.select_files(
154 [apath], show_progress=False)
156 p = pile.Pile()
157 if fns:
158 p.load_files(fns, fileformat='mseed', show_progress=False)
160 self._piles[apath] = p
162 return self._piles[apath]
164 def get_stations(self):
165 stations = []
166 for station_generator in self.station_generators:
167 stations.extend(station_generator.get_stations())
168 return stations
170 def get_distance_range(self, sources):
171 distances = num.array(
172 [sg.get_distance_range(sources)
173 for sg in self.station_generators])
174 return (distances[:, 0].min(), distances[:, 1].max())
176 def get_targets(self):
177 if self._targets:
178 return self._targets
180 for station in self.get_stations():
181 channel_data = []
182 channels = station.get_channels()
183 if channels:
184 for channel in channels:
185 channel_data.append([
186 channel.name,
187 channel.azimuth,
188 channel.dip])
190 else:
191 for c_name in ['BHZ', 'BHE', 'BHN']:
192 channel_data.append([
193 c_name,
194 model.guess_azimuth_from_name(c_name),
195 model.guess_dip_from_name(c_name)])
197 nsl = (station.network, station.station, station.location)
198 for c_name, c_azi, c_dip in channel_data:
199 target = gf.Target(
200 codes=nsl + (c_name,),
201 quantity=self.seismogram_quantity,
202 lat=station.lat,
203 lon=station.lon,
204 north_shift=station.north_shift,
205 east_shift=station.east_shift,
206 depth=station.depth,
207 store_id=self.nsl_to_store_id.get(
208 '.'.join(nsl), self.store_id),
209 optimization='enable',
210 interpolation='nearest_neighbor',
211 azimuth=c_azi,
212 dip=c_dip)
214 self._targets.append(target)
216 return self._targets
218 def get_time_range(self, sources):
219 dmin, dmax = self.get_distance_range(sources)
221 times = num.array([source.time for source in sources],
222 dtype=util.get_time_dtype())
224 tmin_events = num.min(times)
225 tmax_events = num.max(times)
227 tmin = tmin_events + dmin / self.vmax_cut - 10.0 / self.fmin
228 tmax = tmax_events + dmax / self.vmin_cut + 10.0 / self.fmin
230 return tmin, tmax
232 def get_codes_to_deltat(self, engine, sources):
233 deltats = {}
235 targets = self.get_targets()
236 for source in sources:
237 for target in targets:
238 deltats[target.codes] = engine.get_store(
239 target.store_id).config.deltat
241 return deltats
243 def get_useful_time_increment(self, engine, sources):
244 _, dmax = self.get_distance_range(sources)
245 tinc = dmax / self.vmin_cut + 2.0 / self.fmin
247 deltats = set(self.get_codes_to_deltat(engine, sources).values())
248 deltat = reduce(util.lcm, deltats)
249 tinc = int(round(tinc / deltat)) * deltat
250 return tinc
252 def get_relevant_sources(self, sources, tmin, tmax):
253 dmin, dmax = self.get_distance_range(sources)
254 trange = tmax - tmin
255 tmax_pad = trange + tmax + dmin / self.vmax_cut
256 tmin_pad = tmin - (dmax / self.vmin_cut + trange)
258 return [s for s in sources if s.time < tmax_pad and s.time > tmin_pad]
260 def get_waveforms(self, engine, sources, tmin, tmax):
262 sources_relevant = self.get_relevant_sources(sources, tmin, tmax)
263 if not (self.continuous or sources_relevant):
264 return []
266 trs = {}
267 tts = util.time_to_str
269 for nslc, deltat in self.get_codes_to_deltat(engine, sources).items():
270 tr_tmin = round(tmin / deltat) * deltat
271 tr_tmax = (round(tmax / deltat)-1) * deltat
272 nsamples = int(round((tr_tmax - tr_tmin) / deltat)) + 1
274 tr = trace.Trace(
275 *nslc,
276 tmin=tr_tmin,
277 ydata=num.zeros(nsamples),
278 deltat=deltat)
280 self.noise_generator.add_noise(tr)
282 trs[nslc] = tr
284 logger.debug('Forward modelling waveforms between %s - %s...'
285 % (tts(tmin, format='%Y-%m-%d_%H-%M-%S'),
286 tts(tmax, format='%Y-%m-%d_%H-%M-%S')))
288 if not sources_relevant:
289 return list(trs.values())
291 targets = self.get_targets()
292 response = engine.process(sources_relevant, targets)
293 for source, target, res in response.iter_results(
294 get='results'):
296 if isinstance(res, gf.SeismosizerError):
297 logger.warning(
298 'Out of bounds! \nTarget: %s\nSource: %s\n' % (
299 '.'.join(target.codes)), source)
300 continue
302 tr = res.trace.pyrocko_trace()
304 candidate = trs[target.codes]
305 if not candidate.overlaps(tr.tmin, tr.tmax):
306 continue
308 if self.compensate_synthetic_offsets:
309 tr.ydata -= (num.mean(tr.ydata[-3:-1]) +
310 num.mean(tr.ydata[1:3])) / 2.
312 if self.taper:
313 tr.taper(self.taper)
315 candidate.add(tr)
316 trs[target.codes] = candidate
318 return list(trs.values())
320 def get_onsets(self, engine, sources, *args, **kwargs):
322 targets = {t.codes[:3]: t for t in self.get_targets()}
324 markers = []
325 for source in sources:
326 ev = source.pyrocko_event()
327 markers.append(EventMarker(ev))
328 for nsl, target in targets.items():
329 store = engine.get_store(target.store_id)
330 if self.tabulated_phases:
331 tabulated_phases = self.tabulated_phases
333 elif self.tabulated_phases_from_store:
334 tabulated_phases = store.config.tabulated_phases
335 else:
336 tabulated_phases = []
338 for phase in tabulated_phases:
339 t = store.t(phase.id, source, target)
340 if not t:
341 continue
343 noise_scale = self.tabulated_phases_noise_scale
344 if noise_scale != 0.0:
345 t += num.random.normal(scale=noise_scale)
347 t += source.time
348 markers.append(
349 PhaseMarker(
350 phasename=phase.id,
351 tmin=t,
352 tmax=t,
353 event=source.pyrocko_event(),
354 nslc_ids=(nsl+('*',),)
355 )
356 )
357 return markers
359 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
360 self.ensure_waveforms(engine, sources, path, tmin, tmax)
361 self.ensure_responses(path)
363 def ensure_waveforms(self, engine, sources, path, tmin=None, tmax=None):
365 path_waveforms = op.join(path, 'waveforms')
366 util.ensuredir(path_waveforms)
368 p = self._get_pile(path_waveforms)
370 nslc_ids = set(target.codes for target in self.get_targets())
372 def have_waveforms(tmin, tmax):
373 trs_have = p.all(
374 tmin=tmin, tmax=tmax,
375 load_data=False, degap=False,
376 trace_selector=lambda tr: tr.nslc_id in nslc_ids)
378 return any(tr.data_len() > 0 for tr in trs_have)
380 def add_files(paths):
381 p.load_files(paths, fileformat='mseed', show_progress=False)
383 path_traces = op.join(
384 path_waveforms,
385 '%(wmin_year)s',
386 '%(wmin_month)s',
387 '%(wmin_day)s',
388 'waveform_%(network)s_%(station)s_' +
389 '%(location)s_%(channel)s_%(tmin)s_%(tmax)s.mseed')
391 tmin_all, tmax_all = self.get_time_range(sources)
392 tmin = tmin if tmin is not None else tmin_all
393 tmax = tmax if tmax is not None else tmax_all
394 tts = util.time_to_str
396 tinc = self.tinc or self.get_useful_time_increment(engine, sources)
397 tmin = num.floor(tmin / tinc) * tinc
398 tmax = num.ceil(tmax / tinc) * tinc
400 nwin = int(round((tmax - tmin) / tinc))
402 pbar = None
403 try:
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(
413 'Generating waveforms', (nwin-iwin))
415 pbar.update(iwin)
417 trs = self.get_waveforms(engine, sources, tmin_win, tmax_win)
419 try:
420 wpaths = io.save(
421 trs, path_traces,
422 additional=dict(
423 wmin_year=tts(tmin_win, format='%Y'),
424 wmin_month=tts(tmin_win, format='%m'),
425 wmin_day=tts(tmin_win, format='%d'),
426 wmin=tts(tmin_win, format='%Y-%m-%d_%H-%M-%S'),
427 wmax_year=tts(tmax_win, format='%Y'),
428 wmax_month=tts(tmax_win, format='%m'),
429 wmax_day=tts(tmax_win, format='%d'),
430 wmax=tts(tmax_win, format='%Y-%m-%d_%H-%M-%S')))
432 for wpath in wpaths:
433 logger.debug('Generated file: %s' % wpath)
435 add_files(wpaths)
437 except FileSaveError as e:
438 raise ScenarioError(str(e))
440 finally:
441 if pbar is not None:
442 pbar.finish()
444 def ensure_responses(self, path):
445 from pyrocko.io import stationxml
447 path_responses = op.join(path, 'meta')
448 util.ensuredir(path_responses)
450 fn_stationxml = op.join(path_responses, 'waveform_response.xml')
451 i = 0
452 while op.exists(fn_stationxml):
453 fn_stationxml = op.join(
454 path_responses, 'waveform_response-%s.xml' % i)
455 i += 1
457 logger.debug('Writing waveform meta information to StationXML...')
459 stations = self.get_stations()
460 sxml = stationxml.FDSNStationXML.from_pyrocko_stations(stations)
462 sunit = stationxml.quantity_to_units[self.seismogram_quantity]
464 response = stationxml.Response(
465 instrument_sensitivity=stationxml.Sensitivity(
466 value=1.,
467 frequency=1.,
468 input_units=stationxml.Units(sunit),
469 output_units=stationxml.Units('COUNT')),
470 stage_list=[])
472 for net, station, channel in sxml.iter_network_station_channels():
473 channel.response = response
475 sxml.dump_xml(filename=fn_stationxml)
477 def add_map_artists(self, engine, sources, automap):
478 automap.add_stations(self.get_stations())