This will read Thugman's event catalog and dowaload waveforms for each event. ObsPy, numpy, scipy, matplotlib are required.
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 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
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
# check output values
debugOPT = 1
#debugOPT = 0
# show waveforms
plotOPT = 0
use Daniel Trugman's reolocated catalog.
#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())
minimum magnitude is 1.0. Duration is two week from July 4, 2019 at 00:00.00 UTC
# 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())
#print(catalog_select)
#len(catalog_select)
# 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
#len(catalog_select['evla'].index.to_numpy())
#catalog_select['evla']
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
print(stOUT)
if plotOPT:
for tr in stOUT:
_plot = tr.plot()