Source code for pyrocko.scenario.targets.waveform

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
from __future__ import absolute_import, division, print_function

from builtins import range
import hashlib
import math
import logging
import numpy as num

from os import path as op
from functools import reduce

from pyrocko.guts import StringChoice, Float, List, Bool
from pyrocko.gui.marker import PhaseMarker, EventMarker
from pyrocko import gf, model, util, trace, io
from pyrocko.io_common import FileSaveError
from pyrocko import pile

from ..station import StationGenerator, RandomStationGenerator
from .base import TargetGenerator, NoiseGenerator
from ..error import ScenarioError


DEFAULT_STORE_ID = 'global_2s'

logger = logging.getLogger('pyrocko.scenario.targets.waveform')
guts_prefix = 'pf.scenario'


[docs]class WaveformNoiseGenerator(NoiseGenerator): def get_time_increment(self, deltat): return deltat * 1024 def get_intersecting_snippets(self, deltat, codes, tmin, tmax): raise NotImplementedError() def add_noise(self, tr): for ntr in self.get_intersecting_snippets( tr.deltat, tr.nslc_id, tr.tmin, tr.tmax): tr.add(ntr)
[docs]class WhiteNoiseGenerator(WaveformNoiseGenerator): scale = Float.T(default=1e-6) def get_seed_offset2(self, deltat, iw, codes): m = hashlib.sha1(('%e %i %s.%s.%s.%s' % ((deltat, iw) + codes)) .encode('utf8')) return int(m.hexdigest(), base=16) % 10000000 def get_intersecting_snippets(self, deltat, codes, tmin, tmax): tinc = self.get_time_increment(deltat) iwmin = int(math.floor(tmin / tinc)) iwmax = int(math.floor(tmax / tinc)) trs = [] for iw in range(iwmin, iwmax+1): seed_offset = self.get_seed_offset2(deltat, iw, codes) rstate = self.get_rstate(seed_offset) n = int(round(tinc // deltat)) trs.append(trace.Trace( codes[0], codes[1], codes[2], codes[3], deltat=deltat, tmin=iw*tinc, ydata=rstate.normal(loc=0, scale=self.scale, size=n))) return trs
[docs]class WaveformGenerator(TargetGenerator): station_generator = StationGenerator.T( default=RandomStationGenerator.D(), help='The StationGenerator for creating the stations.') noise_generator = WaveformNoiseGenerator.T( default=WhiteNoiseGenerator.D(), help='Add Synthetic noise on the waveforms.') store_id = gf.StringID.T( default=DEFAULT_STORE_ID, help='The GF store to use for forward-calculations.') seismogram_quantity = StringChoice.T( choices=['displacement', 'velocity', 'acceleration', 'counts'], default='displacement') vmin_cut = Float.T( default=2000., help='Minimum velocity to seismic velicty to consider in the model.') vmax_cut = Float.T( default=8000., help='Maximum velocity to seismic velicty to consider in the model.') fmin = Float.T( default=0.01, help='Minimum frequency/wavelength to resolve in the' ' synthetic waveforms.') tabulated_phases = List.T( gf.meta.TPDef.T(), optional=True, help='Define seismic phases to be calculated.') tabulated_phases_from_store = Bool.T( default=False, help='Calculate seismic phase arrivals for all travel-time tables ' 'defined in GF store.') tabulated_phases_noise_scale = Float.T( default=0.0, help='Standard deviation of normally distributed noise added to ' 'calculated phase arrivals.') taper = trace.Taper.T( optional=True, help='Time domain taper applied to synthetic waveforms.') compensate_synthetic_offsets = Bool.T( default=False, help='Center synthetic trace amplitudes using mean of waveform tips.') tinc = Float.T( optional=True, help='Time increment of waveforms.') continuous = Bool.T( default=True, help='Only produce traces that intersect with events.') def __init__(self, *args, **kwargs): super(WaveformGenerator, self).__init__(*args, **kwargs) self._targets = [] self._piles = {} def _get_pile(self, path): apath = op.abspath(path) assert op.isdir(apath) if apath not in self._piles: fns = util.select_files( [apath], show_progress=False) p = pile.Pile() if fns: p.load_files(fns, fileformat='mseed', show_progress=False) self._piles[apath] = p return self._piles[apath] def get_stations(self): return self.station_generator.get_stations()
[docs] def get_targets(self): if self._targets: return self._targets for station in self.get_stations(): channel_data = [] channels = station.get_channels() if channels: for channel in channels: channel_data.append([ channel.name, channel.azimuth, channel.dip]) else: for c_name in ['BHZ', 'BHE', 'BHN']: channel_data.append([ c_name, model.guess_azimuth_from_name(c_name), model.guess_dip_from_name(c_name)]) for c_name, c_azi, c_dip in channel_data: target = gf.Target( codes=( station.network, station.station, station.location, c_name), quantity='displacement', lat=station.lat, lon=station.lon, north_shift=station.north_shift, east_shift=station.east_shift, depth=station.depth, store_id=self.store_id, optimization='enable', interpolation='nearest_neighbor', azimuth=c_azi, dip=c_dip) self._targets.append(target) return self._targets
[docs] def get_time_range(self, sources): dmin, dmax = self.station_generator.get_distance_range(sources) times = num.array([source.time for source in sources], dtype=num.float) tmin_events = num.min(times) tmax_events = num.max(times) tmin = tmin_events + dmin / self.vmax_cut - 10.0 / self.fmin tmax = tmax_events + dmax / self.vmin_cut + 10.0 / self.fmin return tmin, tmax
def get_codes_to_deltat(self, engine, sources): deltats = {} targets = self.get_targets() for source in sources: for target in targets: deltats[target.codes] = engine.get_store( target.store_id).config.deltat return deltats def get_useful_time_increment(self, engine, sources): _, dmax = self.station_generator.get_distance_range(sources) tinc = dmax / self.vmin_cut + 2.0 / self.fmin deltats = set(self.get_codes_to_deltat(engine, sources).values()) deltat = reduce(util.lcm, deltats) tinc = int(round(tinc / deltat)) * deltat return tinc def get_relevant_sources(self, sources, tmin, tmax): dmin, dmax = self.station_generator.get_distance_range(sources) trange = tmax - tmin tmax_pad = trange + tmax + dmin / self.vmax_cut tmin_pad = tmin - (dmax / self.vmin_cut + trange) return [s for s in sources if s.time < tmax_pad and s.time > tmin_pad] def get_waveforms(self, engine, sources, tmin, tmax): sources_relevant = self.get_relevant_sources(sources, tmin, tmax) if not (self.continuous or sources_relevant): return [] trs = {} tts = util.time_to_str for nslc, deltat in self.get_codes_to_deltat(engine, sources).items(): tr_tmin = int(round(tmin / deltat)) * deltat tr_tmax = (int(round(tmax / deltat))-1) * deltat nsamples = int(round((tr_tmax - tr_tmin) / deltat)) + 1 tr = trace.Trace( *nslc, tmin=tr_tmin, ydata=num.zeros(nsamples), deltat=deltat) self.noise_generator.add_noise(tr) trs[nslc] = tr logger.debug('Forward modelling waveforms between %s - %s...' % (tts(tmin, format='%Y-%m-%d_%H-%M-%S'), tts(tmax, format='%Y-%m-%d_%H-%M-%S'))) if not sources_relevant: return list(trs.values()) targets = self.get_targets() response = engine.process(sources_relevant, targets) for source, target, res in response.iter_results( get='results'): if isinstance(res, gf.SeismosizerError): logger.warning( 'Out of bounds! \nTarget: %s\nSource: %s\n' % ( '.'.join(target.codes)), source) continue tr = res.trace.pyrocko_trace() candidate = trs[target.codes] if not candidate.overlaps(tr.tmin, tr.tmax): continue if self.compensate_synthetic_offsets: tr.ydata -= (num.mean(tr.ydata[-3:-1]) + num.mean(tr.ydata[1:3])) / 2. if self.taper: tr.taper(self.taper) resp = self.get_transfer_function(target.codes) if resp: tr = tr.transfer(transfer_function=resp) candidate.add(tr) trs[target.codes] = candidate return list(trs.values()) def get_onsets(self, engine, sources, *args, **kwargs): targets = {t.codes[:3]: t for t in self.get_targets()} markers = [] for source in sources: ev = source.pyrocko_event() markers.append(EventMarker(ev)) for nsl, target in targets.items(): store = engine.get_store(target.store_id) if self.tabulated_phases: tabulated_phases = self.tabulated_phases elif self.tabulated_phases_from_store: tabulated_phases = store.config.tabulated_phases else: tabulated_phases = [] for phase in tabulated_phases: t = store.t(phase.id, source, target) if not t: continue noise_scale = self.tabulated_phases_noise_scale if noise_scale != 0.0: t += num.random.normal(scale=noise_scale) t += source.time markers.append( PhaseMarker( phasename=phase.id, tmin=t, tmax=t, event=source.pyrocko_event(), nslc_ids=(nsl+('*',),) ) ) return markers def get_transfer_function(self, codes): if self.seismogram_quantity == 'displacement': return None elif self.seismogram_quantity == 'velocity': return trace.DifferentiationResponse(1) elif self.seismogram_quantity == 'acceleration': return trace.DifferentiationResponse(2) elif self.seismogram_quantity == 'counts': raise NotImplementedError() def ensure_data(self, engine, sources, path, tmin=None, tmax=None): self.ensure_waveforms(engine, sources, path, tmin, tmax) self.ensure_responses(path) def ensure_waveforms(self, engine, sources, path, tmin=None, tmax=None): path_waveforms = op.join(path, 'waveforms') util.ensuredir(path_waveforms) p = self._get_pile(path_waveforms) nslc_ids = set(target.codes for target in self.get_targets()) def have_waveforms(tmin, tmax): trs_have = p.all( tmin=tmin, tmax=tmax, load_data=False, degap=False, trace_selector=lambda tr: tr.nslc_id in nslc_ids) return any(tr.data_len() > 0 for tr in trs_have) def add_files(paths): p.load_files(paths, fileformat='mseed', show_progress=False) path_traces = op.join( path_waveforms, '%(wmin_year)s', '%(wmin_month)s', '%(wmin_day)s', 'waveform_%(network)s_%(station)s_' + '%(location)s_%(channel)s_%(tmin)s_%(tmax)s.mseed') tmin_all, tmax_all = self.get_time_range(sources) tmin = tmin if tmin is not None else tmin_all tmax = tmax if tmax is not None else tmax_all tts = util.time_to_str tinc = self.tinc or self.get_useful_time_increment(engine, sources) tmin = math.floor(tmin / tinc) * tinc tmax = math.ceil(tmax / tinc) * tinc nwin = int(round((tmax - tmin) / tinc)) pbar = None for iwin in range(nwin): tmin_win = tmin + iwin*tinc tmax_win = tmin + (iwin+1)*tinc if have_waveforms(tmin_win, tmax_win): continue if pbar is None: pbar = util.progressbar('Generating waveforms', (nwin-iwin)) pbar.update(iwin) trs = self.get_waveforms(engine, sources, tmin_win, tmax_win) try: wpaths = io.save( trs, path_traces, additional=dict( wmin_year=tts(tmin_win, format='%Y'), wmin_month=tts(tmin_win, format='%m'), wmin_day=tts(tmin_win, format='%d'), wmin=tts(tmin_win, format='%Y-%m-%d_%H-%M-%S'), wmax_year=tts(tmax_win, format='%Y'), wmax_month=tts(tmax_win, format='%m'), wmax_day=tts(tmax_win, format='%d'), wmax=tts(tmax_win, format='%Y-%m-%d_%H-%M-%S'))) for wpath in wpaths: logger.debug('Generated file: %s' % wpath) add_files(wpaths) except FileSaveError as e: raise ScenarioError(str(e)) if pbar is not None: pbar.finish() def ensure_responses(self, path): from pyrocko.io import stationxml path_responses = op.join(path, 'meta') util.ensuredir(path_responses) fn_stationxml = op.join(path_responses, 'stations.xml') if op.exists(fn_stationxml): return logger.debug('Writing waveform meta information to StationXML...') stations = self.station_generator.get_stations() sxml = stationxml.FDSNStationXML.from_pyrocko_stations(stations) sunit = { 'displacement': 'M', 'velocity': 'M/S', 'acceleration': 'M/S**2', 'counts': 'COUNTS'}[self.seismogram_quantity] response = stationxml.Response( instrument_sensitivity=stationxml.Sensitivity( value=1., frequency=1., input_units=stationxml.Units(sunit), output_units=stationxml.Units('COUNTS')), stage_list=[]) for net, station, channel in sxml.iter_network_station_channels(): channel.response = response sxml.dump_xml(filename=fn_stationxml) def add_map_artists(self, engine, sources, automap): automap.add_stations(self.get_stations())