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

srf and saving config

上级 2429d6b5
......@@ -136,7 +136,7 @@ def deconit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001,
return RFI, rms, it
def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False):
def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False, phase='P'):
"""
Frequency-domain deconvolution using waterlevel method.
......@@ -213,7 +213,7 @@ def deconwater(uin, win, dt, tshift=10., wlevel=0.05, f0=2.0, normalize=False):
gnorm = np.sum(gaussF) * delf * dt
rft = rft.real / gnorm
return rft, rms
return rft, rms, np.nan
def _add_zeros(a, numl, numr):
"""Add zeros at left and rigth side of array a"""
......
......@@ -21,22 +21,13 @@ def rotateZNE(st):
class eq(object):
def __init__(self, pathname, datestr, suffix, switchEN=False, reverseE=False, reverseN=False):
def __init__(self, pathname, datestr, suffix):
self.datestr = datestr
self.st = obspy.read(join(pathname, '*' + datestr + '*' + suffix))
if len(self.st) != 3:
raise ValueError('Sismogram must be in 3 components, but there are {} of {}'.format(len(self.st), datestr))
# if not (self.st[0].stats.npts == self.st[1].stats.npts == self.st[2].stats.npts):
# raise ValueError('Samples are different in 3 components')
if reverseE:
self.st.select(channel='*[E2]')[0].data *= -1
if reverseN:
self.st.select(channel='*[N1]')[0].data *= -1
if switchEN:
chE = self.st.select(channel='*[E2]')[0].stats.channel
chN = self.st.select(channel='*[N1]')[0].stats.channel
self.st.select(channel='*[E2]')[0].stats.channel = chN
self.st.select(channel='*[N1]')[0].stats.channel = chE
self.st.sort()
self.rf = obspy.Stream()
self.timeoffset = 0
......@@ -47,6 +38,17 @@ class eq(object):
self.SArrival = None
self.SRaypara = None
def channel_correct(self, switchEN=False, reverseE=False, reverseN=False):
if reverseE:
self.st.select(channel='*[E2]')[0].data *= -1
if reverseN:
self.st.select(channel='*[N1]')[0].data *= -1
if switchEN:
chE = self.st.select(channel='*[E2]')[0].stats.channel
chN = self.st.select(channel='*[N1]')[0].stats.channel
self.st.select(channel='*[E2]')[0].stats.channel = chN
self.st.select(channel='*[N1]')[0].stats.channel = chE
def __str__(self):
return('Event data class {0}'.format(self.datestr))
......@@ -147,12 +149,13 @@ class eq(object):
pass
if method == 'NE->RT':
self.comp = 'rtz'
self.rf.rotate('NE->RT', back_azimuth=bazi)
elif method == 'RT->NE':
self.rf.rotate('RT->NE', back_azimuth=bazi)
elif method == 'ZNE->LQT':
self.comp = 'lqt'
self.rf.rotate('ZNE->LQT', back_azimuth=bazi, inclination=inc)
self.rf[1].data *= -1
elif method == 'LQT->ZNE':
self.rf.rotate('LQT->ZNE', back_azimuth=bazi, inclination=inc)
else:
......@@ -248,17 +251,7 @@ class eq(object):
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\'')
Q = np.flip(self.rf[1].data, axis=0)
L = np.flip(self.rf[2].data, axis=0)
self.rf[2].data, self.rms, self.it = seispy.decon.deconit(L, -Q, self.rf[1].stats.delta,
self.rf[1].data.shape[0], shift,
f0, itmax, minderr, phase='S')
# self.rf[2].data = seispy.decon.deconv_waterlevel(L, -Q, self.rf[1].stats.delta,
# self.rf[1].data.shape[0], tshift=shift,
# gauss=f0)
# self.rms = [0]
# self.it = 0
# self.rf[2].data = - self.rf[2].data
self.decon_s(**kwargs)
if target_dt is not None:
if self.rf[0].stats.delta != target_dt:
# self.rf.resample(1 / target_dt)
......@@ -267,16 +260,34 @@ class eq(object):
tr.data = resample(tr.data, int((shift + time_after)/target_dt+1))
tr.stats.delta = target_dt
def decon_s(self, tshift, **kwargs):
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 = seispy.decon.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.,
gauss=0, baz=-12345., gcarc=-12345., only_r=False, **kwargs):
if phase == 'P':
if only_r:
loop_lst = [1]
loop_lst = ['R']
else:
loop_lst = [0, 1]
loop_lst = ['R', 'T']
rayp = seispy.geo.srad2skm(self.PArrival.ray_param)
elif phase == 'S':
loop_lst = [2]
if self.comp == 'lqt':
loop_lst = ['L']
else:
loop_lst = ['Z']
rayp = seispy.geo.srad2skm(self.SArrival.ray_param)
else:
raise ValueError('Phase must be in \'P\' or \'S\'')
......@@ -285,19 +296,20 @@ class eq(object):
filename = join(path, self.datestr)
else:
filename = join(path, evtstr)
for i in loop_lst:
for comp in loop_lst:
trrf = self.rf.select(channel='*'+comp)[0]
header = {'evla': evla, 'evlo': evlo, 'evdp': evdp, 'mag': mag, 'baz': baz,
'gcarc': gcarc, 'user0': rayp, 'kuser0': 'Ray Para', 'user1': gauss, 'kuser1': 'G factor'}
for key in kwargs:
header[key] = kwargs[key]
for key, value in header.items():
self.rf[i].stats['sac'][key] = value
tr = SACTrace.from_obspy_trace(self.rf[i])
trrf.stats['sac'][key] = value
tr = SACTrace.from_obspy_trace(trrf)
tr.b = -shift
tr.o = 0
tr.write(filename + '_{0}_{1}.sac'.format(phase, tr.kcmpnm[-1]))
def judge_rf(self, shift, npts, criterion='crust', rmsgate=None):
def judge_rf(self, shift, npts, phase='P', criterion='crust', rmsgate=None):
if not isinstance(criterion, (str, type(None))):
raise TypeError('criterion should be string in [\'crust\', \'mtz\']')
elif criterion is None:
......@@ -306,8 +318,14 @@ class eq(object):
raise ValueError('criterion should be string in [\'crust\', \'mtz\']')
else:
criterion = criterion.lower()
if self.rf[1].stats.npts != npts:
if phase == 'P':
trrf = self.rf.select(channel='*R')[0]
elif phase == 'S' and self.comp == 'lqt':
trrf = self.rf.select(channel='*L')[0]
elif phase == 'S' and self.comp == 'rtz':
trrf = self.rf.select(channel='*Z')[0]
if trrf.stats.npts != npts:
return False
# Final RMS
......@@ -321,8 +339,8 @@ class eq(object):
rmspass = True
# R energy
nt1 = int(np.floor((5+shift)/self.rf[1].stats.delta))
reng = np.sum(np.sqrt(self.rf[1].data[nt1:] ** 2))
nt1 = int(np.floor((5+shift)/trrf.stats.delta))
reng = np.sum(np.sqrt(trrf.data[nt1:] ** 2))
if reng < 0.1:
rengpass = False
else:
......@@ -330,20 +348,20 @@ class eq(object):
# Max amplitude
if criterion == 'crust':
time_P1 = int(np.floor((-2+shift)/self.rf[1].stats.delta))
time_P2 = int(np.floor((4+shift)/self.rf[1].stats.delta))
max_P = np.max(self.rf[1].data[time_P1:time_P2])
if max_P == np.max(np.abs(self.rf[1].data)) and max_P < 1 and rmspass and rengpass:
time_P1 = int(np.floor((-2+shift)/trrf.stats.delta))
time_P2 = int(np.floor((4+shift)/trrf.stats.delta))
max_P = np.max(trrf.data[time_P1:time_P2])
if max_P == np.max(np.abs(trrf.data)) and max_P < 1 and rmspass and rengpass:
return True
else:
return False
elif criterion == 'mtz':
max_deep = np.max(np.abs(self.rf[1].data[int((30 + shift) / self.rf[1].stats.delta):]))
time_P1 = int(np.floor((-5 + shift) / self.rf[1].stats.delta))
time_P2 = int(np.floor((5 + shift) / self.rf[1].stats.delta))
max_P = np.max(self.rf[1].data[time_P1:time_P2])
max_deep = np.max(np.abs(trrf.data[int((30 + shift) / trrf.stats.delta):]))
time_P1 = int(np.floor((-5 + shift) / trrf.stats.delta))
time_P2 = int(np.floor((5 + shift) / trrf.stats.delta))
max_P = np.max(trrf.data[time_P1:time_P2])
if max_deep < max_P * 0.4 and rmspass and rengpass and\
max_P == np.max(np.abs(self.rf[1].data)) and max_P < 1:
max_P == np.max(np.abs(trrf.data)) and max_P < 1:
return True
else:
return False
......
......@@ -38,6 +38,9 @@ class para(object):
self.criterion = 'crust'
self.only_r = False
self.rmsgate = None
self.switchEN=False
self.reverseE=False
self.reverseN=False
def get_para(self):
return self.__dict__
......
......@@ -31,7 +31,7 @@ class StaData():
self.ev_num = len(self.bazi)
self.datar = np.empty([self.ev_num, self.rflength])
self.datat = np.empty([self.ev_num, self.rflength])
for i, idx in enumerate(goodidx[0]):
for i, idx in enumerate(goodidx[0]):
self.datar[i, :] = rrf[idx].data
self.datat[i, :] = trf[idx].data
self.time_axis = np.array([])
......@@ -117,12 +117,15 @@ class RFFigure(Figure):
self.axg.set_xlabel('Distance (\N{DEGREE SIGN})')
def read_sac(self, dt=0.1):
if len(glob.glob(join(self.rfpath, '*_R.sac'))) == 0:
tmp_files = glob.glob(join(self.rfpath, '*_Q.sac'))
self.comp = 'Q'
else:
if len(glob.glob(join(self.rfpath, '*_R.sac'))) != 0:
tmp_files = glob.glob(join(self.rfpath, '*_R.sac'))
self.comp = 'R'
self.comp = 'R'
elif len(glob.glob(join(self.rfpath, '*_L.sac'))) != 0:
tmp_files = glob.glob(join(self.rfpath, '*_L.sac'))
self.comp = 'L'
else:
tmp_files = glob.glob(join(self.rfpath, '*_Z.sac'))
self.comp = 'Z'
self.log.RFlog.info('Reading PRFs from {}'.format(self.rfpath))
self.filenames = [basename(sac_file).split('_')[0] for sac_file in sorted(tmp_files)]
self.rrf = obspy.read(join(self.rfpath, '*_{}.sac'.format(self.comp))).sort(['starttime']).resample(1/dt)
......
from genericpath import exists
import obspy
from obspy import UTCDateTime
from obspy.io.sac import SACTrace
from obspy.taup import TauPyModel
import re
from seispy.io import wsfetch
from os.path import join
from seispy.io import wsfetch
from seispy.para import para
import seispy
from seispy import distaz
from seispy.eq import eq
from seispy.setuplog import setuplog
from seispy.sviewerui import MatplotlibWidget
import glob
import numpy as np
from datetime import timedelta
import pandas as pd
from obspy.taup import TauPyModel
import configparser
import argparse
import sys
import matplotlib.pyplot as plt
from PyQt5.QtWidgets import QApplication
from seispy.sviewerui import MatplotlibWidget
import pickle
def pickphase(eqs, para):
......@@ -53,13 +56,13 @@ def read_catalog(logpath, b_time, e_time, stla, stlo, magmin=5.5, magmax=10, dis
lines = f.readlines()
for line in lines:
line_sp = line.strip().split()
date_now = obspy.UTCDateTime.strptime('.'.join(line_sp[0:3]) + 'T' + '.'.join(line_sp[4:7]),
date_now = UTCDateTime.strptime('.'.join(line_sp[0:3]) + 'T' + '.'.join(line_sp[4:7]),
'%Y.%m.%dT%H.%M.%S')
evla = float(line_sp[7])
evlo = float(line_sp[8])
evdp = float(line_sp[9])
mw = float(line_sp[10])
dis = seispy.distaz(stla, stlo, evla, evlo).delta
dis = distaz(stla, stlo, evla, evlo).delta
# bazi = seispy.distaz(stla, stlo, evla, evlo).getBaz()
if b_time <= date_now <= e_time and magmin <= mw <= magmax and dismin <= dis <= dismax:
this_data = pd.DataFrame([[date_now, evla, evlo, evdp, mw]], columns=col)
......@@ -79,7 +82,7 @@ def load_station_info(pathname, ref_comp, suffix):
def match_eq(eq_lst, pathname, stla, stlo, ref_comp='Z', suffix='SAC', offset=None,
tolerance=210, dateformat='%Y.%j.%H.%M.%S', switchEN=False, reverseE=False, reverseN=False):
tolerance=210, dateformat='%Y.%j.%H.%M.%S'):
pattern = datestr2regex(dateformat)
ref_eqs = glob.glob(join(pathname, '*{0}*{1}'.format(ref_comp, suffix)))
if len(ref_eqs) == 0:
......@@ -91,7 +94,7 @@ def match_eq(eq_lst, pathname, stla, stlo, ref_comp='Z', suffix='SAC', offset=No
except IndexError:
raise IndexError('Error data format of {} in {}'.format(dateformat, ref_sac))
if isinstance(offset, (int, float)):
sac_files.append([datestr, obspy.UTCDateTime.strptime(datestr, dateformat), -offset])
sac_files.append([datestr, UTCDateTime.strptime(datestr, dateformat), -offset])
elif offset is None:
try:
tr = obspy.read(ref_sac)[0]
......@@ -100,7 +103,7 @@ def match_eq(eq_lst, pathname, stla, stlo, ref_comp='Z', suffix='SAC', offset=No
sac_files.append([datestr, tr.stats.starttime-tr.stats.sac.b, tr.stats.sac.o])
else:
raise TypeError('offset should be int or float type')
new_col = ['dis', 'bazi', 'data']
new_col = ['dis', 'bazi', 'data', 'datestr']
eq_match = pd.DataFrame(columns=new_col)
for datestr, b_time, offs in sac_files:
date_range_begin = b_time + timedelta(seconds=offs - tolerance)
......@@ -109,14 +112,14 @@ def match_eq(eq_lst, pathname, stla, stlo, ref_comp='Z', suffix='SAC', offset=No
if len(results) != 1:
continue
try:
this_eq = eq(pathname, datestr, suffix, switchEN=switchEN, reverseE=reverseE, reverseN=reverseN)
this_eq = eq(pathname, datestr, suffix)
except Exception as e:
continue
this_eq.get_time_offset(results.iloc[0]['date'])
daz = seispy.distaz(stla, stlo, results.iloc[0]['evla'], results.iloc[0]['evlo'])
this_df = pd.DataFrame([[daz.delta, daz.baz, this_eq]], columns=new_col, index=results.index.values)
daz = distaz(stla, stlo, results.iloc[0]['evla'], results.iloc[0]['evlo'])
this_df = pd.DataFrame([[daz.delta, daz.baz, this_eq, datestr]], columns=new_col, index=results.index.values)
eq_match = eq_match.append(this_df)
ind = eq_match.index.drop_duplicates(False)
ind = eq_match.index.drop_duplicates(keep=False)
eq_match = eq_match.loc[ind]
return pd.concat([eq_lst, eq_match], axis=1, join='inner')
......@@ -153,9 +156,9 @@ def CfgParser(cfg_file):
for sec in sections:
for key, value in cf.items(sec):
if key == 'date_begin':
pa.__dict__[key] = obspy.UTCDateTime(value)
pa.__dict__[key] = UTCDateTime(value)
elif key == 'date_end':
pa.__dict__[key] = obspy.UTCDateTime(value)
pa.__dict__[key] = UTCDateTime(value)
elif key == 'offset':
try:
pa.__dict__[key] = float(value)
......@@ -275,7 +278,7 @@ class RF(object):
self.eqs = match_eq(self.eq_lst, self.para.datapath, self.stainfo.stla, self.stainfo.stlo,
ref_comp=self.para.ref_comp, suffix=self.para.suffix,
offset=self.para.offset, tolerance=self.para.tolerance,
dateformat=self.para.dateformat, switchEN=switchEN, reverseE=reverseE, reverseN=reverseN)
dateformat=self.para.dateformat)
except Exception as e:
self.logger.RFlog.error('{0}'.format(e))
raise e
......@@ -285,37 +288,44 @@ class RF(object):
else:
self.logger.RFlog.info('{0} earthquakes matched'.format(self.eqs.shape[0]))
'''
def save(self, path=''):
if path == '':
path = '{0}.{1}.h5'.format(self.stainfo.network, self.stainfo.station)
d = {'para': self.para.__dict__, 'stainfo': self.stainfo.__dict__, 'eq_lst': self.eq_lst, 'eqs': self.eqs}
path = '{0}.{1}_matched.pkl'.format(self.stainfo.network, self.stainfo.station)
try:
self.logger.RFlog.info('Saving project to {0}'.format(path))
dd.io.save(path, d)
except PerformanceWarning:
pass
save_eqs = self.eqs[['date', 'evla', 'evlo', 'evdp', 'mag', 'bazi', 'dis', 'datestr']]
with open(path, 'wb') as f:
pickle.dump({'para': self.para, 'eqs': save_eqs}, f, -1)
except Exception as e:
self.logger.RFlog.error('{0}'.format(e))
raise IOError(e)
def load(self, path):
try:
self.logger.RFlog.info('Loading {0}'.format(path))
fdd = dd.io.load(path)
except Exception as e:
self.logger.RFlog.error('{0}'.format(e))
raise IOError('Cannot read {0}'.format(path))
try:
self.para.__dict__.update(fdd['para'])
self.stainfo.__dict__.update(fdd['stainfo'])
self.eq_lst = fdd['eq_lst']
self.eqs = fdd['eqs']
except Exception as e:
raise ValueError(e)
'''
with open(path, 'rb') as f:
rfdata = pickle.load(f)
self.para = rfdata['para']
if not exists(self.para.datapath):
self.logger.RFlog.error('Data path {} was not found'.format(self.para.datapath))
sys.exit(1)
self.eqs = rfdata['eqs']
eqs = []
for i, row in self.eqs.iterrows():
try:
this_eq = eq(self.para.datapath, row['datestr'], self.para.suffix)
except Exception as e:
self.logger.RFlog.warning('Cannot read {}, skipping'.format(row['datestr']))
continue
this_eq.get_time_offset(row['date'])
eqs.append(this_eq)
self.eqs.insert(1, 'data', eqs)
def channel_correct(self, switchEN=False, reverseE=False, reverseN=False):
self.logger.RFlog.info('Correct components with switchEN: {}, reverseE: {}, reverseN: {}'.format(switchEN, reverseE, reverseN))
self.para.switchEN = switchEN
self.para.reverseE = reverseE
self.para.reverseN = reverseN
for _, row in self.eqs.iterrows():
row['data'].channel_correct(switchEN, reverseE, reverseN)
def detrend(self):
self.logger.RFlog.info('Detrend all data')
......@@ -367,19 +377,14 @@ class RF(object):
self.logger.RFlog.info('Average {:.1f} deg offset in back-azimuth'.format(baz_shift))
self.eqs['bazi'] = np.mod(self.eqs['bazi'] + baz_shift, 360)
# def _baz_confirm(self, offset, ampt_all):
# y = np.arange(self.eqs.shape[0])
# x = np.linspace(-offset, offset, ampt_all.shape[1])
# fig = _plotampt(x, y, ampt_all)
# if messagebox.askyesno('Back-azimuth correction',
# 'Would you like to keep this correction?'):
# fig.close()
# else:
# self.logger.RFlog.error('Manual interruption.')
# fig.close()
# sys.exit(1)
def rotate(self, method='NE->RT', search_inc=False):
def rotate(self, search_inc=False):
targ_comp = ''.join(sorted(self.para.comp.upper()))
if targ_comp == 'RTZ':
method='NE->RT'
elif targ_comp == 'LQT':
method='ZNE->LQT'
else:
raise ValueError('comp must be in RTZ or LQT.')
self.logger.RFlog.info('Rotate {0} phase to {1}'.format(self.para.phase, method))
drop_idx = []
for i, row in self.eqs.iterrows():
......@@ -413,14 +418,8 @@ class RF(object):
self.logger.RFlog.info('{0} events left after virtual checking'.format(self.eqs.shape[0]))
def deconv(self):
if self.para.phase == 'P':
shift = self.para.time_before
time_after = self.para.time_after
elif self.para.phase == 'S':
shift = self.para.time_after
time_after = self.para.time_before
else:
pass
shift = self.para.time_before
time_after = self.para.time_after
drop_lst = []
count = 0
......@@ -458,7 +457,7 @@ class RF(object):
else:
self.logger.RFlog.info('Save RFs with and criterion of {}'.format(criterion))
for i, row in self.eqs.iterrows():
if row['data'].judge_rf(shift, npts, criterion=criterion, rmsgate=self.para.rmsgate):
if row['data'].judge_rf(shift, npts, phase=self.para.phase, criterion=criterion, rmsgate=self.para.rmsgate):
row['data'].saverf(self.para.rfpath, evtstr=row['date'].strftime('%Y.%j.%H.%M.%S'), phase=self.para.phase, shift=shift,
evla=row['evla'], evlo=row['evlo'], evdp=row['evdp'], baz=row['bazi'],
mag=row['mag'], gcarc=row['dis'], gauss=self.para.gauss, only_r=self.para.only_r)
......@@ -503,7 +502,8 @@ def prf():
pjt = RF(cfg_file=arg.cfg_file)
pjt.load_stainfo()
pjt.search_eq(local=arg.islocal)
pjt.match_eq(switchEN=arg.isswitch, reverseN=reverseN, reverseE=reverseE)
pjt.match_eq()
pjt.channel_correct(switchEN=arg.isswitch, reverseN=reverseN, reverseE=reverseE)
pjt.detrend()
pjt.filter()
pjt.cal_phase()
......@@ -515,13 +515,7 @@ def prf():
else:
pass
pjt.trim()
targ_comp = ''.join(sorted(pjt.para.comp.upper()))
if targ_comp == 'RTZ':
pjt.rotate(method='NE->RT')
elif targ_comp == 'LQT':
pjt.rotate(method='ZNE->LQT')
else:
raise ValueError('comp must be in RTZ or LQT.')
pjt.rotate()
pjt.deconv()
pjt.saverf()
......
......@@ -39,12 +39,15 @@ class RPickFigure(RFFigure):
self.axg = self.axb.twiny()
def read_sac(self, dt=0.1):
if len(glob.glob(join(self.rfpath, '*_R.sac'))) == 0:
tmp_files = glob.glob(join(self.rfpath, '*_Q.sac'))
self.comp = 'Q'
else:
if len(glob.glob(join(self.rfpath, '*_R.sac'))) != 0:
tmp_files = glob.glob(join(self.rfpath, '*_R.sac'))
self.comp = 'R'
self.comp = 'R'
elif len(glob.glob(join(self.rfpath, '*_L.sac'))) != 0:
tmp_files = glob.glob(join(self.rfpath, '*_L.sac'))
self.comp = 'L'
else:
tmp_files = glob.glob(join(self.rfpath, '*_Z.sac'))
self.comp = 'Z'
self.log.RFlog.info('Reading PRFs from {}'.format(self.rfpath))
self.filenames = [basename(sac_file).split('_')[0] for sac_file in sorted(tmp_files)]
self.rrf = obspy.read(join(self.rfpath, '*_{}.sac'.format(self.comp))).sort(['starttime']).resample(1/dt)
......
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册