Waveform request¶
This will dowaload waveforms for one event and will save them in SAC format. ObsPy is required.
Import ObsPy module¶
In [1]:
from obspy import read
from obspy import UTCDateTime
from obspy import read, read_inventory
from obspy.clients.fdsn import Client
from obspy.io.sac.sactrace import SACTrace
import obspy as ob
print("# obspy version = ",ob.__version__)
# to make directory
import os
Plot Option¶
if 1, plotting waveforms etc
In [2]:
plotOPT = 1 # plotting waveforms etc
Set client (Data Center) for event search¶
This example uses USGS. We can use other dataceneter (e.g., NCEDC, IRIS)
In [3]:
#clientEQ = Client("IRIS")
clientEQ = Client("USGS")
Event search¶
Set time window (st, et) and range of magunitude (minmag, maxmag). This example will find the 2018 Mw4.4 Berkeley earthquake with obspy get_events function
In [4]:
#https://earthquake.usgs.gov/earthquakes/eventpage/nc72948801/executive
#M 4.4 - 2km SE of Berkeley, CA
#2018-01-04 10:39:37 (UTC)37.855°N 122.257°W12.3 km depth
st = UTCDateTime("2018-01-04T10:39:30") #
et = UTCDateTime("2018-01-04T10:39:40") #
minmag = 4.3
maxmag = 4.5
catalog = clientEQ.get_events(starttime=st , endtime=et,
minmagnitude=minmag, maxmagnitude=maxmag)
print(catalog)
if plotOPT:
_plot = catalog.plot()
Set client (Data Center)¶
This example will use NCEDC. We can use other dataceneter (e.g., SCEDC, IRIS...)
In [5]:
server_nameOUT = "NCEDC" # data from NCEDC
#server_nameOUT = "SCEDC" # data from SCEDC
#server_nameOUT = "IRIS" # data from IRIS
client = Client(server_nameOUT, timeout=600) # data from NCEDC
Distance threshold to search for stations¶
This example will find all stations within 1.0 deg (~111km) from the epicenter
In [6]:
distdeg_from_eq = 1.00 # ~111km from epicenter
Time window¶
This example will extract data starting from -90s before the origin time and its duration is 5 minutes
In [7]:
pre_tw_sec = 90 # 90 seconds before the origin time
tw_sec = 5*60 # 5 min duration
Which components?¶
In many cases: HN data are strong-motion data (highest sampling at each station) & HH data are broadband data (highest sampling at each station). This example asks HN? data (strong-motion data) only
In [8]:
com_listOUT = "HN?" # strong-motion
#com_listOUT = "HH?" # broadband
#com_listOUT = "HH?,HN?" # broadband & strong-motion
Which netowrk (network code)?¶
This example will use a wildcard. Some examples are BK (Berkeley Seismological Laboratory), NC (USGS-Menlo), CE (CGS), NP (National Ground Motion), PB (PBO). more info. https://www.fdsn.org/networks/
In [9]:
#net_listOUT = "BK,NC,CE,NP,PB" #
net_listOUT = "*" # wildcard (all networks availanle at data center)
Which station (station code)?¶
Usually this should be a wildcard as we are searching all stations with the distance threhould. but for this example, we already identify which stations are needed
In [10]:
sta_list = "1688,1722,1737,1828,1856,58340,58347,58349,58393,58406,58497,58665,C002,C006,C007,C011,C013,C021,C024,C033,C048,C064,C067,CMC,CNI,CPM,CSL,CSPB,P181,RFSB,VAK"
#sta_list = "VAK" # test
# sta_list = "*" # wildcard
Frequency passband¶
Which frequency band will be used when the instrument response is corrected? This example uses a 0.5-40 Hz passband
In [11]:
# f1<f2<f3<f4. effective band is between f2 and f3
f1 = 0.1 # 0.1 Hz
f2 = 0.5 # 0.5 Hz
f3 = 40 # 40 Hz
f4 = 50 # 50 Hz
Event information¶
This example will only read the first event in catalog
In [12]:
event = catalog[0]
origin = event.origins[0]
origin_time = origin.time
evlat = origin.latitude
evlon = origin.longitude
evdp_km = origin.depth / 1000
evmag = event.magnitudes[0].mag
evyearOUT = origin_time.year
evjdayOUT = origin_time.julday
evhourOUT = origin_time.hour
evminOUT = origin_time.minute
evsecOUT = origin_time.second
evid = event.origins[0]['extra']['dataid']['value']
Output directory¶
Based on event-id and filtering parameter etc. This example will make 2018.004.103937_M4.38_nc72948801_fl0.5_fh40_dist1.0deg
In [13]:
# need for output directory
evyearOUT2 = (str)(evyearOUT)
evjdayOUT2 = (str)(evjdayOUT)
if evjdayOUT < 100:
evjdayOUT2 = "0"+(str)(evjdayOUT)
if evjdayOUT < 10:
evjdayOUT2 = "00"+(str)(evjdayOUT)
evhourOUT2 = (str)(evhourOUT)
if evhourOUT < 10:
evhourOUT2 = "0"+(str)(evhourOUT)
evminOUT2 = (str)(evminOUT)
if evminOUT < 10:
evminOUT2 = "0"+(str)(evminOUT)
evsecOUT2 = (str)(evsecOUT)
if evsecOUT < 10:
evsecOUT2 = "0"+(str)(evsecOUT)
#print("# evyearOUT2 = ",evyearOUT2," evjdayOUT2 = ",evjdayOUT2," evhourOUT2 = ",evhourOUT2," evminOUT2 = ",evminOUT2," evsecOUT2 = ",evsecOUT2)
evmseedid = evyearOUT2+"."+evjdayOUT2+"."+evhourOUT2+""+evminOUT2+""+evsecOUT2
#print("# evmseedid = "+evmseedid)
mseedid = evyearOUT2+"."+evjdayOUT2+"."+evhourOUT2+""+evminOUT2+""+evsecOUT2
# name for output directory
pwd_dir = os.getcwd()
#sacdir= pwd_dir +"/"+ evmseedid +"_M"+(str)(evmag)+"_"+(str)(eventid_ncedc)+"_fl"+(str)(f2)+"_fh"+str(f3)+"_dist"+(str)(distkm_from_eq)+"km"
sacdir= pwd_dir +"/"+ evmseedid +"_M"+(str)(evmag)+"_"+(str)(evid)+"_fl"+(str)(f2)+"_fh"+str(f3)+"_dist"+(str)(distdeg_from_eq)+"deg"
print("# sacdir = ",sacdir)
# create output directory
if not os.path.exists(sacdir):
os.makedirs(sacdir)
Extracting inventory (instrument response) and channel list based on search parameters¶
In [14]:
latitude = evlat
longitude = evlon
minradius = 0
maxradius = distdeg_from_eq
starttime = UTCDateTime(origin_time-pre_tw_sec)
endtime = UTCDateTime(origin_time+tw_sec)
invOUT = client.get_stations(network=net_listOUT, station=sta_list, channel=com_listOUT,
starttime=starttime, endtime=endtime,
latitude=latitude, longitude=longitude,
minradius=minradius, maxradius=maxradius,
level="response")
channelOUT = client.get_stations(network=net_listOUT, station=sta_list, channel=com_listOUT,
starttime=starttime, endtime=endtime,
latitude=latitude, longitude=longitude,
minradius=minradius, maxradius=maxradius,
level="channel")
Station list from our search¶
In [15]:
# get station 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
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)
Extracting waveforms based on station list¶
In [16]:
# get seismic waveforms nased on station & network lists
stOUT = client.get_waveforms(network=net_listOUT, station=sta_listOUT, channel=com_listOUT, location="*",
starttime=starttime, endtime=endtime)
stOUT2 = stOUT.copy()
#st2 = stOUT2.select(network="NC")
# check if data are gappy. short segment (less than 1-sec-long) will be removed.
for trace in stOUT2:
nptsOUT = trace.stats.npts
deltaOUT = trace.stats.delta
sampHZOUT = 1.0/deltaOUT
#print("sampHZOUT = ", sampHZOUT, "nptsOUT = ", nptsOUT)
if nptsOUT < sampHZOUT:
stOUT2.remove(trace)
print(stOUT2.__str__(extended=True))
Correcting instrument response¶
This example will provide ground acceleration (m/s). If velocity data (m/s) is needed -> "ACC", if displacement data (m) is needed -> "DISP"
In [17]:
# remove instrument response
#output_unut = "VEL" # m/s
output_unit = "ACC" # m/s*2
#output_unit = "DISP" # m
stOUT2.remove_response(inventory=invOUT, pre_filt=(f1, f2, f3, f4), output=output_unit)
Out[17]:
Plot one waveform¶
In [18]:
if plotOPT:
#_plot = stOUT2.plot()
_plot = stOUT2[0].plot()
Save waveforms in SAC format¶
In [19]:
#st2 = stOUT2.select(network="NC",station="N*", location="10")
#print(st2.__str__(extended=True))
#print("# remove_response server_nameOUT = ",server_nameOUT)
#st2.remove_response(inventory=invOUT,
# pre_filt=(f1, f2, f3, f4),
# output="VEL")
# save waveforms in SAC format
for i in range(len(stOUT2)):
#print (st_ncedc2[i])
sac = SACTrace.from_obspy_trace(stOUT2[i])
sac.lcalda = True
netOUT = stOUT2[i].stats.network
comOUT = stOUT2[i].stats.channel
staOUT = stOUT2[i].stats.station
locOUT = stOUT2[i].stats.location
# st_starttime is different from each sacfile, but mseedid is used origin_time - pre_tw_sec
#st_starttime = st_ncedc2[i].stats.starttime
#print("# st_starttime = ",st_starttime)
stid = netOUT+"."+staOUT+"."+locOUT+"."+comOUT
print("# staOUT = ",staOUT," netOUT = ",netOUT," comOUT = ",comOUT," locOUT = ",locOUT," stid = ",stid)
#print("# stid = "+stid)
sta_coordinate = invOUT.get_coordinates(stid, starttime)
sac.stla = sta_coordinate['latitude']
sac.stlo = sta_coordinate['longitude']
sac.stel = sta_coordinate['elevation']
sac.evla = evlat
sac.evlo = evlon
sac.evdp = evdp_km
sac.kuser0 = "ncss"
sac.kevnm = (str)(evid)
sac.mag = evmag
outsacfi = sacdir +"/" + staOUT+"."+netOUT+"."+comOUT+"."+locOUT+".D."+mseedid+"."+(str)(evid)+".sac"
#print("# outsacfi = "+outsacfi)
#st_ncedc[i].write(outsacfi, format="SAC")
sac.write(outsacfi)
In [ ]:
In [ ]: