Ridgecrest Waveform request

This will read Thugman's event catalog and dowaload waveforms for each event. ObsPy, numpy, scipy, matplotlib are required.

Import ObsPy module

In [ ]:
from obspy import read
from obspy import UTCDateTime
from obspy import read, read_inventory
from obspy.clients.fdsn import Client

import obspy as ob
print("# obspy version = ",ob.__version__)

Import SciPy, NumPy, matplotlib, pandas module

In [ ]:
import numpy as np
import scipy as sp
import matplotlib as mpl
import pandas as pd

print("# numpy version = ",np.__version__)
print("# scipy version = ",sp.__version__)
print("# matplotlib version = ",mpl.__version__)
print("# pandas version = ",pd.__version__)

import matplotlib.pyplot as plt

import os

function: request_waveform

In [ ]:
def request_waveform(eventid, request_server, distdeg_from_eq, catalog, debugOPT):

    #checkOPT = 1
    checkOPT = 0
    
    # name for output directory
    pwd_dir = os.getcwd() 
    # v1
    #outdir = pwd_dir +"/"+(str)(eventid) +"_dist"+(str)(distdeg_from_eq)+"deg"
    # v2
    outdir = pwd_dir +"/"+(str)(eventid)
    
    print("# outdir = ",outdir)

    # create output directory
    if not os.path.exists(outdir):
        os.makedirs(outdir)

    # get event info. based on eventif
    event_select = catalog [ (catalog ['evid'] == eventid) ]
    print("# event_select = ", event_select)

    index = event_select['evla'].index.to_numpy()
    if debugOPT:
        print("# index = ", index)

    latitude = event_select['evla'][index].to_numpy()
    longitude = event_select['evlo'][index].to_numpy()
    if debugOPT:
        print("# latitude = ", latitude," longitude = ", longitude)


    
    origin_time = event_select['origin_time'][index].to_numpy()
    #print("# origin_time = ",  origin_time)
    
    origin_time = UTCDateTime( origin_time[0]  )
    if debugOPT:
        print("# origin_time = ",  origin_time)

    starttime_tw = UTCDateTime(origin_time-pre_tw_sec)
    endtime_tw = UTCDateTime(origin_time+tw_sec)

    if debugOPT:
        print("# starttime_tw = ",starttime_tw," endtime_tw = ",endtime_tw)

           
    ### v2
    yearOUT = origin_time.year
    jdayOUT = origin_time.julday
    hourOUT = origin_time.hour
    minOUT = origin_time.minute
    secOUT = origin_time.second

    yearOUT2 = (str)(yearOUT)

    jdayOUT2 = (str)(jdayOUT)
    if jdayOUT < 100:
        jdayOUT2 = "0"+(str)(jdayOUT)


    if jdayOUT < 10:
        jdayOUT2 = "00"+(str)(jdayOUT)
        

    hourOUT2 = (str)(hourOUT)
    if hourOUT < 10:
        hourOUT2 = "0"+(str)(hourOUT)


    minOUT2 = (str)(minOUT)
    if minOUT < 10:
        minOUT2 = "0"+(str)(minOUT)


    secOUT2 = (str)(secOUT)
    if secOUT < 10:
        secOUT2 = "0"+(str)(secOUT)


    mseedid = yearOUT2+"."+jdayOUT2+"."+hourOUT2+""+minOUT2+""+secOUT2

    if debugOPT:
        print("# yearOUT2 = ",yearOUT2," jdayOUT2 = ",jdayOUT2," hourOUT2 = ",hourOUT2," minOUT2 = ",minOUT2," secOUT2 = ",secOUT2)
        print("# mseedid = "+mseedid)
    ### v2
        
    

    # first set
    # Without nodal first (network 3J)
    net_listOUT = "*"
    sta_list = "*"
    
    if checkOPT:
        net_listOUT = "CI"
        sta_list = "WMF"

    #com_listOUT = "HHZ,HHZ,DP1,EP1"
    #HNZ HLZ 
    #B? 20 or 40sps 
    #M>=3 or 4 EQ -> only B stations 
    # v1 vertical only
    com_listOUT = "H?Z,E?Z,D?Z,E?Z,C?Z"
    # v2 includ all components
    com_listOUT = "H??,E??,D??,E??,C??"
    if checkOPT:
        com_listOUT = "HHZ"
        
    #com_listOUT = "H?Z,H?1,E?Z,E?1,D?Z,D?1,E?Z,E?1,C?Z,C?1"

    if request_server == 1:
        server_nameOUT = "iris"
        clientOUT = Client("IRIS", timeout=600)
    
    elif request_server == 2:
        server_nameOUT = "scedc"
        clientOUT = Client("SCEDC", timeout=600)

    elif request_server == 3:
        server_nameOUT = "ncedc"
        clientOUT = Client("NCEDC", timeout=600)

    if debugOPT:
        print(clientOUT)
    

    # search area
    minradius = 0
    maxradius = distdeg_from_eq
    
    # inventory 
    invOUT = clientOUT.get_stations(network=net_listOUT, station=sta_list, channel=com_listOUT, 
                                    starttime=starttime_tw, endtime=endtime_tw,
                                    latitude=latitude, longitude=longitude,
                                    minradius=minradius, maxradius=maxradius,
                                    level="response")
    # channel info.
    channelOUT = clientOUT.get_stations(network=net_listOUT, station=sta_list, channel=com_listOUT, 
                                    starttime=starttime_tw, endtime=endtime_tw,
                                    latitude=latitude, longitude=longitude,
                                    minradius=minradius, maxradius=maxradius,
                                    level="channel")
    #print(invOUT)
        
    # save station list
    #inv_fiOUT = outdir + "/" + (str)(eventid)+".sta."+(str)(distdeg_from_eq)+"deg."+server_nameOUT+".txt"
    inv_fiOUT = outdir + "/" + (str)(eventid)+"."+server_nameOUT+".txt"

    if debugOPT:
        print("# inv_fiOUT = "+inv_fiOUT)
    invOUT.write(inv_fiOUT, format='STATIONTXT')  
    
    # save station xml
    #inv_fiOUT = outdir + "/" + (str)(eventid)+".sta."+(str)(distdeg_from_eq)+"deg."+server_nameOUT+".xml"
    inv_fiOUT = outdir + "/" + (str)(eventid)+"."+server_nameOUT+".xml"

    if debugOPT:
        print("# inv_fiOUT = "+inv_fiOUT)
    invOUT.write(inv_fiOUT, format='STATIONXML')  

    # save sac pz
    #inv_fiOUT = outdir + "/" + (str)(eventid)+".sta."+(str)(distdeg_from_eq)+"deg."+server_nameOUT+".sacpz"
    inv_fiOUT = outdir + "/" + (str)(eventid)+"."+server_nameOUT+".sacpz"

    if debugOPT:
        print("# inv_fiOUT = "+inv_fiOUT)
    invOUT.write(inv_fiOUT, format='SACPZ')  

    
    # get sncl list
    ch_listOUT = ""
    j2=0
    sncl_list = []
    for i in range(len(channelOUT)):        
        netOUT = channelOUT[i].code
        #print (inv_ncedc[i].code)
        for ii in range(len(channelOUT[i])):
            #print(ii)
            #print(inv[i][ii].code)    
            staOUT = channelOUT[i][ii].code
            for iii in range(len(channelOUT[i][ii])):
                chaOUT = channelOUT[i][ii][iii].code 
                locOUT = channelOUT[i][ii][iii].location_code
                snclOUT = netOUT+"."+staOUT+"."+locOUT+"."+chaOUT
                #print("# j2 = ",j2," snclOUT = ",snclOUT)
                sncl_list.append(snclOUT)
                j2+=1

                

    # get station list
    sta_listOUT = ""
    j=0
    for i in range(len(invOUT)):
        #print (inv_ncedc[i].code)
        for ii in range(len(invOUT[i])):
            #print(ii)
            #print(inv[i][ii].code)
            if j==0:
                    sta_listOUT = invOUT[i][ii].code
            else:
                    sta_listOUT = sta_listOUT+","+invOUT[i][ii].code
            j+=1



    print("# sta_listOUT = ",sta_listOUT)


    # get seismic waveforms nased on station & network lists
    stOUT = clientOUT.get_waveforms(network=net_listOUT, station=sta_listOUT, channel=com_listOUT, location="*",
                                    starttime=starttime_tw, endtime=endtime_tw)
    #CI.WMF..HHZ
    #stOUT = clientOUT.get_waveforms(network="CI", station="WMF", channel="BHZ", location="*",
    #                                starttime=starttime_tw, endtime=endtime_tw)    
    
    
    if debugOPT:
        print(stOUT)
    


    #for i in range(len(stOUT)):
    for tr in stOUT:
        if debugOPT:
            print(tr)

        netOUT = tr.stats.network
        comOUT = tr.stats.channel
        staOUT = tr.stats.station
        locOUT = tr.stats.location

        if netOUT == "3J":
            # nodal data will be excluded
            continue

        stid = netOUT+"."+staOUT+"."+locOUT+"."+comOUT
        sncl_id = staOUT+"."+netOUT+"."+comOUT+"."+locOUT

        # v1
        #mseed_fiOUT = outdir+"/"+sncl_id+".ms"
        # v2
        mseed_fiOUT = outdir+"/"+sncl_id+"."+mseedid+"."+str(eventid)+".ms"
        
        if debugOPT:
            print("# staOUT = ",staOUT," netOUT = ",netOUT," comOUT = ",comOUT," locOUT = ",locOUT)
            print("# mseed_fiOUT = ", mseed_fiOUT)
        #print("# stid = "+stid)
        tr.write(mseed_fiOUT, format="MSEED") 
        
    return stOUT
    
    

debugOPT

In [ ]:
# check output values 
debugOPT = 1
#debugOPT = 0

plotOPT

In [ ]:
# show waveforms
plotOPT = 0

Reading relocated catalog

use Daniel Trugman's reolocated catalog.

In [ ]:
#cat Dataset\ S1.txt | awk '{print $1","$2","$3","$4","$5","$6","$7","$8}' | less >! Dataset_S1.csv
#catalog_fi = "/Users/taira/Downloads/bssa-2020009_supplement_datasets+s1+and+s2/Dataset_S1.csv"
catalog_fi = "./Dataset_S1.csv"
print("# catalog_fi = ", catalog_fi)

#Dataset S1. Relocated earthquake catalog. Each row corresponds to an earthquake, and there are 7 columns in total: 
#1.	Event ID
#2.	Origin Time
#3.	Magnitude (SCEDC preferred)
#4.	Longitude
#5.	Latitude
#6.	Depth (km)
#7.	Relocated (0 = No, 1 = Yes) -> locid


catalog = pd.read_csv(catalog_fi, skiprows=None,
                       sep=",",names=["evid", "origin_time1","origin_time2", "scedc_mag", "evlo", "evla", "dep", "locid"],header=None)

catalog['origin_time'] = catalog['origin_time1']+"T"+catalog['origin_time2']
catalog['time'] = pd.to_datetime(catalog['origin_time'])
print(catalog.head())

Selecting events

minimum magnitude is 1.0. Duration is two week from July 4, 2019 at 00:00.00 UTC

In [ ]:
# minimum magnitude
minmag = 1.0 

# two-week time window 
start_time = "2019-07-04T00:00:00"
end_time = "2019-07-18T00:00:00"


# half hour time window
# beta test 04.14.21
start_time = "2019-07-13T01:00:00"
end_time = "2019-07-13T01:30:00"

# select events based on magnitude and origin time
catalog_select = catalog[ (catalog['time']  >=  start_time) & (catalog['time']  <=  end_time) & (catalog['scedc_mag']  >=  minmag)  ]

print(catalog_select.head())
In [ ]:
#print(catalog_select)
In [ ]:
#len(catalog_select)
In [ ]:
 

Waveform request parameters

In [ ]:
# request server list
# 1 -> IRIS
# 2 -> SCEDC
# 3 -> NCEDC
request_servers = [3, 1, 2] # NCEDC -> IRIS -> SCEDC
# SCEDC test
#request_servers = [2] # SCEDC

# distance from hypocenter in degrees
distdeg_from_eq = 1.0 # -> ~110km

# time window from origin time to extract waveform
# from BE
#You might consider increasing the length to 2 minutes (10 before to 110 after), but I don’t feel strongly about it.
# from RA -> 90 seconds
# coda?  55min! 4.0<=M 300km

# time window from event origin time
pre_tw_sec = 0 # 0 seconds before the origin time
tw_sec = 90 # 90 sec after the origin time
In [ ]:
#len(catalog_select['evla'].index.to_numpy())
#catalog_select['evla']
In [ ]:
testOPT = 1 # -> only one event
testOPT = 0

selectEQ_index = catalog_select['evla'].index.to_numpy()
print("# selectEQ_index = ", selectEQ_index)

for i in range(len(selectEQ_index)):
    i2 = selectEQ_index[i]
    eventid = catalog_select['evid'][i2]
    print("# eventid = ", eventid)
    
    # RC mainshock 
    #eventid = 38457511 # 

    for request_server in request_servers:
        print("# request_server = ", request_server)
        stOUT = request_waveform(eventid, request_server, distdeg_from_eq, catalog, debugOPT)

    if (i == 0) and testOPT:
        break
        

Showing waveforms from the last request server

In [ ]:
print(stOUT)
In [ ]:
if plotOPT: 
    for tr in stOUT:
        _plot = tr.plot()
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: