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

updata makeRF4station

上级 0d9dcf81
文件已添加
此差异已折叠。
#!/usr/bin/env python
def Usage():
print('./makeRF4station -Sstation/slat/slon -Yyear1/month1/day1/year2/month2/day2 -Ccomp [-M] [-T+N+E] [-r] para.cfg')
print(" -S Specified station name, station latitude, station longitude")
print(" Note: station name must be same as folder name")
print(" -Y Set date range when search earthquake from GMCT Catalog")
print(" -C Specified componet formats separated by dot in sac-file name, such as \"BHZ\"")
print('./makeRF4station -Sstation -Yyear1/month1/day1/year2/month2/day2 -Ccomp [-M] [-T+N+E] [-r] para.cfg')
import numpy as np
import obspy
......@@ -44,7 +40,6 @@ change_E = 0
change_N = 0
ismtz = 0
ist = 1
isdat = False
for op, value in opts:
if op == "-S":
staname = value
......@@ -222,9 +217,9 @@ for thiseq in eq:
this_seis[1].stats.channel = "T"
this_seis[2].stats.channel = "Z"
# ------- Bandpass Filter ------------#
R_filter = bandpass(this_seis[0].data, freqmin, freqmax, 1/dt, corners=4, zerophase=True)
T_filter = bandpass(this_seis[1].data, freqmin, freqmax, 1/dt, corners=4, zerophase=True)
Z_filter = bandpass(this_seis[2].data, freqmin, freqmax, 1/dt, corners=4, zerophase=True)
R_filter = bandpass(this_seis[0].data, freqmin, freqmax, 1/dt, corners=3, zerophase=True)
T_filter = bandpass(this_seis[1].data, freqmin, freqmax, 1/dt, corners=3, zerophase=True)
Z_filter = bandpass(this_seis[2].data, freqmin, freqmax, 1/dt, corners=3, zerophase=True)
# ------- Calculate P arrival time ---#
arrivals = model.get_travel_times(source_depth_in_km=dep, distance_in_degree=dis, phase_list=["P"])
ttime_P = arrivals[0].time - O
......@@ -259,10 +254,12 @@ for thiseq in eq:
cti = RMS[-1] < 0.2 and max_deep < max_P*0.3 and max_P == np.max(np.abs(RF)) and max_P < 1
if cti:
print("----"+staname+"-----"+date_name+"-----")
tRF, tRMS, tit = seispy.decov.decovit(T, Z, dt, RFlength, time_before, gauss, 400, 0.001)
if ist:
tRF, tRMS, tit = seispy.decov.decovit(T, Z, dt, RFlength, time_before, gauss, 400, 0.001)
if RFlength != newlength:
RF = resample(RF, newlength)
tRF = resample(tRF, newlength)
if ist:
tRF = resample(tRF, newlength)
this_seis[0].data = this_seis[0].data[extbegin:extend+1]
this_seis[1].data = this_seis[1].data[extbegin:extend+1]
this_seis[2].data = this_seis[2].data[extbegin:extend+1]
......@@ -281,14 +278,14 @@ for thiseq in eq:
this_seis[i].stats.sac.mag = mag
if not os.path.exists(out_path):
os.makedirs(out_path)
this_seis[0].write(os.path.join(out_path, date_name+'.'+staname+'.R.SAC'), 'SAC')
this_seis[1].write(os.path.join(out_path, date_name+'.'+staname+'.T.SAC'), 'SAC')
this_seis[2].write(os.path.join(out_path, date_name+'.'+staname+'.Z.SAC'), 'SAC')
RF_R = this_seis[0].copy()
RF_R.stats.sac.user1 = gauss
RF_R.stats.sac.b = -time_before
RF_R.stats.sac.a = 0
RF_R.data = RF
this_seis[0].write(os.path.join(out_path, date_name+'.'+staname+'.R.SAC'), 'SAC')
this_seis[1].write(os.path.join(out_path, date_name+'.'+staname+'.T.SAC'), 'SAC')
this_seis[2].write(os.path.join(out_path, date_name+'.'+staname+'.Z.SAC'), 'SAC')
RF_R.write(os.path.join(RF_path, date_name+'_P_R.sac'), 'SAC')
if ist:
RF_T = this_seis[1].copy()
......
文件模式从 100644 更改为 100755
......@@ -27,18 +27,24 @@ def Usage():
print("Usage: python pickrf_R.py -Sstation_name para.cfg")
def get_sac():
if sys.argv[1:] == []:
Usage()
sys.exit(1)
for o in sys.argv[1:]:
if os.path.isfile(o):
head = o
break
try:
opts, args = getopt.getopt(sys.argv[1:], "S:")
opts, args = getopt.getopt(sys.argv[1:], "S:h")
except:
print("invalid argument")
sys.exit(1)
for op, value in opts:
if op == "-S":
staname = value
elif op == "-h":
Usage()
sys.exit(1)
else:
print("invalid argument")
sys.exit(1)
......@@ -98,8 +104,8 @@ class plotrffig():
self.plotwave()
self.plotbaz()
self.fig.suptitle("%s (Latitude: %5.2f\N{DEGREE SIGN}, Longitude: %5.2f\N{DEGREE SIGN})" % (opts.staname, opts.stla, opts.stlo), fontsize=20)
ax.set_ylim(rfidx[self.ipage][0], rfidx[self.ipage][-1]+1)
ax.set_yticks(np.arange(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+1))
ax.set_ylim(rfidx[self.ipage][0], rfidx[self.ipage][-1]+2)
ax.set_yticks(np.arange(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+2))
ylabels = opts.filenames[rfidx[self.ipage][0]::]
ylabels.insert(0, '')
ax.set_yticklabels(ylabels)
......@@ -203,12 +209,12 @@ class plotrffig():
if self.ipage < 0:
self.ipage = 0
return
ax.set_ylim(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+1)
ax.set_yticks(np.arange(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+1))
ax.set_ylim(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+2)
ax.set_yticks(np.arange(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+2))
ylabels = opts.filenames[self.rfidx[self.ipage][0]::]
ylabels.insert(0, '')
ax.set_yticklabels(ylabels)
ax_baz.set_ylim(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+1)
ax_baz.set_ylim(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+2)
ax_baz.set_yticks(ax.get_yticks())
self.azi_label = ['%5.2f' % opts.baz[i] for i in self.rfidx[self.ipage]]
self.azi_label.insert(0, "")
......@@ -233,7 +239,7 @@ class plotrffig():
self.azi_label = ['%5.2f' % opts.baz[i] for i in self.rfidx[self.ipage]]
else:
ax.set_ylim(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+2)
ax.set_yticks(np.arange(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+1))
ax.set_yticks(np.arange(self.rfidx[self.ipage][0], self.rfidx[self.ipage][-1]+2))
ylabels = opts.filenames[self.rfidx[self.ipage][0]::]
ylabels.insert(0, '')
ax.set_yticklabels(ylabels)
......@@ -247,8 +253,8 @@ class plotrffig():
def main():
opts = initopts.opts()
opts.maxidx = 20
opts.enf = 5
opts.xlim = [-2, 30]
opts.enf = 7
opts.xlim = [-2, 80]
opts.ylim = [0, 22]
opts.rffiles, opts.filenames, opts.path, opts.image_path, opts.cut_path, opts.staname = get_sac()
opts.evt_num = len(opts.rffiles)
......
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册