Source code for seismic.synthetics.backends.backend_tws

#!/usr/bin/env python
"""
Backend for making synthetic seismograms using Telewavesim.
"""

import numpy as np
import scipy.signal as sig
from telewavesim.utils import Model, run_plane
import obspy

from seismic.synthetics.backends.synthesizer_base import Synthesizer
from seismic.model_properties import LayerProps
from seismic.units_utils import KM_PER_DEG
from seismic.stream_processing import zne_order

# pylint: disable=invalid-name


[docs]def synthesizer(): """Getter for backend Synthesizer class :return: Class name :rtype: SynthesizerMatrixPropagator """ return SynthesizerMatrixPropagator
# end func
[docs]class SynthesizerMatrixPropagator(Synthesizer): """ Class to synthesize seismograms from a 1D model description using Kennett's matrix propagator method. """ def __init__(self, station_latlon, layerprops): """ Initialization :param station_latlon: See documentation for :func:`~seismic.synthetics.backends.synthesizer_base.Synthesizer.synthesize` :param layerprops: List of LayerProps. Last layer should be mantle properties. :type layerprops: list(seismic.model_properties.LayerProps) """ super().__init__(station_latlon) Vp = [layer.Vp for layer in layerprops] Vs = [layer.Vs for layer in layerprops] self.model = Model([layer.H for layer in layerprops], [layer.rho for layer in layerprops], Vp, Vs, ['iso']*len(layerprops)) self._kappa = np.array(Vp)/np.array(Vs) # end func @property def kappa(self): """Return ratio of Vp/Vs for each layer :return: *k* value per layer :rtype: numpy.array """ return self._kappa # end func
[docs] def synthesize(self, src_latlon, fs, time_window): """ See documentation for :func:`~seismic.synthetics.backends.synthesizer_base.Synthesizer.synthesize` """ duration = time_window[1] - time_window[0] npts = int(np.ceil(duration*fs)) dt = 1.0/fs stream_all = obspy.Stream() for src_lat, src_lon in src_latlon: stream = self._synthone(src_lat, src_lon, dt, npts, time_window[0]) stream.traces = sorted(stream.traces, key=zne_order) for tr in stream: tr.stats.starttime = tr.stats.onset + time_window[0] # end for stream_all += stream # end for return stream_all
# end func def _synthone(self, src_lat, src_lon, dt, npts, time_init): """ Synthesize single event """ event_id_base = 'hiperseis:PRM/evid=' stats = self.compute_event_stats(src_lat, src_lon, event_id_base) stats.pop('tt_model') ray_param_sec_per_km = stats['slowness']/KM_PER_DEG traces_zne = run_plane(self.model, ray_param_sec_per_km, npts, dt, baz=stats['back_azimuth']) # Taper beginning to deal with spectral artefacts traces_zne.taper(0.05, side='left') # Synthetics here have been observed to have high noise at Nyquist freq. Filter it out. kernel = sig.firwin(9, 0.5) for tr in traces_zne: tr.data = sig.filtfilt(kernel, 1, tr.data) # end for # Extract time location of max Z amplitude and use as the onset time marker. tr_z = traces_zne.select(component='Z')[0] onset_marker = tr_z.times()[(tr_z.data == np.max(tr_z.data))][0] times_rel_mantle = tr_z.times()[0] - onset_marker for tr in traces_zne: stats['channel'] = '**' + tr.stats['channel'][-1] tr.stats.update(stats) if time_init < times_rel_mantle: n_leading = int((times_rel_mantle - time_init)/dt) data_new = np.append(np.zeros(n_leading), tr.data[:(npts - n_leading)]) assert data_new.shape == tr.data.shape tr.data = data_new elif time_init > times_rel_mantle: n_trailing = int((time_init - times_rel_mantle)/dt) data_new = np.append(tr.data[(npts - n_trailing):], np.zeros(n_trailing)) assert data_new.shape == tr.data.shape tr.data = data_new # end if # end for traces_zne.differentiate() return traces_zne
# end func # end class if __name__ == '__main__': example = SynthesizerMatrixPropagator((-20, 160), [LayerProps(6.4, 4.2, 2.7, 35.0), LayerProps(8.2, 6.8, 3.3, np.nan)]) test_strm = example.synthesize(((30, 150),), 100.0, (-10, 30)) print(test_strm) # end if