marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

Included AICPicker of class Picker

Browse Source
This commit is contained in:
Ludger Küperkoch 2014-12-05 16:32:55 +01:00
parent ca809c4673
commit 479058a41e
1 changed files with 173 additions and 155 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

328
pylot/core/pick/run_makeCF.py
View File

@ -10,18 +10,20 @@ from obspy.core import read
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from CharFuns import * from CharFuns import *
from Picker import *
import glob import glob
import argparse import argparse
import pdb
def run_makeCF(project, database, event, iplot, station=None): def run_makeCF(project, database, event, iplot, station=None):
#parameters for CF calculation #parameters for CF calculation
t2 = 7 #length of moving window for HOS calculation [sec] t2 = 7 #length of moving window for HOS calculation [sec]
p = 4 #order of statistics p = 4 #order of statistics
cuttimes = [10, 40] #start and end time vor CF calculation cuttimes = [10, 40] #start and end time vor CF calculation
bpz = [2, 30] #corner frequencies of bandpass filter, vertical component bpz = [2, 30] #corner frequencies of bandpass filter, vertical component
bph = [2, 15] #corner frequencies of bandpass filter, horizontal components bph = [2, 15] #corner frequencies of bandpass filter, horizontal components
tdetz= 1.2 #length of AR-determination window [sec], vertical component tdetz= 1.2 #length of AR-determination window [sec], vertical component
tdeth= 0.8 #length of AR-determination window [sec], horizontal components tdeth= 0.8 #length of AR-determination window [sec], horizontal components
tpredz = 0.4 #length of AR-prediction window [sec], vertical component tpredz = 0.4 #length of AR-prediction window [sec], vertical component
tpredh = 0.4 #length of AR-prediction window [sec], horizontal components tpredh = 0.4 #length of AR-prediction window [sec], horizontal components
addnoise = 0.001 #add noise to seismogram for stable AR prediction addnoise = 0.001 #add noise to seismogram for stable AR prediction
@ -29,165 +31,181 @@ def run_makeCF(project, database, event, iplot, station=None):
arhorder = 4 #chosen order of AR process, horizontal components arhorder = 4 #chosen order of AR process, horizontal components
#get waveform data #get waveform data
if station: if station:
dpz = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*EHZ.msd' % (project, database, event, station) dpz = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*EHZ.msd' % (project, database, event, station)
dpe = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*EHE.msd' % (project, database, event, station) dpe = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*EHE.msd' % (project, database, event, station)
dpn = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*EHN.msd' % (project, database, event, station) dpn = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*EHN.msd' % (project, database, event, station)
#dpz = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*_z.gse' % (project, database, event, station)
#dpe = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*_e.gse' % (project, database, event, station)
#dpn = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*_n.gse' % (project, database, event, station)
else: else:
dpz = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/*EHZ.msd' % (project, database, event) dpz = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/*EHZ.msd' % (project, database, event)
dpe = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/*EHE.msd' % (project, database, event) dpe = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/*EHE.msd' % (project, database, event)
dpn = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/*EHN.msd' % (project, database, event) dpn = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/*EHN.msd' % (project, database, event)
wfzfiles = glob.glob(dpz) wfzfiles = glob.glob(dpz)
wfefiles = glob.glob(dpe) wfefiles = glob.glob(dpe)
wfnfiles = glob.glob(dpn) wfnfiles = glob.glob(dpn)
if wfzfiles: if wfzfiles:
for i in range(len(wfzfiles)): for i in range(len(wfzfiles)):
print 'Vertical component data found ...' print 'Vertical component data found ...'
print wfzfiles[i] print wfzfiles[i]
st = read('%s' % wfzfiles[i]) st = read('%s' % wfzfiles[i])
st_copy = st.copy() st_copy = st.copy()
#filter and taper data #filter and taper data
tr_filt = st[0].copy() tr_filt = st[0].copy()
tr_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False) tr_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False)
tr_filt.taper(max_percentage=0.05, type='hann') tr_filt.taper(max_percentage=0.05, type='hann')
st_copy[0].data = tr_filt.data st_copy[0].data = tr_filt.data
############################################################## ##############################################################
#calculate HOS-CF using subclass HOScf of class CharacteristicFunction #calculate HOS-CF using subclass HOScf of class CharacteristicFunction
hoscf = HOScf(st_copy, cuttimes, t2, p) #instance of HOScf hoscf = HOScf(st_copy, cuttimes, t2, p) #instance of HOScf
############################################################## #get corresponding time array
#calculate AIC-HOS-CF using subclass AICcf of class CharacteristicFunction thoscf = np.arange(0, len(hoscf.getCF()) / tr_filt.stats.sampling_rate, tr_filt.stats.delta) + cuttimes[0]
#class needs stream object => build it ##############################################################
tr_aic = tr_filt.copy() #get onset time from HOS-CF using class Picker
tr_aic.data = hoscf.getCF() #hospick = PragPicker(hoscf.getCF(), thoscf, 2, 70, [1, 0.5, 0.2], 2, 0.001, 0.2)
st_copy[0].data = tr_aic.data #pdb.set_trace()
aiccf = AICcf(st_copy, cuttimes, t2, p) #instance of AICcf ##############################################################
############################################################## #calculate AIC-HOS-CF using subclass AICcf of class CharacteristicFunction
#calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction #class needs stream object => build it
#get stream object of filtered data tr_aic = tr_filt.copy()
st_copy[0].data = tr_filt.data tr_aic.data = hoscf.getCF()
arzcf = ARZcf(st_copy, cuttimes, tpredz, arzorder, tdetz, addnoise) #instance of ARZcf st_copy[0].data = tr_aic.data
############################################################## aiccf = AICcf(st_copy, cuttimes, t2, p) #instance of AICcf
#calculate AIC-ARZ-CF using subclass AICcf of class CharacteristicFunction #get corresponding time array
#class needs stream object => build it taiccf = np.arange(0, len(aiccf.getCF()) / tr_filt.stats.sampling_rate, tr_filt.stats.delta) + cuttimes[0]
tr_arzaic = tr_filt.copy() ##############################################################
tr_arzaic.data = arzcf.getCF() #get onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking
st_copy[0].data = tr_arzaic.data aicpick = AICPicker(aiccf.getCF(), taiccf, tr_filt.stats.delta, 2, 70, [1, 0.5, 0.2], 2, 0.001, 2.5)
araiccf = AICcf(st_copy, cuttimes, t2, p) #instance of AICcf ##############################################################
#calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction
#get stream object of filtered data
st_copy[0].data = tr_filt.data
arzcf = ARZcf(st_copy, cuttimes, tpredz, arzorder, tdetz, addnoise) #instance of ARZcf
##############################################################
#calculate AIC-ARZ-CF using subclass AICcf of class CharacteristicFunction
#class needs stream object => build it
tr_arzaic = tr_filt.copy()
tr_arzaic.data = arzcf.getCF()
st_copy[0].data = tr_arzaic.data
araiccf = AICcf(st_copy, cuttimes, t2, p) #instance of AICcf
elif not wfzfiles: elif not wfzfiles:
print 'No vertical component data found!' print 'No vertical component data found!'
if wfefiles and wfnfiles: if wfefiles and wfnfiles:
for i in range(len(wfefiles)): for i in range(len(wfefiles)):
print 'Horizontal component data found ...' print 'Horizontal component data found ...'
print wfefiles[i] print wfefiles[i]
print wfnfiles[i] print wfnfiles[i]
#merge streams #merge streams
H = read('%s' % wfefiles[i]) H = read('%s' % wfefiles[i])
H += read('%s' % wfnfiles[i]) H += read('%s' % wfnfiles[i])
H_copy = H.copy() H_copy = H.copy()
#filter and taper data #filter and taper data
trH1_filt = H[0].copy() trH1_filt = H[0].copy()
trH2_filt = H[1].copy() trH2_filt = H[1].copy()
trH1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) trH1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False)
trH2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) trH2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False)
trH1_filt.taper(max_percentage=0.05, type='hann') trH1_filt.taper(max_percentage=0.05, type='hann')
trH2_filt.taper(max_percentage=0.05, type='hann') trH2_filt.taper(max_percentage=0.05, type='hann')
H_copy[0].data = trH1_filt.data H_copy[0].data = trH1_filt.data
H_copy[1].data = trH2_filt.data H_copy[1].data = trH2_filt.data
############################################################## ##############################################################
#calculate ARH-CF using subclass ARHcf of class CharcteristicFunction #calculate ARH-CF using subclass ARHcf of class CharcteristicFunction
arhcf = ARHcf(H_copy, cuttimes, tpredh, arhorder, tdeth, addnoise) #instance of ARHcf arhcf = ARHcf(H_copy, cuttimes, tpredh, arhorder, tdeth, addnoise) #instance of ARHcf
############################################################## ##############################################################
#create stream with 3 traces #create stream with 3 traces
#merge streams #merge streams
AllC = read('%s' % wfefiles[i]) AllC = read('%s' % wfefiles[i])
AllC += read('%s' % wfnfiles[i]) AllC += read('%s' % wfnfiles[i])
AllC += read('%s' % wfzfiles[i]) AllC += read('%s' % wfzfiles[i])
#filter and taper data #filter and taper data
All1_filt = AllC[0].copy() All1_filt = AllC[0].copy()
All2_filt = AllC[1].copy() All2_filt = AllC[1].copy()
All3_filt = AllC[2].copy() All3_filt = AllC[2].copy()
All1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) All1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False)
All2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) All2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False)
All3_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False) All3_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False)
All1_filt.taper(max_percentage=0.05, type='hann') All1_filt.taper(max_percentage=0.05, type='hann')
All2_filt.taper(max_percentage=0.05, type='hann') All2_filt.taper(max_percentage=0.05, type='hann')
All3_filt.taper(max_percentage=0.05, type='hann') All3_filt.taper(max_percentage=0.05, type='hann')
AllC[0].data = All1_filt.data AllC[0].data = All1_filt.data
AllC[1].data = All2_filt.data AllC[1].data = All2_filt.data
AllC[2].data = All3_filt.data AllC[2].data = All3_filt.data
#calculate AR3C-CF using subclass AR3Ccf of class CharacteristicFunction #calculate AR3C-CF using subclass AR3Ccf of class CharacteristicFunction
ar3ccf = AR3Ccf(AllC, cuttimes, tpredz, arhorder, tdetz, addnoise) #instance of AR3Ccf ar3ccf = AR3Ccf(AllC, cuttimes, tpredz, arhorder, tdetz, addnoise) #instance of AR3Ccf
############################################################## ##############################################################
if iplot: if iplot:
#plot vertical trace #plot vertical trace
plt.figure() plt.figure()
tr = st[0] tr = st[0]
tstepz = tpredz / 16 tstepz = tpredz / 16
tdata = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta) tdata = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta)
thoscf = np.arange(0, len(hoscf.getCF()) / tr.stats.sampling_rate, tr.stats.delta) + cuttimes[0] tarzcf = np.arange(0, len(arzcf.getCF()) * tstepz, tstepz) + cuttimes[0] + tdetz +tpredz
taiccf = np.arange(0, len(aiccf.getCF()) / tr.stats.sampling_rate, tr.stats.delta) + cuttimes[0] taraiccf = np.arange(0, len(araiccf.getCF()) * tstepz, tstepz) + cuttimes[0] +tdetz + tpredz
tarzcf = np.arange(0, len(arzcf.getCF()) * tstepz, tstepz) + cuttimes[0] + tdetz +tpredz p1 = plt.plot(tdata, tr_filt.data/max(tr_filt.data), 'k')
taraiccf = np.arange(0, len(araiccf.getCF()) * tstepz, tstepz) + cuttimes[0] +tdetz + tpredz p2 = plt.plot(thoscf, hoscf.getCF()/max(hoscf.getCF()), 'r')
p1 = plt.plot(tdata, tr_filt.data/max(tr_filt.data), 'k') p3 = plt.plot(taiccf, aiccf.getCF()/max(aiccf.getCF()), 'b')
p2 = plt.plot(thoscf, hoscf.getCF()/max(hoscf.getCF()), 'r') p4 = plt.plot(tarzcf, arzcf.getCF()/max(arzcf.getCF()), 'g')
p3 = plt.plot(taiccf, aiccf.getCF()/max(aiccf.getCF()), 'b') p5 = plt.plot(taraiccf, araiccf.getCF()/max(araiccf.getCF()), 'y')
p4 = plt.plot(tarzcf, arzcf.getCF()/max(arzcf.getCF()), 'g') plt.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1], 'b')
p5 = plt.plot(taraiccf, araiccf.getCF()/max(araiccf.getCF()), 'y') plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [1, 1], 'b')
plt.yticks([]) plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [-1, -1], 'b')
plt.xlabel('Time [s]') plt.yticks([])
plt.ylabel('Normalized Counts') plt.xlabel('Time [s]')
plt.title([tr.stats.station, tr.stats.channel]) plt.ylabel('Normalized Counts')
plt.suptitle(tr.stats.starttime) plt.title([tr.stats.station, tr.stats.channel])
plt.legend([p1, p2, p3, p4, p5], ['Data', 'HOS-CF', 'HOSAIC-CF', 'ARZ-CF', 'ARZAIC-CF']) plt.suptitle(tr.stats.starttime)
#plot horizontal traces plt.legend([p1, p2, p3, p4, p5], ['Data', 'HOS-CF', 'HOSAIC-CF', 'ARZ-CF', 'ARZAIC-CF'])
plt.figure(2) #plot horizontal traces
plt.subplot(211) plt.figure(2)
tsteph = tpredh / 4 plt.subplot(211)
th1data = np.arange(0, trH1_filt.stats.npts / trH1_filt.stats.sampling_rate, trH1_filt.stats.delta) tsteph = tpredh / 4
th2data = np.arange(0, trH2_filt.stats.npts / trH2_filt.stats.sampling_rate, trH2_filt.stats.delta) th1data = np.arange(0, trH1_filt.stats.npts / trH1_filt.stats.sampling_rate, trH1_filt.stats.delta)
tarhcf = np.arange(0, len(arhcf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdeth +tpredh th2data = np.arange(0, trH2_filt.stats.npts / trH2_filt.stats.sampling_rate, trH2_filt.stats.delta)
p21 = plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k') tarhcf = np.arange(0, len(arhcf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdeth +tpredh
p22 = plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r') p21 = plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k')
plt.yticks([]) p22 = plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r')
plt.ylabel('Normalized Counts') plt.yticks([])
plt.title([trH1_filt.stats.station, trH1_filt.stats.channel]) plt.ylabel('Normalized Counts')
plt.suptitle(trH1_filt.stats.starttime) plt.title([trH1_filt.stats.station, trH1_filt.stats.channel])
plt.legend([p21, p22], ['Data', 'ARH-CF']) plt.suptitle(trH1_filt.stats.starttime)
plt.subplot(212) plt.legend([p21, p22], ['Data', 'ARH-CF'])
p23 = plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k') plt.subplot(212)
p24 = plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r') p23 = plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k')
plt.title([trH2_filt.stats.station, trH2_filt.stats.channel]) p24 = plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r')
plt.yticks([]) plt.title([trH2_filt.stats.station, trH2_filt.stats.channel])
plt.xlabel('Time [s]') plt.yticks([])
plt.ylabel('Normalized Counts') plt.xlabel('Time [s]')
#plot 3-component window plt.ylabel('Normalized Counts')
plt.figure(3) #plot 3-component window
tar3ccf = np.arange(0, len(ar3ccf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdetz +tpredz plt.figure(3)
plt.subplot(311) tar3ccf = np.arange(0, len(ar3ccf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdetz +tpredz
p31 = plt.plot(tdata, tr_filt.data/max(tr_filt.data), 'k') plt.subplot(311)
p32 = plt.plot(tar3ccf, ar3ccf.getCF()/max(ar3ccf.getCF()), 'r') p31 = plt.plot(tdata, tr_filt.data/max(tr_filt.data), 'k')
plt.yticks([]) p32 = plt.plot(tar3ccf, ar3ccf.getCF()/max(ar3ccf.getCF()), 'r')
plt.xticks([]) plt.yticks([])
plt.ylabel('Normalized Counts') plt.xticks([])
plt.title([tr.stats.station, tr.stats.channel]) plt.ylabel('Normalized Counts')
plt.legend([p31, p32], ['Data', 'AR3C-CF']) plt.title([tr.stats.station, tr.stats.channel])
plt.subplot(312) plt.legend([p31, p32], ['Data', 'AR3C-CF'])
plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k') plt.subplot(312)
plt.plot(tar3ccf, ar3ccf.getCF()/max(ar3ccf.getCF()), 'r') plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k')
plt.yticks([]) plt.plot(tar3ccf, ar3ccf.getCF()/max(ar3ccf.getCF()), 'r')
plt.xticks([]) plt.yticks([])
plt.ylabel('Normalized Counts') plt.xticks([])
plt.title([trH1_filt.stats.station, trH1_filt.stats.channel]) plt.ylabel('Normalized Counts')
plt.subplot(313) plt.title([trH1_filt.stats.station, trH1_filt.stats.channel])
plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k') plt.subplot(313)
plt.plot(tar3ccf, ar3ccf.getCF()/max(ar3ccf.getCF()), 'r') plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k')
plt.yticks([]) plt.plot(tar3ccf, ar3ccf.getCF()/max(ar3ccf.getCF()), 'r')
plt.ylabel('Normalized Counts') plt.yticks([])
plt.title([trH2_filt.stats.station, trH2_filt.stats.channel]) plt.ylabel('Normalized Counts')
plt.xlabel('Time [s]') plt.title([trH2_filt.stats.station, trH2_filt.stats.channel])
plt.show() plt.xlabel('Time [s]')
raw_input() plt.show()
plt.close() raw_input()
plt.close()
parser = argparse.ArgumentParser() parser = argparse.ArgumentParser()
parser.add_argument('--project', type=str, help='project name (e.g. Insheim)') parser.add_argument('--project', type=str, help='project name (e.g. Insheim)')
parser.add_argument('--database', type=str, help='event data base (e.g. 2014.09_Insheim)') parser.add_argument('--database', type=str, help='event data base (e.g. 2014.09_Insheim)')