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

add back-azimuth correction

上级 bc4de82c
......@@ -51,6 +51,7 @@ class eq(object):
def detrend(self):
self.st.detrend(type='linear')
self.st.detrend(type='constant')
self.fix_channel_name()
def filter(self, freqmin=0.05, freqmax=1, order=4):
self.st.filter('bandpass', freqmin=freqmin, freqmax=freqmax, corners=order, zerophase=True)
......@@ -84,7 +85,41 @@ class eq(object):
real_inc = inc_range[real_inc_idx]
return real_inc
def rotate(self, bazi, inc=None, method='NE->RT', phase='P', inv=None):
def search_baz(self, bazi, time_b=10, time_e=20, offset=90):
p_arr, _ = self.arr_correct(write_to_sac=False)
this_st = self.st.copy()
this_st.filter('bandpass', freqmin=0.03, freqmax=0.5)
this_st.trim(this_st[0].stats.starttime+p_arr-time_b, this_st[0].stats.starttime+p_arr+time_e)
bazs = np.arange(bazi-offset, bazi+offset)
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)
idx = seispy.geo.extrema(ampt, opt='min')
if len(idx) == 0:
return np.nan
elif len(idx) > 1:
return None
else:
return bazs[seispy.geo.extrema(ampt, opt='min')[0]] - bazi
def fix_channel_name(self):
if self.st.select(channel='??1') and hasattr(self.st.select(channel='*1')[0].stats.sac, 'cmpaz'):
if self.st.select(channel='*1')[0].stats.sac.cmpaz == 0:
self.st.select(channel='*1')[0].stats.channel = self.st.select(channel='*1')[0].stats.channel[:-1] + 'N'
self.st.select(channel='*2')[0].stats.channel = self.st.select(channel='*2')[0].stats.channel[:-1] + 'E'
elif self.st.select(channel='*1')[0].stats.sac.cmpaz != 0:
rotateZNE(self.st)
self.st.sort(['channel'])
elif self.st.select(channel='1'):
self.st.select(channel='1')[0].stats.channel = 'Z'
self.st.select(channel='2')[0].stats.channel = 'N'
self.st.select(channel='3')[0].stats.channel = 'E'
self.st.sort(['channel'])
else:
pass
def rotate(self, bazi, inc=None, method='NE->RT', phase='P'):
if phase not in ('P', 'S'):
raise ValueError('phase must be in [\'P\', \'S\']')
......@@ -102,23 +137,6 @@ class eq(object):
else:
pass
if inv is not None:
self.rf.rotate('->ZNE', inventory=inv)
elif self.rf.select(channel='??1') and hasattr(self.rf.select(channel='*1')[0].stats.sac, 'cmpaz'):
if self.rf.select(channel='*1')[0].stats.sac.cmpaz == 0:
self.rf.select(channel='*1')[0].stats.channel = self.rf.select(channel='*1')[0].stats.channel[:-1] + 'N'
self.rf.select(channel='*2')[0].stats.channel = self.rf.select(channel='*2')[0].stats.channel[:-1] + 'E'
elif self.rf.select(channel='*1')[0].stats.sac.cmpaz != 0:
rotateZNE(self.rf)
self.rf.sort(['channel'])
elif self.rf.select(channel='1'):
self.rf.select(channel='1')[0].stats.channel = 'Z'
self.rf.select(channel='2')[0].stats.channel = 'N'
self.rf.select(channel='3')[0].stats.channel = 'E'
self.rf.sort(['channel'])
else:
pass
if method == 'NE->RT':
self.rf.rotate('NE->RT', back_azimuth=bazi)
elif method == 'RT->NE':
......@@ -237,16 +255,18 @@ class eq(object):
tr.o = 0
tr.write(filename + '_{0}_{1}.sac'.format(phase, tr.kcmpnm[-1]))
def judge_rf(self, shift, criterion='crust'):
def judge_rf(self, shift, npts, criterion='crust'):
if not isinstance(criterion, (str, type(None))):
raise TypeError('criterion should be string in [\'crust\', \'mtz\']')
elif criterion is None:
pass
elif criterion.lower() not in ['crust', 'mtz', 'lab']:
elif criterion.lower() not in ['crust', 'mtz']:
raise ValueError('criterion should be string in [\'crust\', \'mtz\']')
else:
criterion = criterion.lower()
if self.rf[1].stats.npts != npts:
return False
if criterion == 'crust':
time_P1 = int(np.floor((-2+shift)/self.rf[1].stats.delta))
time_P2 = int(np.floor((2+shift)/self.rf[1].stats.delta))
......
......@@ -97,3 +97,14 @@ class para(object):
if not isinstance(value, str):
raise TypeError('catalogpath should be \'str\' type not \'{0}\''.format(type(value)))
self._catalogpath = value
@property
def criterion(self):
return self._criterion
@criterion.setter
def criterion(self, value):
if value == '':
self._criterion = None
else:
self._criterion = value
......@@ -10,11 +10,13 @@ import seispy
from seispy.eq import eq
from seispy.setuplog import setuplog
import glob
import numpy as np
from datetime import timedelta, datetime
import pandas as pd
from obspy.taup import TauPyModel
import configparser
import argparse
import sys
def datestr2regex(datestr):
......@@ -294,6 +296,19 @@ class RF(object):
row['data'].get_arrival(self.model, row['evdp'], row['dis'])
# row['data'].get_raypara(self.model, row['evdp'], row['dis'])
def baz_correct(self, time_b=10, time_e=20, offset=90):
self.logger.RFlog.info('correct back-azimuth')
shift_all = []
for i, row in self.eqs.iterrows():
shift_all.append(row['data'].search_baz(row['bazi'], time_b=time_b, time_e=time_e, offset=offset))
if None in shift_all:
self.logger.RFlog.error('Range of searching bazi is too small.')
sys.exit(1)
shift_all = np.array(shift_all)
baz_shift = np.mean(shift_all[np.where(np.logical_not(np.isnan(shift_all)))])
self.logger.RFlog.info('Average {:.4f} deg offset in back-azimuth'.format(baz_shift))
self.eqs['bazi'] = self.eqs['bazi'] + baz_shift
def rotate(self, method='NE->RT'):
self.logger.RFlog.info('Rotate {0} phase {1}'.format(self.para.phase, method))
drop_idx = []
......@@ -350,6 +365,7 @@ class RF(object):
def saverf(self):
if self.para.phase == 'P':
shift = self.para.time_before
npts = int((self.para.time_before + self.para.time_after)/self.para.target_dt+1)
criterion = self.para.criterion
elif self.para.phase == 'S':
shift = self.para.time_after
......@@ -360,7 +376,7 @@ class RF(object):
self.logger.RFlog.info('Save RFs with criterion of {}'.format(criterion))
for i, row in self.eqs.iterrows():
if row['data'].judge_rf(shift, criterion=criterion):
if row['data'].judge_rf(shift, npts, criterion=criterion):
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)
......@@ -373,15 +389,40 @@ def prf():
parser = argparse.ArgumentParser(description="Calculating RFs for single station")
parser.add_argument('cfg_file', type=str, help='Path to RF configure file')
parser.add_argument('-l', help="use local catalog. Default is false", dest='islocal', action='store_true')
parser.add_argument('-r', help='Reverse components: EN, E or N', default=None, type=str)
parser.add_argument('-b', help='Correct back-azimuth via minimal energy of T component', default=None, type=str)
arg = parser.parse_args()
if arg.r is not None:
arg.r = arg.r.upper()
if arg.r == 'NE' or arg.r == 'EN':
reverseE = True
reverseN = True
elif arg.r == 'E':
reverseE = True
reverseN = False
elif arg.r == 'N':
reverseE = False
reverseN = True
else:
raise ValueError('component name must be in EN, E or N')
else:
reverseN = False
reverseE = False
if arg.b is not None:
try:
baz_corr = [float(val) for val in arg.b.split('/')]
except:
raise ValueError('Format of baz correction must be as 10/20/45')
pjt = RF(cfg_file=arg.cfg_file)
pjt.load_stainfo()
pjt.search_eq(local=arg.islocal)
pjt.match_eq()
pjt.match_eq(reverseN=reverseN, reverseE=reverseE)
pjt.detrend()
pjt.filter()
pjt.cal_phase()
pjt.drop_eq_snr()
if arg.b is not None:
pjt.baz_correct(baz_corr[0], baz_corr[1], baz_corr[2])
pjt.trim()
pjt.rotate()
pjt.deconv()
......
......@@ -19,7 +19,7 @@ def convertinfo(info):
def ndkparse(ndk_str):
find_re = re.compile(
r'[A-Z]+\s(\d+)/(\d+)/(\d+)\s+(\d+):(\d+):(\d+.\d)\s+(.+?)\s+(.+?)\s+(.+?)\s+.+?\s+(.+?)\s+.+?\\n', re.DOTALL)
r'[A-Z]+\s(\d+)/(\d+)/(\d+)\s+(\d+):(\d+):(\d+.\d)\s+(.+?)\s+(.+?)\s+(.+?)\s+.+?\s+(.+?)\s+.+?\n', re.DOTALL)
ndk_lst = [convertinfo(info) for info in find_re.findall(ndk_str)]
return ndk_lst
......@@ -59,7 +59,9 @@ def main():
parser = argparse.ArgumentParser(description="Update CMT Catalog")
parser.add_argument('-i', help='Input Catalog', dest='inlog',
type=str, default=join(dirname(__file__), 'data', 'EventCMT.dat'))
parser.add_argument('-o', help='Onput Catalog', dest='outlog', type=str, default='')
parser.add_argument('-o', help='Onput Catalog', dest='outlog', type=str,
default=join(dirname(__file__), 'data', 'EventCMT.dat'))
parser.add_argument('-u', help='url of ndk', type=str)
arg = parser.parse_args()
fetch_cata(inlog=arg.inlog, outlog=arg.outlog)
......@@ -77,7 +79,10 @@ def ndk2dat():
except Exception as e:
raise TypeError('Error type for ndk file\n{}'.format(e))
with open(arg.outlog, 'w+') as f:
f.writelines(log.lst)
for year, mon, day, hour, min, sec, lat, lon, dep, mw in log_lst:
evt_time = datetime(year, mon, day, hour, min, int(sec))
f.write('%d %d %d %s %d %d %d %6.2f %6.2f %s %s\n' % (
year, mon, day, evt_time.strftime('%j'), hour, min, int(sec), lat, lon, dep, mw))
if __name__ == '__main__':
......
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册