Waveform plot¶
This will download continuous seismic waveforms & plot them and requires ObsPy
Import ObsPy module¶
In [1]:
from obspy import read
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.signal.invsim import paz_2_amplitude_value_of_freq_resp
import obspy as ob
print("# obspy version = ",ob.__version__)
Import Math and NumPy module¶
In [2]:
import math as M
import numpy as np
print("# numpy version = ",np.__version__)
Making SmartSolo instrument response¶
In [3]:
# https://smartsolo.com/wp-content/uploads/2020/02/smartsolo_brochure_web_en_feb2020-v2.pdf
G_nodal= 76.7 # V/m/s
fc_nodal = 5.0 # corner frequency
damp_nodal = 0.7 # damping
In [4]:
# pole & zero
poles_nodal = [-(damp_nodal + M.sqrt(1 - damp_nodal ** 2) * 1j) * 2 * np.pi * fc_nodal,
-(damp_nodal - M.sqrt(1 - damp_nodal ** 2) * 1j) * 2 * np.pi * fc_nodal]
# should be 2 zeros for geophone
zeros_nodal = [0+0j, 0+0j]
In [5]:
# paz_nodalTMP. This is for computing A0. gain=1, sensitivity=1
paz_nodalTMP = {'poles': poles_nodal, 'zeros': zeros_nodal,'gain': 1, 'sensitivity':1 }
print(paz_nodalTMP)
In [6]:
# A0 is a normalization factor
normalized_frequency = 20 # should be above the coner period
A0_nodal = 1/(paz_2_amplitude_value_of_freq_resp(paz_nodalTMP, normalized_frequency))
print(A0_nodal)
Total sensitivity¶
currently the total sensitivity sets to be 1000 to match the BKS amplitude
In [7]:
pre_amp = 63.095734 # 36 dB # http://www.sengpielaudio.com/calculator-amplification.htm
FIR_gain = 1
#datalogger_gain = 16.9524751916 # (1069.62886533/ 63.09573/)
datalogger_gain = 16.9524751916 # (1069.62886533/ 63.09573/)
total_sensitivity = A0_nodal * G_nodal * pre_amp * datalogger_gain * FIR_gain
#total_sensitivity = A0_nodal * pre_amp * datalogger_gain * FIR_gain
### total_sensitivity is fixed to be 1000!
total_sensitivity = 1000
#print(total_sensitivity )
print("# pre_amp = ",pre_amp )
print("# G_nodal = ",G_nodal )
print("# datalogger_gain = ",datalogger_gain )
print("# FIR_gain = ",FIR_gain )
print("# total_sensitivity = ",total_sensitivity )
SmartSolo instrument response¶
In [8]:
paz_nodal = {'poles': poles_nodal, 'zeros': zeros_nodal,'gain': G_nodal, 'sensitivity':total_sensitivity }
print(paz_nodal)
Set client (Data Center)¶
This example uses IRIS. We can use other dataceneter (e.g., NCEDC, SCEDC)
In [9]:
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)?
In [10]:
# use Z, N or E
channel = "Z"
#channel = "N"
#channel = "E"
In [11]:
# BKS HNZ (EpiSensor) data as reference
sta = "BKS" # station
net = "BK" # network
com = "HN"+channel # component
loc = "00" # location
In [12]:
# separately added Nodal seismic data
sta2 = "1568" # station
net2 = "AA" # network # not used
com2 = "DP"+channel # component
loc2 = "00" # location # not used
#client2 = Client("IRIS") # data from IRIS
Set time window¶
This example uses 2-hour data for the 2020 M6.8 Philippines earthquake
In [13]:
#M 2.5 - 1km S of Pleasant Hill, CA
#2019-10-15 05:23:58 (UTC)37.938°N 122.060°W13.5 km depth
#https://earthquake.usgs.gov/earthquakes/eventpage/nc73291865/executive
start_day = "2019-10-15T05:23:58.50"
end_day = "2019-10-15T05:24:28.50"
mseed_time = "2019-10-15" # this needs for nodal mseed data.
starttime = UTCDateTime(start_day)
endtime = UTCDateTime(end_day)
Download BK.BKS EpiSensor data¶
use get_waveforms to download data and do st.plot() for plotting
In [14]:
# st is Cascadia data
st = client.get_waveforms(network=net, station=sta, location=loc, channel=com,
starttime=starttime, endtime=endtime,
attach_response=True)
_plot = st.plot()
Download SmartSolo seismic data¶
data are at ftp site
In [15]:
# st2 is nodal seismic data
#st_nodal_fi = dirTMP+"/B1482.DP"+channel+".2019-09-30T00:00:00.000.mseed"
#st_nodal_fi = dirTMP+"/B1482.DP"+channel+".2019-09-27T00:00:00.000.mseed"
dirTMP = "http://ncedc.org/ftp/outgoing/taira/Nodal/bk_test"
st2 = read(dirTMP+"/B"+sta2+"."+com2+"."+mseed_time+"T00:00:00.000.mseed")
# trim
st2.trim(starttime, endtime)
_plot = st2.plot()
Correct instrument response for BK.BKS¶
use remove_response to correct the instrument response. We can select output unit (displacement, velocity or accerelation)
In [16]:
st.detrend() # remove liner trend
st.detrend("demean") # demean
st.taper(0.05) # cosin taper
st = st.remove_response( output="VEL" ) # get velocity data (m/s)
#st = st.remove_response( output="DISP" ) # get displacement data (m)
#st = st.remove_response( output="ACC" ) # get acceleration data (m/s^2)
_plot = st.plot()
Correct instrument response for SmartSolo seismic data¶
use "decimate" (500sps -> 100sps) and "simulate" from obspy
In [17]:
# decimate to 100Hz
st2.decimate(factor=5, strict_length=False)
st2.simulate(paz_remove=paz_nodal, paz_simulate=None)
Out[17]:
Filtering and plot¶
Example uses 1-15 Hz passband
In [18]:
st_all = st2.copy() + st.copy()
#st_all = st2.copy() + st.copy()
st_all.detrend() # remove liner trend
st_all.detrend("demean") # demean
st_all.taper(0.05) # cosin taper
fl = 1 # in Hz
fh = 15 # in Hz
# do bandpass filter
st_all.filter(type='bandpass', freqmin=fl, freqmax=fh, corners=6, zerophase=True)
t1 = UTCDateTime("2019-10-15T05:24:00.00")
t2 = UTCDateTime("2019-10-15T05:24:15.00")
#_plot = st_all.plot(size=(800, 400))
_plot = st_all.plot(size=(1000, 400), starttime=t1, endtime=t2)
In [19]:
for tr in st_all:
print("# max. amplitude = ",tr.max())
tr.plot(starttime=t1, endtime=t2)
#tr.plot(size=(1200,300))
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: