PSD (Power Spectral Density) measurement¶
This will compuete PSD for continuous seismic waveforms and requires ObsPy and Matplotlib
Import ObsPy module¶
In [69]:
from obspy import read, read_inventory
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal import PPSD
from obspy.imaging.cm import pqlx
import obspy as ob
print("# obspy version = ",ob.__version__)
Set size of figure¶
In [70]:
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 10,10
Set font size¶
In [71]:
# font size
SMALL_SIZE = 12
MEDIUM_SIZE = 14
BIGGER_SIZE = 16
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
Set client (Data Center)¶
This example uses NCEDC. We can use other dataceneter (e.g., SCEDC, IRIS, GEOFON...)
In [72]:
client = Client("NCEDC") # data from NCEDC
#client = Client("SCEDC") # data from SCEDC
#client = Client("IRIS") # data from IRIS
Set SNCL¶
Which SNCL (Station, Network, Component, Location)? This example uses BKS.BK.HHZ.00 data
In [73]:
# BKS BHZ data
sta = "BKS" # station
com = "HHZ" # componnet
net = "BK" # network
loc = "00" # location "--" for blank location code
Set time window¶
This example uses data collected between 2020/06/12/-2020/06/13
In [74]:
start_day = "2020-06-01"
end_day = "2020-06-03"
starttime = UTCDateTime(start_day)
endtime = UTCDateTime(end_day)
PSD time window¶
This example uses 7200-sec-long data
In [75]:
ppsd_length = 7200 # unit is seconds
Get instrument¶
use get_stations with level = "response"
In [76]:
inv = client.get_stations(network=net, station=sta,
starttime=starttime,
endtime=endtime, level="response")
print(inv)
Download seismic data¶
use get_waveforms
In [77]:
st = client.get_waveforms(net, sta, loc, com,
starttime, endtime)
print(st)
Plot seismic data¶
check waveforms
In [78]:
st.plot() # check waveform
Out[78]:
PSD estimate¶
use PPSD function from obsPy
In [79]:
# psd estimate
min_db = -200 # min. db
max_db = -50 # max db
ddb = 1 # 1 db increment
In [80]:
tr = st[0] # first segment only. if data have gaps, only first segment is stored
ppsd = PPSD(stats=tr.stats, metadata=inv, db_bins=(min_db, max_db, ddb), ppsd_length=ppsd_length)
In [81]:
ppsd.add(tr)
Out[81]:
PSD plot¶
In [82]:
period_low = 0.05 # 0.05 sec
period_max = 1000 # 1000 sec
In [83]:
# plot
ppsd.plot(cmap=pqlx, period_lim=(period_low, period_max), show_mean = True, )
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: