In [18]:
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
#plt.style.use("bmh")
# figure size
plt.rcParams['figure.figsize'] = 15,15
In [19]:
# font size
SMALL_SIZE = 22
MEDIUM_SIZE = 24
BIGGER_SIZE = 32

#SMALL_SIZE = 32
#MEDIUM_SIZE = 32
#BIGGER_SIZE = 36

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
In [20]:
from obspy import read, read_inventory
In [21]:
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal import PPSD
from obspy.imaging.cm import pqlx
In [22]:
# data from NCEDC
client = Client("NCEDC")
In [23]:
# BKS BHZ data
sta = "BKS" 
com = "BHZ"
net = "BK"

start_day = "2020-03-01"
end_day = "2020-03-02"
In [24]:
starttime = UTCDateTime(start_day)
endtime = UTCDateTime(end_day)

# png (plot) file
png_name = net+"."+sta+"."+com+"."+start_day+"."+end_day+".png"
print(png_name)

BK.BKS.BHZ.2020-03-01.2020-03-02.png
In [25]:
# get response file
inv = client.get_stations(network=net, station=sta,
                                starttime=starttime,
                                endtime=endtime, level="response")
In [26]:
# download waveform
st = client.get_waveforms(net,  sta, "*", com,
                          starttime,  endtime)
print(st)

1 Trace(s) in Stream:
BK.BKS.00.BHZ | 2020-03-01T00:00:00.019539Z - 2020-03-01T23:59:59.994539Z | 40.0 Hz, 3456000 samples
In [27]:
# check waveform
st.plot()
Out[27]:
In [28]:
ppsd_length = 7200 #  2h
In [29]:
# first segment only. if data have gaps, only first segment is stored
tr = st[0]
In [30]:
# psd estimate
min_db = -200
max_db = -50
ddb = 1 # 1 db increment
ppsd = PPSD(stats=tr.stats, metadata=inv, db_bins=(min_db, max_db, ddb), ppsd_length=ppsd_length)
In [31]:
ppsd.add(tr)
Out[31]:
True
In [32]:
period_low = 0.05 # 0.05 sec
period_max = 1000 # 1000 sec
In [33]:
# plot 
ppsd.plot(cmap=pqlx, period_lim=(period_low, period_max), show_mean = True, )

/anaconda3/envs/netops/lib/python3.7/site-packages/obspy/signal/spectral_estimation.py:1830: MatplotlibDeprecationWarning: 
The set_clim function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use ScalarMappable.set_clim instead.
  cb.set_clim(*fig.ppsd.color_limits)
In [34]:
# save png file
ppsd.plot(png_name, cmap=pqlx, period_lim=(period_low, period_max), show_mean = True, )
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: