testst.py

from pyrtt.evstream import read_event
from obspy.taup import TauPyModel
from obspy import read_inventory
import pprint
import matplotlib.pyplot as plt
import glob
import objgraph
model = TauPyModel(model='iasp91')
print('created model')
#st = read_event('tstEvents/20110101_0956*/processed/*HZ',model=model)
#st = read_event('tstEvents/20130417_120331.a/processed/*HZ',model=model)
#st = read_event('tstEvents/20130419_030552.a/processed/*HZ',model=model)
#st = read_event('tstEvents/20130419_195841.a/p*/*HZ',model=model)
#st = read_event('tstEvents/20130420_131252.a/processed/*HZ',model=model)
#st = read_event('tstEvents/20130421_134830.a/processed/*HZ',model=model)  #Doesn't work well
#st = read_event('tstEvents/20130422_011632.a/processed/*HZ',model=model)   #Works eh
#st = read_event('tstEvents/20130422_043639.a/p*/*HZ',model=model) #Removes all
evs = glob.glob('tstEvents/20120[45]*')

def initial_process(st):
    st.clean_duplicates()
    st.detrend(type='demean')
    st.detrend(type='simple')
    st.taper(0.05)
    st.interpolate(20.0)

for ev in evs:
    print(ev)
    st = read_event(f'{ev}/processed/*HZ',model=model)

    fig,ax = plt.subplots(4,sharex=True)
    initial_process(st)
    #Ray theory filt
    st.filter_rtt(.1,2)
    print('filtered')
    st.envelope_rtt()
    print('envelopes')
    st.trim_rtt(10,10)
    print('trimmed 10 before t0 and 10 after t0')
    for tr in st:
        if len(tr.data) != 400:
            st.remove(tr)
     
    st.stack()
    st.snr()
    print(f'SNR: {st.info.stream_snr}')
    for i,tr in enumerate(st):
        tmp = tr.window_around_arrival(window=(5,5),arrival='t0')
        ax[0].plot(tmp)

   #ax[0].plot(tr.times(),tr.data)
    print(f'Before env Rejection: {len(st)}')
    st.envelope_rejection(8)
    print(f'After Env Rejection/Before Iter: {len(st)}')
    st.iter_env_rejection(1.5)
    print(f'After iter env: {len(st)}')
    for i,tr in enumerate(st):
        tmp = tr.window_around_arrival(window=(5,5),arrival='t0')
        ax[1].plot(tmp)
     
    st.iccs(cc_cutoff=0.78)

    if len(st) < 5:
        continue
    print(f'After iccs: {len(st)}')
#fig,ax = plt.subplots(2,sharex=True)
    for i,tr in enumerate(st):
        tmpcc = tr.window_around_arrival(window=(5,5),arrival='t1')
        tmp = tr.window_around_arrival(window=(5,5),arrival='t0')
    #print(i,tr.stats['envmax'])
#   ax[0].plot(st.info.stack_envelope)
        ax[2].plot(tmp)
        ax[3].plot(tmpcc)
    #print(tr.stats.t0, tr.stats.t1)
    #ax[0].axvline(tr.stats.t1)
    ax[0].set_title('Input Data')
    ax[1].set_title('After Envelope Rejection')
    ax[2].set_title('Theoretical Arrival')
    ax[3].set_title('Iccs around new arrival')
#vline = st[3].stats.t1
#ax[3].axvline(vline/st[3].stats.delta)
#plt.show()

    #plt.show()
    plt.title(f'{st[0].stats.origintime}  {st[0].stats.mag}')
    if len(st) > 5:
        plt.savefig(f'{ev}/processed/rejections.pdf')
        st.write_sac_for_mccc(path=f'{ev}/processed/',name='low',arrival_header='t1')

    fig,ax = plt.subplots(3,sharex=True)
    for tr in st:
        tmp = tr.window_around_arrival(window=(5,5),arrival='t1')
        ax[0].plot(tmp)

    plt.savefig(f'{ev}/processed/FF_rej.pdf')
    del st
    del fig
    del ax

objgraph.show_most_common_types()
#objgraph.show_backrefs(random.choice(objgraph.by_type()))