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

add RFTrace for decon

上级 791a09fa
......@@ -138,7 +138,7 @@ class CCPProfile():
self.idxs.append(np.where(distaz(bin_info[0], bin_info[1], self.stalst[:, 0], self.stalst[:, 1]).delta <= dis)[0])
elif self.cpara.width is not None and self.cpara.shape == 'rect':
self.logger.CCPlog.info('Select stations within {} km perpendicular to the profile'.format(self.cpara.width))
self.idxs = self._proj_sta()
self.idxs = self._proj_sta(self.cpara.width)
self._write_sta()
else:
if self.cpara.width is not None:
......
......@@ -49,10 +49,9 @@ def phaseshift(x, nfft, dt, tshift):
return x
def deconit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001, info=False, phase='P'):
def deconit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001, phase='P'):
"""
Created on Wed Sep 10 14:21:38 2014
[RFI, rms, it]=makeRFitdecon(uin,win,dt,nt,tshift,f0,itmax,minderr)
In:
uin = numerator (radial for PdS)
......@@ -213,61 +212,7 @@ def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False, p
gnorm = np.sum(gaussF) * delf * dt
rft = rft.real / gnorm
return rft, rms, np.nan
def _add_zeros(a, numl, numr):
"""Add zeros at left and rigth side of array a"""
return np.hstack([np.zeros(numl), a, np.zeros(numr)])
def _acorrt(a, num):
"""
Not normalized auto-correlation of signal a.
Sample 0 corresponds to zero lag time. Auto-correlation will consist of
num samples. Correlation is performed in time domain with scipy.
:param a: Data
:param num: Number of returned data points
:return: autocorrelation
"""
return correlate(_add_zeros(a, 0, num - 1), a, 'valid')
def _xcorrt(a, b, num, zero_sample=0):
"""
Not normalized cross-correlation of signals a and b.
:param a,b: data
:param num: The cross-correlation will consist of num samples.\n
The sample with 0 lag time will be in the middle.
:param zero_sample: Signals a and b are aligned around the middle of their
signals.\n
If zero_sample != 0 a will be shifted additionally to the left.
:return: cross-correlation
"""
if zero_sample > 0:
a = _add_zeros(a, 2 * abs(zero_sample), 0)
elif zero_sample < 0:
a = _add_zeros(a, 0, 2 * abs(zero_sample))
dif = len(a) - len(b) + 1 - num
if dif > 0:
b = _add_zeros(b, (dif + 1) // 2, dif // 2)
else:
a = _add_zeros(a, (-dif + 1) // 2, (-dif) // 2)
return correlate(a, b, 'valid')
def decontime(uin, win, dt, tshift=10, spiking=1., normalize=True):
nt = uin.size
nshift = int(tshift * dt)
STS = _acorrt(win, nt)
STS = STS / STS[0]
STS[0] += spiking
STR = _xcorrt(uin, win, nt, nshift)
assert len(STR) == len(STS)
rf = solve_toeplitz(STS, STR)
if normalize:
norm = 1 / np.max(np.abs(rf))
rf *= norm
return rf
return rft, rms
def deconvolute(uin, win, dt, method='iter', **kwargs):
......@@ -276,7 +221,30 @@ def deconvolute(uin, win, dt, method='iter', **kwargs):
elif method.lower() == 'water':
return deconwater(uin, win, dt, **kwargs)
else:
raise ValueError('method must be in the \'iter\' or \'water\'')
raise ValueError('method must be \'iter\' or \'water\'')
class RFTrace(obspy.Trace):
def __init__(self, data=..., header=None):
super().__init__(data=data, header=header)
@classmethod
def deconvolute(cls, utr, wtr, method='iter', **kwargs):
header = utr.stats.__getstate__()
for key, value in kwargs.items():
header[key] = value
if method.lower() == 'iter':
rf, rms, it = deconit(utr.data, wtr.data, utr.stats.delta, **kwargs)
header['rms'] = rms
header['iter'] = it
elif method.lower() == 'water':
rf, rms = deconwater(utr.data, wtr.data, utr.stats.delta, **kwargs)
header['rms'] = rms
header['iter'] = np.nan
else:
raise ValueError('method must be \'iter\' or \'water\'')
return cls(rf, header)
if __name__ == '__main__':
......
......@@ -6,7 +6,7 @@ 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
from seispy.decon import deconvolute
from seispy.decon import RFTrace
from seispy.geo import snr, srad2skm, rotateSeisENtoTR, skm2srad, \
rssq, extrema
from obspy.signal.trigger import recursive_sta_lta
......@@ -267,14 +267,12 @@ class eq(object):
raise ValueError('method must be in \'iter\' or \'water\'')
if phase == 'P':
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]
self.rf[1] = RFTrace.deconvolute(self.rf[1], self.rf[2], **kwargs)
self.rms = self.rf[1].stats.rms
if method == 'iter':
self.it = decon_out_r[2]
self.it = self.rf[1].stats.iter
if not only_r:
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]
self.rf[0] = RFTrace.deconvolute(self.rf[0], self.rf[2], **kwargs)
else:
# TODO: if 'Q' not in self.rf[1].stats.channel or 'L' not in self.rf[2].stats.channel:
# raise ValueError('Please rotate component to \'LQT\'')
......@@ -291,16 +289,15 @@ class eq(object):
if self.comp == 'lqt':
win = self.rf.select(channel='*Q')[0]
uin = self.rf.select(channel='*L')[0]
wdat = win.data
else:
win = self.rf.select(channel='*R')[0]
uin = self.rf.select(channel='*Z')[0]
wdat = -win.data
udat = uin.data
wdat[0:int((tshift-4)/win.stats.delta)] = 0
srf, self.rms, self.it = deconvolute(udat, wdat, win.stats.delta,
phase='S', tshift=tshift, **kwargs)
uin.data = np.flip(srf)
win.data *= -1
win.data[0:int((tshift-4)/win.stats.delta)] = 0
uin = RFTrace.deconvolute(uin, win, phase='S', tshift=tshift, **kwargs)
self.rms = uin.stats.rms
self.it = uin.stats.iter
uin.data = np.flip(uin.data)
def saverf(self, path, evtstr=None, phase='P', shift=0, evla=-12345., evlo=-12345., evdp=-12345., mag=-12345.,
gauss=0, baz=-12345., gcarc=-12345., only_r=False, **kwargs):
......
import numpy as np
from numpy import pi, mod
from seispy import distaz
from scipy import interpolate
def sind(deg):
rad = np.radians(deg)
......@@ -82,12 +80,15 @@ def skm2srad(skm):
def rot3D(bazi, inc):
"""
:param bazi:
:param inc:
:return:
Rotate components from ZRT to LQT
M = [cos(inc) -sin(inc)*sin(bazi) -sin(inc)*cos(bazi);
sin(inc) cos(inc)*sin(bazi) cos(inc)*cos(bazi);
0 -cos(bazi) sin(bazi)];
:param bazi: back-azimuth of station-event pair
:param inc: Incidence angle
:return: Coefficient matrix M
"""
if isinstance(inc, float) or isinstance(inc, int):
......@@ -181,6 +182,6 @@ def extrema(x, opt='max'):
elif opt == 'min':
idx = np.intersect1d(np.where(np.diff(x) < 0)[0]+1, np.where(np.diff(x) > 0)[0])
else:
raise ImportError('Wrong Options!!!')
raise ValueError('opt must be \'max\' or \'min\'')
return idx
......@@ -50,9 +50,9 @@ class RFFigure(Figure):
self.maxidx = 20
self.xlim = xlim
def init_canvas(self):
def init_canvas(self, order='baz'):
self.init_figure(width=self.width, height=self.height, dpi=self.dpi)
self.read_sac()
self.read_sac(order=order)
self.set_figure()
self.set_page()
self.init_variables()
......@@ -116,7 +116,13 @@ class RFFigure(Figure):
self.axb.set_xlabel("Backazimuth (\N{DEGREE SIGN})")
self.axg.set_xlabel('Distance (\N{DEGREE SIGN})')
def read_sac(self, dt=0.1):
def read_sac(self, dt=0.1, order='baz'):
if not isinstance(order, str):
raise TypeError('The order must be str type')
elif not order in ['baz', 'dis']:
raise ValueError('The order must be \'baz\' or \'dis\'')
else:
pass
if len(glob.glob(join(self.rfpath, '*_R.sac'))) != 0:
tmp_files = glob.glob(join(self.rfpath, '*_R.sac'))
self.comp = 'R'
......@@ -134,18 +140,24 @@ class RFFigure(Figure):
self.evt_num = len(self.rrf)
self.log.RFlog.info('A total of {} PRFs loaded'.format(self.evt_num))
self.baz = np.array([tr.stats.sac.baz for tr in self.rrf])
self.sort_baz_()
self.gcarc = np.array([tr.stats.sac.gcarc for tr in self.rrf])
self._sort(order)
self.axpages, self.rfidx = indexpags(self.evt_num, self.maxidx)
self.staname = (self.rrf[0].stats.network+'.'+self.rrf[0].stats.station).strip('.')
self.fig.suptitle("%s (Latitude: %5.2f\N{DEGREE SIGN}, Longitude: %5.2f\N{DEGREE SIGN})" % (self.staname, self.rrf[0].stats.sac.stla, self.rrf[0].stats.sac.stlo), fontsize=20)
def sort_baz_(self):
idx = np.argsort(self.baz)
def _sort(self, order):
if order == 'baz':
idx = np.argsort(self.baz)
elif order == 'dis':
idx = np.argsort(self.gcarc)
else:
pass
self.baz = self.baz[idx]
self.gcarc = self.gcarc[idx]
self.rrf = [self.rrf[i] for i in idx]
self.trf = [self.trf[i] for i in idx]
self.gcarc = [self.rrf[i].stats.sac.gcarc for i in range(self.evt_num)]
# self.gcarc = [self.rrf[i].stats.sac.gcarc for i in range(self.evt_num)]
self.filenames = [self.filenames[i] for i in idx]
def plotwave(self):
......
......@@ -17,16 +17,17 @@ import glob
class MyMplCanvas(FigureCanvas):
def __init__(self, parent=None, rfpath='', only_r=False, width=21, height=11, dpi=100, xlim=[-2, 30]):
def __init__(self, parent=None, rfpath='', only_r=False, width=21, height=11,
dpi=100, xlim=[-2, 30], order='baz'):
plt.rcParams['axes.unicode_minus'] = False
if only_r:
self.rffig = RPickFigure(rfpath, width=width, height=height, dpi=dpi, xlim=xlim)
self.rffig.init_canvas()
self.rffig.init_canvas(order=order)
else:
self.rffig = RFFigure(rfpath, width=width, height=height, dpi=dpi, xlim=xlim)
self.rffig.init_canvas()
self.rffig.init_canvas(order=order)
FigureCanvas.__init__(self, self.rffig.fig)
self.setParent(parent)
......@@ -38,15 +39,17 @@ class MyMplCanvas(FigureCanvas):
class MatplotlibWidget(QMainWindow):
def __init__(self, rfpath, only_r=False, xlim=[-2, 30], parent=None):
def __init__(self, rfpath, only_r=False, xlim=[-2, 30],
order='baz', parent=None):
super(MatplotlibWidget, self).__init__(parent)
self.only_r = only_r
self.initUi(rfpath, only_r, xlim)
self.initUi(rfpath, only_r, xlim, order=order)
def initUi(self, rfpath, only_r, xlim):
def initUi(self, rfpath, only_r, xlim, order='baz'):
self.layout = QVBoxLayout()
self.add_btn()
self.mpl = MyMplCanvas(self, rfpath=rfpath, only_r=only_r, width=21, height=11, dpi=100, xlim=xlim)
self.mpl = MyMplCanvas(self, rfpath=rfpath, only_r=only_r, width=21, height=11,
dpi=100, xlim=xlim, order=order)
self.layout.addWidget(self.mpl, 2)
self.mpl.mpl_connect('button_press_event', self.on_click)
......@@ -168,7 +171,9 @@ class MatplotlibWidget(QMainWindow):
def main():
parser = argparse.ArgumentParser(description="User interface for picking PRFs")
parser.add_argument('rf_path', type=str, help='Path to PRFs')
parser.add_argument('-x', help="Set the x limits of the current axes, defaults to 30s for RT, 85s for R.",
parser.add_argument('-a', help='Arrangement of RFs, defaults to \'baz\'', dest='order',
default='baz', type=str, metavar='baz|dis')
parser.add_argument('-x', help="Set x limits of the current axes, defaults to 30s for RT, 85s for R.",
dest='xlim', default=None, type=float, metavar='xmax')
arg = parser.parse_args()
rfpath = arg.rf_path
......@@ -183,7 +188,7 @@ def main():
if arg.xlim is not None:
xlim = arg.xlim
app = QApplication(sys.argv)
ui = MatplotlibWidget(rfpath, only_r=only_r, xlim=[-2, xlim])
ui = MatplotlibWidget(rfpath, only_r=only_r, xlim=[-2, xlim], order=arg.order)
ui.show()
sys.exit(app.exec_())
......
import glob
import os
import obspy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.figure import Figure
from seispy.pickfigure import RFFigure, indexpags
from os.path import join, basename
from seispy.setuplog import setuplog
from seispy.plotR import init_figure, set_fig, plot_waves
......@@ -29,8 +27,8 @@ class RPickFigure(RFFigure):
super().__init__(rfpath, width=width, height=height, dpi=dpi, xlim=xlim)
self.enf = 8
def init_canvas(self):
return super().init_canvas()
def init_canvas(self, order='baz'):
return super().init_canvas(order=order)
def init_figure(self, width, height, dpi):
self.fig = Figure(figsize=(width, height), dpi=dpi)
......@@ -38,7 +36,13 @@ class RPickFigure(RFFigure):
self.axb = self.fig.add_axes([0.755, 0.05, 0.22, 0.84])
self.axg = self.axb.twiny()
def read_sac(self, dt=0.1):
def read_sac(self, dt=0.1, order='baz'):
if not isinstance(order, str):
raise TypeError('The order must be str type')
elif not order in ['baz', 'dis']:
raise ValueError('The order must be \'baz\' or \'dis\'')
else:
pass
if len(glob.glob(join(self.rfpath, '*_R.sac'))) != 0:
tmp_files = glob.glob(join(self.rfpath, '*_R.sac'))
self.comp = 'R'
......@@ -55,17 +59,23 @@ class RPickFigure(RFFigure):
self.evt_num = len(self.rrf)
self.log.RFlog.info('A total of {} PRFs loaded'.format(self.evt_num))
self.baz = np.array([tr.stats.sac.baz for tr in self.rrf])
self.sort_baz_()
self.gcarc = np.array([tr.stats.sac.gcarc for tr in self.rrf])
self._sort(order)
self.axpages, self.rfidx = indexpags(self.evt_num, self.maxidx)
self.staname = (self.rrf[0].stats.network+'.'+self.rrf[0].stats.station).strip('.')
self.fig.suptitle("%s (Latitude: %5.2f\N{DEGREE SIGN}, Longitude: %5.2f\N{DEGREE SIGN})" % (self.staname, self.rrf[0].stats.sac.stla, self.rrf[0].stats.sac.stlo), fontsize=20)
def sort_baz_(self):
idx = np.argsort(self.baz)
def _sort(self, order):
if order == 'baz':
idx = np.argsort(self.baz)
elif order == 'dis':
idx = np.argsort(self.gcarc)
else:
pass
self.baz = self.baz[idx]
self.gcarc = self.gcarc[idx]
self.rrf = [self.rrf[i] for i in idx]
self.gcarc = [self.rrf[i].stats.sac.gcarc for i in range(self.evt_num)]
# self.gcarc = [self.rrf[i].stats.sac.gcarc for i in range(self.evt_num)]
self.filenames = [self.filenames[i] for i in idx]
def set_figure(self):
......
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册