Source code for seispy.harmonics

import numpy as np
from scipy.sparse.linalg import lsqr
from seispy.geo import cosd, sind
import matplotlib.pyplot as plt
from obspy.io.sac import SACTrace
from os.path import join


[docs] class Harmonics(): def __init__(self, rfsta, tmin=-5, tmax=10, bin_stack=True) -> None: """ Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs :param rfsta: data class of RFStation :type rfsta: :class:`seispy.rfcorrect.RFStation` :param tmin: Start time relative to P :type tmin: float :param tmax: End time relative to P :type tmax: float """ self.rfsta = rfsta self.tmin = tmin self.tmax = tmax if bin_stack: val = 10 stacked_data = self.rfsta.bin_stack(lim=[0, 360], val=val) self.datar = stacked_data['data_prime'] self.datat = stacked_data['datat'] self.bazi = np.arange(0, 360, val)+val/2 else: self.datar = self.rfsta.data_prime self.datat = self.rfsta.datat self.bazi = self.rfsta.bazi self.trace_num = self.datar.shape[0] self.cut_trace()
[docs] def cut_trace(self): """Trim RFs from tmin to tmax """ nb = int((self.tmin + self.rfsta.shift) / self.rfsta.sampling) ne = int((self.tmax + self.rfsta.shift) / self.rfsta.sampling) self.datar = self.datar[:, nb:ne+1] self.datat = self.datat[:, nb:ne+1] self.nsamp = ne - nb + 1 self.time_axis = np.linspace(self.tmin, self.tmax, self.nsamp)
[docs] def harmo_trans(self): """Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs """ self.traces = np.vstack((self.datar, self.datat)) harmonic = np.zeros((self.trace_num*2, 5)) unmodel = np.zeros((self.trace_num*2, 5)) self.harmonic_trans = np.zeros((5, self.nsamp)) self.unmodel_trans = np.zeros((5, self.nsamp)) for i, baz in enumerate(self.bazi): harmonic[i] = np.array([1, cosd(baz), sind(baz), cosd(baz*2), sind(baz*2)]) harmonic[i+self.trace_num] = np.array([0, cosd(baz+90), sind(baz+90), cosd(baz*2+90), sind(baz*2+90)]) unmodel[i] = np.array([1, cosd(baz), sind(baz), cosd(baz*2), sind(baz*2)]) unmodel[i+self.trace_num] = np.array([0, cosd(baz-90), sind(baz-90), cosd(baz*2-90), sind(baz*2-90)]) for i in range(self.nsamp): b = self.traces[:, i] self.harmonic_trans[:, i] = lsqr(harmonic, b, damp=0)[0] self.unmodel_trans[:, i] = lsqr(unmodel, b, damp=0)[0]
[docs] def write_constant(self, out_sac_path='./'): """ Write constant component to SAC file. :param out_sac_path: Output path, defaults to './' :type out_sac_path: str, optional """ sac = SACTrace(data=self.harmonic_trans[0, :]) sac.b = self.tmin sac.delta = self.rfsta.sampling sac.user0 = 0.0 sac.user1 = self.rfsta.f0[0] sac.stla = self.rfsta.stla sac.stlo = self.rfsta.stlo sac.stel = self.rfsta.stel sac.write(join(out_sac_path, '{}_constant_R.sac'.format(self.rfsta.staname)))
[docs] def plot(self, outpath='./', enf=2.): """Plot harmonic and unmodeled components :param outpath: Output path, defaults to './' :type outpath: str, optional :param enf: Amplification factor, defaults to 2.0 :type enf: float, optional """ plt.style.use("bmh") plt.rc('grid', color='white', linestyle='-', linewidth=0.7) plt.rcParams["axes.grid.axis"] = "x" fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True) mtx = [self.harmonic_trans, self.unmodel_trans] bound = np.zeros(self.nsamp) titles = ['Dipping/Anisotropic', 'Unmodeled'] for iax, ax in enumerate(axes): for i in range(5): data = mtx[iax][5-i-1] * enf + (i + 1) ax.fill_between(self.time_axis, data, bound + i+1, where=data > i+1, facecolor='red', alpha=0.7) ax.fill_between(self.time_axis, data, bound + i+1, where=data < i+1, facecolor='#1193F4', alpha=0.7) ax.plot([0, 0], [0, 6], color='k', lw=1.2) ax.set_xlim([-1, self.tmax]) ax.set_xlabel('Time after P (s)', fontsize=13) ax.set_ylim([0, 6]) ax.set_title(titles[iax]) axes[0].set_yticks([1, 2, 3, 4, 5]) axes[0].set_yticklabels(['sin2${\\theta}$', 'cos2${\\theta}$', 'sin${\\theta}$', 'cos${\\theta}$', 'Constant'], fontsize=13) fig.savefig(join(outpath, '{}_harmonic_trans.png'.format(self.rfsta.staname)), dpi=300, bbox_inches='tight')