提交 be581fa5 编辑于 作者: Mijian Xu's avatar Mijian Xu 😷
浏览文件

pick UI for S wave

上级 a37e4b41
import numpy as np
from numpy.core.getlimits import _KNOWN_TYPES
from numpy.lib.arraysetops import isin
import obspy
from obspy.io.sac import SACTrace
from obspy.signal.rotate import rotate2zne, rotate_zne_lqt
from scipy.signal import resample
from os.path import dirname, join, expanduser
import seispy
from seispy.decon import deconvolute
from seispy.geo import snr, srad2skm, rotateSeisENtoTR, skm2srad, \
rssq, extrema
from obspy.signal.trigger import recursive_sta_lta
def rotateZNE(st):
try:
......@@ -37,6 +42,7 @@ class eq(object):
self.PRaypara = None
self.SArrival = None
self.SRaypara = None
self.trigger_shift = 0
def channel_correct(self, switchEN=False, reverseE=False, reverseN=False):
if reverseE:
......@@ -97,10 +103,10 @@ class eq(object):
bazs = np.arange(-offset, offset) + bazi
ampt = np.zeros_like(bazs)
for i, b in enumerate(bazs):
t, _ = seispy.geo.rotateSeisENtoTR(this_st[0].data, this_st[1].data, b)
ampt[i] = seispy.geo.rssq(t)
t, _ = rotateSeisENtoTR(this_st[0].data, this_st[1].data, b)
ampt[i] = rssq(t)
ampt = ampt / np.max(ampt)
idx = seispy.geo.extrema(ampt, opt='min')
idx = extrema(ampt, opt='min')
if len(idx) == 0:
corr_baz = np.nan
elif len(idx) > 1:
......@@ -129,9 +135,6 @@ class eq(object):
if phase not in ('P', 'S'):
raise ValueError('phase must be in [\'P\', \'S\']')
if self.rf == obspy.Stream():
raise ValueError('Please cut out 3 components before')
if inc is None:
if self.PArrival is None or self.SArrival is None:
raise ValueError('inc must be specified')
......@@ -150,14 +153,14 @@ class eq(object):
if method == 'NE->RT':
self.comp = 'rtz'
self.rf.rotate('NE->RT', back_azimuth=bazi)
self.st.rotate('NE->RT', back_azimuth=bazi)
elif method == 'RT->NE':
self.rf.rotate('RT->NE', back_azimuth=bazi)
self.st.rotate('RT->NE', back_azimuth=bazi)
elif method == 'ZNE->LQT':
self.comp = 'lqt'
self.rf.rotate('ZNE->LQT', back_azimuth=bazi, inclination=inc)
self.st.rotate('ZNE->LQT', back_azimuth=bazi, inclination=inc)
elif method == 'LQT->ZNE':
self.rf.rotate('LQT->ZNE', back_azimuth=bazi, inclination=inc)
self.st.rotate('LQT->ZNE', back_azimuth=bazi, inclination=inc)
else:
pass
......@@ -165,15 +168,15 @@ class eq(object):
st_noise = self.trim(length, 0, phase=phase, isreturn=True)
st_signal = self.trim(0, length, phase=phase, isreturn=True)
try:
snr_E = seispy.geo.snr(st_signal[0].data, st_noise[0].data)
snr_E = snr(st_signal[0].data, st_noise[0].data)
except IndexError:
snr_E = 0
try:
snr_N = seispy.geo.snr(st_signal[1].data, st_noise[1].data)
snr_N = snr(st_signal[1].data, st_noise[1].data)
except IndexError:
snr_N = 0
try:
snr_Z = seispy.geo.snr(st_signal[2].data, st_noise[2].data)
snr_Z = snr(st_signal[2].data, st_noise[2].data)
except IndexError:
snr_Z = 0
return snr_E, snr_N, snr_Z
......@@ -204,23 +207,39 @@ class eq(object):
return Pcorrect_time, Scorrect_time
def trim(self, time_before, time_after, phase='P', isreturn=False):
"""
offset = sac.b - real o
"""
def _get_time(self, time_before, time_after, phase='P'):
if phase not in ['P', 'S']:
raise ValueError('Phase must in \'P\' or \'S\'')
P_arr, S_arr = self.arr_correct(write_to_sac=False)
time_dict = dict(zip(['P', 'S'], [P_arr, S_arr]))
t1 = self.st[2].stats.starttime + (time_dict[phase] - time_before)
t2 = self.st[2].stats.starttime + (time_dict[phase] + time_after)
if isreturn:
return self.st.copy().trim(t1, t2)
t1 = self.st[2].stats.starttime + (time_dict[phase] + self.trigger_shift - time_before)
t2 = self.st[2].stats.starttime + (time_dict[phase] + self.trigger_shift + time_after)
return t1, t2
def phase_trigger(self, time_before, time_after, phase='S', stl=5, ltl=10):
t1, t2 = self._get_time(time_before, time_after, phase)
self.t1_pick, self.t2_pick = t1, t2
self.st_pick = self.st.copy().trim(t1, t2)
if phase == 'P':
tr = self.st_pick.select(channel='*Z')[0]
else:
self.rf = self.st.copy().trim(t1, t2)
tr = self.st_pick.select(channel='*T')[0]
df = tr.stats.sampling_rate
cft = recursive_sta_lta(tr.data, int(stl*df), int(ltl*df))
n_trigger = np.argmax(np.diff(cft)[int(ltl*df):])+int(ltl*df)
self.t_trigger = t1 + n_trigger/df
self.trigger_shift = n_trigger/df - time_before
def trim(self, time_before, time_after, phase='P'):
"""
offset = sac.b - real o
"""
t1, t2 = self._get_time(time_before, time_after, phase)
self.st.trim(t1, t2)
def deconvolute(self, shift, time_after, f0=2, phase='P', method='iter', only_r=False, itmax=400, minderr=0.001, wlevel=0.05, target_dt=None):
self.rf = self.st.copy()
if phase not in ['P', 'S']:
raise ValueError('Phase must in \'P\' or \'S\'')
if self.rf == obspy.Stream():
......@@ -240,13 +259,13 @@ class eq(object):
raise ValueError('method must be in \'iter\' or \'water\'')
if phase == 'P':
decon_out_r = seispy.decon.deconvolute(self.rf[1].data, self.rf[2].data, self.rf[1].stats.delta, **kwargs)
decon_out_r = deconvolute(self.rf[1].data, self.rf[2].data, self.rf[1].stats.delta, **kwargs)
self.rf[1].data = decon_out_r[0]
self.rms = decon_out_r[1]
if method == 'iter':
self.it = decon_out_r[2]
if not only_r:
decon_out_t = seispy.decon.deconvolute(self.rf[0].data, self.rf[2].data, self.rf[1].stats.delta, **kwargs)
decon_out_t = deconvolute(self.rf[0].data, self.rf[2].data, self.rf[1].stats.delta, **kwargs)
self.rf[0].data = decon_out_t[0]
else:
# TODO: if 'Q' not in self.rf[1].stats.channel or 'L' not in self.rf[2].stats.channel:
......@@ -271,8 +290,8 @@ class eq(object):
wdat = -win.data
udat = uin.data
wdat[0:int((tshift-4)/win.stats.delta)] = 0
srf, self.rms, self.it = seispy.decon.deconvolute(udat, wdat, win.stats.delta,
phase='S', tshift=tshift, **kwargs)
srf, self.rms, self.it = deconvolute(udat, wdat, win.stats.delta,
phase='S', tshift=tshift, **kwargs)
uin.data = np.flip(srf)
def saverf(self, path, evtstr=None, phase='P', shift=0, evla=-12345., evlo=-12345., evdp=-12345., mag=-12345.,
......@@ -282,13 +301,13 @@ class eq(object):
loop_lst = ['R']
else:
loop_lst = ['R', 'T']
rayp = seispy.geo.srad2skm(self.PArrival.ray_param)
rayp = srad2skm(self.PArrival.ray_param)
elif phase == 'S':
if self.comp == 'lqt':
loop_lst = ['L']
else:
loop_lst = ['Z']
rayp = seispy.geo.srad2skm(self.SArrival.ray_param)
rayp = srad2skm(self.SArrival.ray_param)
else:
raise ValueError('Phase must be in \'P\' or \'S\'')
......@@ -330,7 +349,7 @@ class eq(object):
# Final RMS
if rmsgate is not None:
if type(self.rms) == np.ndarray:
if isin(self.rms, np.ndarray):
rms = self.rms[-1]
else:
rms = self.rms
......
......@@ -413,9 +413,14 @@ class RF(object):
for _, row in self.eqs.iterrows():
row['data'].trim(self.para.time_before, self.para.time_after, self.para.phase)
def pick(self):
pickphase(self.eqs, self.para)
def pick(self, prepick=True, stl=5, ltl=10):
if prepick:
self.logger.RFlog.info('Pre-pick {} arrival using STA/LTA method'.format(self.para.phase))
for _, row in self.eqs.iterrows():
row['data'].phase_trigger(self.para.time_before, self.para.time_after,
self.para.phase, stl=stl, ltl=ltl)
self.logger.RFlog.info('{0} events left after virtual checking'.format(self.eqs.shape[0]))
pickphase(self.eqs, self.para)
def deconv(self):
shift = self.para.time_before
......@@ -474,7 +479,7 @@ def prf():
metavar='N|E|NE', default=None, type=str)
parser.add_argument('-s', help='Switch the East and North components', dest='isswitch', action='store_true')
parser.add_argument('-b', help='Correct back-azimuth. \nIf "baz" is specified, the corr_baz = raw_baz + baz. \n'
'If no arguments following -b Back-azimuth will be corrected with minimal '
'If there is no argument, the back-azimuth will be corrected with minimal '
'energy of T component. The searching range is raw_baz +/- 90',
dest='baz', nargs='?', const=0, type=float)
arg = parser.parse_args()
......@@ -514,8 +519,8 @@ def prf():
pjt.baz_correct()
else:
pass
pjt.trim()
pjt.rotate()
pjt.trim()
pjt.deconv()
pjt.saverf()
......
import sys
import os
import argparse
# matplotlib.use("Qt5Agg")
from PyQt5.QtCore import Qt
from PyQt5.QtGui import QIcon, QKeySequence
from PyQt5.QtWidgets import QApplication, QMainWindow, QVBoxLayout, \
QSizePolicy, QWidget, QDesktopWidget, \
......@@ -12,6 +10,8 @@ from seispy.setuplog import setuplog
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.figure import Figure
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
......@@ -19,6 +19,8 @@ class SFigure(Figure):
def __init__(self, eqs, para, width=21, height=11, dpi=100):
self.eqs = eqs
self.row_num = eqs.shape[0]
self.picker_time = pd.DataFrame({'trigger_shift': np.array([eqs.iloc[i]['data'].trigger_shift for i in range(self.row_num)])})
self.picker_time.set_index(eqs.index, inplace=True)
self.para = para
self.log = setuplog()
self.idx = 0
......@@ -29,6 +31,14 @@ class SFigure(Figure):
def init_figure(self, width=21, height=11, dpi=100):
self.fig, self.axs = plt.subplots(3, 1, sharex=True, figsize=(width, height), dpi=dpi, tight_layout=True)
def set_cross_hair_visible(self, visible):
need_redraw = self.cursor[0].get_visible() != visible
# self.horizontal_line.set_visible(visible)
for cursor in self.cursor:
cursor.set_visible(visible)
# self.text.set_visible(visible)
return need_redraw
def set_properties(self):
date = self.eqs.iloc[self.idx]['date'].strftime("%Y.%m.%dT%H:%M:%S")
......@@ -39,6 +49,7 @@ class SFigure(Figure):
self.fig.suptitle(title)
for ax in self.axs:
ax.set(xlim=[-self.para.time_before, self.para.time_after])
ax.autoscale(enable=True,axis='y', tight=True)
ax.minorticks_on()
self.axs[2].set_xlabel('Time (s)')
......@@ -47,11 +58,15 @@ class SFigure(Figure):
ax.cla()
def plot(self, lc='tab:blue'):
st = self.eqs.iloc[self.idx]['data'].rf
st = self.eqs.iloc[self.idx]['data'].st_pick
# t1, t2 = self.eqs.iloc[self.idx]['data'].t1_pick, st = self.eqs.iloc[self.idx]['data'].t2_pick
self.time_axis = st[0].times()-self.para.time_before
self.picker = []
for i in range(3):
self.axs[i].plot(self.time_axis, st[i].data, color=lc)
self.axs[i].axvline(x=0, lw=1, ls='--', color='r')
self.picker.append(self.axs[i].axvline(self.picker_time.iloc[self.idx]['trigger_shift'], color='k', lw=1.5, visible=True))
# self.axs[i].axvline(x=trigger, lw=1, alpha=0.3, color='g')
self.axs[i].axvline(x=0, lw=1.5, alpha=0.3, color='r')
self.axs[i].set_ylabel(st[i].stats.channel)
self.set_properties()
......@@ -64,6 +79,31 @@ class SFigure(Figure):
self.plot(lc=self.drop_color)
else:
self.plot()
def on_mouse_move(self, event):
if not event.inaxes:
need_redraw = self.set_cross_hair_visible(False)
if need_redraw:
for ax in self.axs:
ax.figure.canvas.draw()
else:
self.set_cross_hair_visible(True)
x = event.xdata
# update the line positions
# self.horizontal_line.set_ydata(y)
for cursor in self.cursor:
cursor.set_xdata(x)
# ax.figure.canvas.draw()
# self.cursor[1].set_xdata(x)
# self.axs[1].figure.canvas.draw()
def on_click(self, event):
if not event.inaxes:
return
self.picker_time.iloc[self.idx]['trigger_shift'] = event.xdata
for picker in self.picker:
picker.set_xdata(self.picker_time.iloc[self.idx]['trigger_shift'])
picker.set_visible(True)
def back_action(self):
self.clear()
......@@ -89,6 +129,8 @@ class SFigure(Figure):
def finish(self):
self.eqs.drop(self.drop_lst, inplace=True)
for i, row in self.eqs.iterrows():
row['data'].trigger_shift = self.picker_time['trigger_shift'][i]
class MyMplCanvas(FigureCanvas):
def __init__(self, parent=None, eqs=None, para=None, width=21, height=11, dpi=100):
......@@ -111,6 +153,9 @@ class MatplotlibWidget(QMainWindow):
def initUi(self, eqs, para):
self.mpl = MyMplCanvas(self, eqs=eqs, para=para, width=21, height=11, dpi=100)
# self.mpl.mpl_connect('motion_notify_event', self.on_mouse_move)
self.mpl.mpl_connect('button_press_event', self.on_click)
self.setCursor(Qt.CrossCursor)
self.layout = QVBoxLayout()
self.layout.addWidget(self.mpl, 2)
......@@ -126,6 +171,14 @@ class MatplotlibWidget(QMainWindow):
def exit_app(self):
self.close()
def on_click(self, event):
self.mpl.sfig.on_click(event)
self.mpl.draw()
def on_mouse_move(self, event):
self.mpl.sfig.on_mouse_move(event)
self.mpl.draw()
def next_connect(self):
self.mpl.sfig.next_action()
self.mpl.draw()
......@@ -144,7 +197,8 @@ class MatplotlibWidget(QMainWindow):
def finish_connect(self):
self.mpl.sfig.finish()
QApplication.quit()
self.close()
# QApplication.quit()
def _define_global_shortcuts(self):
self.key_c = QShortcut(QKeySequence('c'), self)
......
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册