makeRF4station.py 12.1 KB
Newer Older
Mijian Xu's avatar
Mijian Xu committed
1
2
3
#!/usr/bin/env python

def Usage():
Mijian Xu's avatar
Mijian Xu committed
4
   print('./makeRF4station -Sstation/slat/slon -Yyear1/month1/day1/year2/month2/day2 -Ccomp [-M] [-T+N+E] [-r] para.cfg')
xumj's avatar
xumj committed
5
6
7
8
   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\"")
Mijian Xu's avatar
Mijian Xu committed
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

import numpy as np
import obspy
import re
import os
import glob
import getopt
import sys
import datetime
import seispy
from obspy.taup import TauPyModel
from scipy.signal import detrend
from obspy.signal.filter import bandpass
from scipy.signal import resample
try:
    import configparser
    config = configparser.ConfigParser()
except:
    import ConfigParser
    config = ConfigParser.ConfigParser()

#################################
# ------- Get arguments --------#
#################################
try:
Mijian Xu's avatar
Mijian Xu committed
34
   opts, args = getopt.getopt(sys.argv[1:], "S:Y:C:T:Mr")
Mijian Xu's avatar
Mijian Xu committed
35
36
37
38
39
40
41
42
43
44
except:
    print('Arguments are not found!')
    Usage()
    sys.exit(1)
if opts == []:
    Usage()
    sys.exit(1)

change_E = 0
change_N = 0
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
45
46
ismtz = 0
ist = 1
xumj's avatar
xumj committed
47
isdat = False
Mijian Xu's avatar
Mijian Xu committed
48
49
for op, value in opts:
    if op == "-S":
Mijian Xu's avatar
Mijian Xu committed
50
51
       staname = value
#       station = value
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
52
53
    elif op == "-M":
       ismtz = 1
Mijian Xu's avatar
Mijian Xu committed
54
55
56
57
58
59
60
61
    elif op == "-Y":
       date_str = value
    elif op == "-C":
        comp = value
    elif op == "-T":
        change_chan = value
        change_E = change_chan.find("E")
        change_N = change_chan.find("N")
Mijian Xu's avatar
Mijian Xu committed
62
    elif op == "-r":
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
63
        ist = 0
Mijian Xu's avatar
Mijian Xu committed
64
65
66
67
68
69
70
71
72
73
74
75
76
77
    else:
       Usage()
       sys.exit(1)

for o in sys.argv[1:]:
    if os.path.isfile(o):
        head = o
        config.read(head)
        break

#################################
# --------- Get Para -----------#
#################################
model = TauPyModel(model="iasp91")
Mijian Xu's avatar
Mijian Xu committed
78
79
80
#staname = station.split('/')[0]
#stalat = float(station.split('/')[1])
#stalon = float(station.split('/')[2])
Mijian Xu's avatar
Mijian Xu committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
yearrange1 = int(date_str.split('/')[0])
monthrange1 = int(date_str.split('/')[1])
dayrange1 = int(date_str.split('/')[2])
yearrange2 = int(date_str.split('/')[3])
monthrange2 = int(date_str.split('/')[4])
dayrange2 = int(date_str.split('/')[5])
daterange1 = datetime.datetime(yearrange1, monthrange1, dayrange1)
daterange2 = datetime.datetime(yearrange2, monthrange2, dayrange2)
data_path = config.get('path', 'data_path')
out_path = config.get('path', 'out_path')
out_path = os.path.join(out_path, staname)
RF_path = config.get('path', 'RF_path')
RF_path = os.path.join(RF_path, staname)
evt_list = config.get('path', 'evt_list')
image_path = config.get('path', 'image_path')
if not os.path.exists(image_path):
    os.makedirs(image_path)
gate_mw = float(config.get('para', 'gate_mw'))
gate_dis1 = float(config.get('para', 'gate_dis1'))
gate_dis2 = float(config.get('para', 'gate_dis2'))
time_before = float(config.get('para', 'time_before'))
time_after = float(config.get('para', 'time_after'))
Mijian Xu's avatar
Mijian Xu committed
103
tolerance = float(config.get('para', 'tolerance'))
Mijian Xu's avatar
Mijian Xu committed
104
105
106
107
108
109
110
111
112
113
114
115
116
117
offset = float(config.get('para', 'offset'))
gate_noise = float(config.get('para', 'gate_noise'))
gauss = float(config.get('para', 'gauss'))
freqmin = float(config.get('para', 'freqmin'))
freqmax = float(config.get('para', 'freqmax'))
sampling = float(config.get('para', 'sampling'))
newlength = ((time_before+time_after)/sampling)+1
image_path = os.path.join(image_path, staname+'_R_gauss'+str(gauss)+'.ps')
fid_evtlist = open(evt_list, 'r')

#################################
# --------- Search eq ----------#
#################################
eq_lst = []
Mijian Xu's avatar
Mijian Xu committed
118
ex_sac = obspy.read(glob.glob(os.path.join(data_path, staname, '*.'+comp+'.*[Ss][Aa][Cc]'))[0])[0]
Mijian Xu's avatar
Mijian Xu committed
119
120
stalat = ex_sac.stats.sac.stla
stalon = ex_sac.stats.sac.stlo
Mijian Xu's avatar
Mijian Xu committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
for evt in fid_evtlist.readlines():
    evt = evt.strip('\n')
    year_evt = int(evt.split()[0])
    mon_evt = int(evt.split()[1])
    day_evt = int(evt.split()[2])
    hour_evt = int(evt.split()[4])
    min_evt = int(evt.split()[5])
    sec_evt = int(evt.split()[6])
    lat_evt = float(evt.split()[7])
    lon_evt = float(evt.split()[8])
    dis_evt = seispy.distaz(stalat, stalon, lat_evt, lon_evt).delta
    bazi_evt = seispy.distaz(stalat, stalon, lat_evt, lon_evt).getBaz()
    dep_evt = float(evt.split()[9])
    mw_evt = float(evt.split()[10])
    date_evt = datetime.datetime(year_evt, mon_evt, day_evt, hour_evt, min_evt, sec_evt)
    if daterange1 <= date_evt <= daterange2 and gate_dis1 <= dis_evt <= gate_dis2 and mw_evt > gate_mw:
        matchedeq = [date_evt, lat_evt, lon_evt, dep_evt, dis_evt, bazi_evt, mw_evt]
        eq_lst.append(matchedeq)

#################################
# -------- Assign Files --------#
#################################
eq = []
for nowevt in eq_lst:
    find_result = []
Mijian Xu's avatar
Mijian Xu committed
146
    for sac in glob.glob(os.path.join(data_path, staname, '*.'+comp+'.*[Ss][Aa][Cc]')):
Mijian Xu's avatar
Mijian Xu committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
        sacname = os.path.basename(sac)
        date_name = re.search("\d{4}\D\d{3}\D\d{2}\D\d{2}\D\d{2}", sacname).group()
        date_name_sp = re.split('\D', date_name)
        year_sac = int(date_name_sp[0])
        jday_sac = int(date_name_sp[1])
        foo = datetime.datetime(year_sac, 1, 1) + datetime.timedelta(jday_sac - 1)
        mon_sac = foo.month
        day_sac = foo.day
        hour_sac = int(date_name_sp[2])
        min_sac = int(date_name_sp[3])
        sec_sac = int(date_name_sp[4])
        if sec_sac == 60:
            sec_sac = 0
            date_sac = datetime.datetime(year_sac, mon_sac, day_sac, hour_sac, min_sac, sec_sac) + datetime.timedelta(seconds=60)
        else:
            date_sac = datetime.datetime(year_sac, mon_sac, day_sac, hour_sac, min_sac, sec_sac)
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
163
        
Mijian Xu's avatar
Mijian Xu committed
164
165
        # if nowevt[0] + datetime.timedelta(seconds=offset) - datetime.timedelta(seconds=tolerance) <= date_sac <= nowevt[0] + datetime.timedelta(seconds=offset) + datetime.timedelta(seconds=tolerance):
        if date_sac - datetime.timedelta(seconds=offset) - datetime.timedelta(seconds=tolerance) <= nowevt[0] <= date_sac - datetime.timedelta(seconds=offset) + datetime.timedelta(seconds=tolerance):
Mijian Xu's avatar
Mijian Xu committed
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
            print(date_sac, nowevt[0])
            ocr = (date_sac - nowevt[0])
            ocr = ocr.total_seconds()
            find_result.append([nowevt[0], date_name, nowevt[1], nowevt[2], nowevt[3], nowevt[4], nowevt[5], nowevt[6], ocr])
    if len(find_result) == 1:
        eq.append(find_result[0])
print(len(eq))

#################################
# -------- Process RFs ---------#
#################################
if not os.path.exists(RF_path):
    os.makedirs(RF_path)
for thiseq in eq:
    date_name = thiseq[1]
    bazi = thiseq[6]
    dis = thiseq[5]
    dep = thiseq[4]
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
184
    mag = thiseq[7]
Mijian Xu's avatar
Mijian Xu committed
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
    O = thiseq[-1]
    try:
        this_seis = obspy.core.read(os.path.join(data_path, staname, '*'+date_name+'*.[Ss][Aa][Cc]'))
    except:
        continue
# --------- Calculate dt -----------#
    try:
        dt_E = np.mean(np.diff(this_seis[0].times()))
        dt_N = np.mean(np.diff(this_seis[1].times()))
        dt_Z = np.mean(np.diff(this_seis[2].times()))
    except:
        continue
    if not (dt_E == dt_N == dt_Z and len(this_seis[0].data) == len(this_seis[1].data) == len(this_seis[2].data)):
        continue
    dt = dt_E
# -------- Detrend -----------------#
    if change_E:
        this_seis[0].data = this_seis[0].data * -1
    if change_N:
        this_seis[1].data = this_seis[1].data * -1
    this_seis[0].data = detrend(this_seis[0].data, type='linear')
    this_seis[0].data = detrend(this_seis[0].data, type='constant')
    this_seis[1].data = detrend(this_seis[1].data, type='linear')
    this_seis[1].data = detrend(this_seis[1].data, type='constant')
    this_seis[2].data = detrend(this_seis[2].data, type='linear')
    this_seis[2].data = detrend(this_seis[2].data, type='constant')
# ------- Rotate ENZ to RTZ ---------#
    if comp == 1 or comp == 2 or comp == 3:
        E = this_seis[2].data
        N = this_seis[1].data
        Z = this_seis[0].data
    else:
        E = this_seis[0].data
        N = this_seis[1].data
        Z = this_seis[2].data
    (this_seis[0].data, this_seis[1].data, this_seis[2].data) = seispy.geo.rotateSeisENZtoTRZ(E, N, Z, bazi)
    this_seis[0].stats.channel = "R"
    this_seis[1].stats.channel = "T"
    this_seis[2].stats.channel = "Z"
# ------- Bandpass Filter ------------#
xumj's avatar
xumj committed
225
226
227
    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)
Mijian Xu's avatar
Mijian Xu committed
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
# ------- 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
    rayp = seispy.geo.srad2skm(arrivals[0].ray_param)
    # this_seis[0].stats.sac.t1 = ttime_P
    # this_seis[1].stats.sac.t1 = ttime_P
    # this_seis[2].stats.sac.t1 = ttime_P
# ----------- filter by snr ----------#
    snrbegin = np.floor((ttime_P-50)/dt)
    snrend = np.floor((ttime_P+50)/dt)
    snro = np.floor(ttime_P/dt)
    snr_R = seispy.geo.snr(R_filter[snro:snrend], R_filter[snrbegin:snro])
    snr_T = seispy.geo.snr(T_filter[snro:snrend], T_filter[snrbegin:snro])
    snr_Z = seispy.geo.snr(Z_filter[snro:snrend], Z_filter[snrbegin:snro])
    print(date_name, snr_R, snr_Z, ttime_P, O)
    if snr_R > gate_noise and snr_Z > gate_noise:
        extbegin = np.floor((ttime_P-time_before)/dt)
        extend = np.floor((ttime_P+time_after)/dt)
        R = R_filter[extbegin:extend+1]
        T = T_filter[extbegin:extend+1]
        Z = Z_filter[extbegin:extend+1]
        RFlength = R.shape[0]
# --- calculate receiver functions ----#
        (RF, RMS, it) = seispy.decov.decovit(R, Z, dt, RFlength, time_before, gauss, 400, 0.001)
        max_deep = np.max(np.abs(RF[np.floor((25+time_before)/dt):-1]))
        time_P1 = np.floor((-2+time_before)/dt)
        time_P2 = np.floor((2+time_before)/dt)
        max_P = np.max(RF[time_P1:time_P2])
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
256
257
258
259
260
        if ismtz == 0:
           cti = max_P == np.max(np.abs(RF)) and max_P < 1
        else:
           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:
Mijian Xu's avatar
Mijian Xu committed
261
            print("----"+staname+"-----"+date_name+"-----")
Mijian Xu's avatar
Mijian Xu committed
262
            tRF, tRMS, tit = seispy.decov.decovit(T, Z, dt, RFlength, time_before, gauss, 400, 0.001)
Mijian Xu's avatar
Mijian Xu committed
263
264
            if RFlength != newlength:
                RF = resample(RF, newlength)
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
265
                tRF = resample(tRF, newlength)
Mijian Xu's avatar
Mijian Xu committed
266
267
268
269
270
271
272
273
274
275
276
            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]
            timeaxis = np.linspace(-time_before, time_after, RFlength)
            for i in range(3):
                this_seis[i].stats.sac.evla = thiseq[2]
                this_seis[i].stats.sac.evlo = thiseq[3]
                this_seis[i].stats.sac.evdp = dep
                this_seis[i].stats.sac.baz = bazi
                this_seis[i].stats.sac.gcarc = dis
                this_seis[i].stats.sac.user0 = rayp
Mijian Xu's avatar
Mijian Xu committed
277
278
279
                this_seis[i].stats.sac.user2 = freqmin
                this_seis[i].stats.sac.user3 = freqmax
                this_seis[i].stats.sac.a = time_before
Mijian Xu's avatar
Mijian Xu committed
280
                this_seis[i].stats.sac.ka = 'P'
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
281
                this_seis[i].stats.sac.mag = mag
Mijian Xu's avatar
Mijian Xu committed
282
283
            if not os.path.exists(out_path):
                os.makedirs(out_path)
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
284
285
286
            RF_R = this_seis[0].copy()
            RF_R.stats.sac.user1 = gauss
            RF_R.stats.sac.b = -time_before
Mijian Xu's avatar
Mijian Xu committed
287
            RF_R.stats.sac.a = 0
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
288
            RF_R.data = RF
xumj's avatar
xumj committed
289
290
291
            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')
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
292
293
            RF_R.write(os.path.join(RF_path,  date_name+'_P_R.sac'), 'SAC')
            if ist:
Mijian Xu's avatar
Mijian Xu committed
294
                RF_T = this_seis[1].copy()
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
295
296
                RF_T.stats.sac.user1 = gauss
                RF_T.stats.sac.b = -time_before
Mijian Xu's avatar
Mijian Xu committed
297
                RF_T.stats.sac.a = 0
Mijian Xu's avatar
Add TRF    
Mijian Xu committed
298
299
                RF_T.data = tRF
                RF_T.write(os.path.join(RF_path,  date_name+'_P_T.sac'), 'SAC')