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

import obspy as ob
print("# obspy version = ",ob.__version__)
# obspy version =  1.2.2

Import Math and NumPy module

In [2]:
import math as M
import numpy as np
print("# numpy version = ",np.__version__)
# numpy version =  1.18.5

Set client (Data Center)

This example uses IRIS. We can use other dataceneter (e.g., NCEDC, SCEDC)

In [3]:
#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 [4]:
#http://ds.iris.edu/mda/NV/NSMTC/
sta = "NSMTC" # station
net = "NV" # network
com = "C?Z" # to get all SA-ULN and Geophone Z data
loc = "*" # use wildcard to get all location but later I only use CNZ.B2 for the metadata check
In [5]:
# separately added PGC data for comparison
#http://ds.iris.edu/mda/CN/PGC/--/HHZ/?starttime=2017-08-10T18:10:00&endtime=2599-12-31T23:59:59
sta2 = "PGC" # station
net2 = "CN" # network
com2 = "H?Z" # component # HHZ or HNZ. HHZ is broadband and HNZ is accelerometer
loc2 = "--" # location
client2 = Client("IRIS") # data from IRIS

Set time window

This example uses 2-hour data for the 2020 M6.8 Philippines earthquake

In [6]:
#M 8.2 - 101km SSW of Tres Picos, Mexico
#2017-09-08 04:49:19 (UTC)15.022°N 93.899°W47.4 km depth

start_day = "2017-09-08T04:49:19"
end_day = "2017-09-08T05:49:19"

starttime = UTCDateTime(start_day)
endtime = UTCDateTime(end_day)

Download seismic data collected at Cascadia experiment

use get_waveforms to download data and do st.plot() for plotting

In [7]:
# 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()

Copy original data

In [8]:
# st_mod is for an update metadata. we will use st_mod below
st_mod = st.copy()
st_mod2 = st.copy() # for geophone  fc=4.5 h=0.34 G=32V/m/s
st_mod3 = st.copy() # for geophone  fc=1.0 h=0.34 G=1.5V/m/s

Download seismic data collected at CN.PGC data

In [9]:
# st2 is referect PGC data
st2 = client2.get_waveforms(network=net2, station=sta2, location=loc2, channel=com2,
                     starttime=starttime, endtime=endtime, 
                     attach_response=True)
_plot = st2.plot()

Correct instrument response for NV.NSMTC

use remove_response to correct the instrument response. We can select output unit (displacement, velocity or accerelation)

In [10]:
# select B2.CNZ data here to compere the metadata
#st = st.select(id="NV.NSMTC.B2.CNZ")
st = st.select(id="NV.NSMTC.G2.CHZ")

st.detrend() # remove liner trend
st.detrend("demean") # demean
st.taper(0.05) # cosin taper
# decimate to 100Hz
st.decimate(factor=5, strict_length=False)

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()

Update channel name

In [11]:
print(st)
1 Trace(s) in Stream:
NV.NSMTC.G2.CHZ | 2017-09-08T04:49:19.000000Z - 2017-09-08T05:49:19.000000Z | 100.0 Hz, 360001 samples
In [12]:
# updated channel name so that when we plot data we can see the sensor gain value info.
#st[0].stats.channel = "CNZ-60V/g"
st[0].stats.channel = "CHZ-81.82V/m/s x -1 fc=8.0 h=0.7"

Correct instrument response for CN.PGC

In [13]:
st2.detrend() # remove liner trend
st2.detrend("demean") # demean
st2.taper(0.05) # cosin taper
st2 = st2.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 = st2.plot()

Select G2.CHZ data

In [14]:
# select G2.CHZ data only 
st_mod2 = st_mod2.select(id="NV.NSMTC.G2.CHZ")
print(st_mod2)
st_mod3 = st_mod3.select(id="NV.NSMTC.G2.CHZ")
print(st_mod3)
1 Trace(s) in Stream:
NV.NSMTC.G2.CHZ | 2017-09-08T04:49:19.000000Z - 2017-09-08T05:49:19.000000Z | 500.0 Hz, 1800001 samples
1 Trace(s) in Stream:
NV.NSMTC.G2.CHZ | 2017-09-08T04:49:19.000000Z - 2017-09-08T05:49:19.000000Z | 500.0 Hz, 1800001 samples

Making new pole & zero

"mod2" uses a set of parameters from the spechseet

In [15]:
# this should match -> http://ds.iris.edu/NRL/sensors/oyo_geospace/RESP.XX.NS387..SHZ.GS11D.4_5.380.NONE.32
G_mod2 = 32 # V/m/s (Open-Circuit Sensitivity)
fc_mod2 = 4.5 # corner frequency
damp_mod2 = 0.34 # damping (Open-Circuit Damping 34%)
In [16]:
# pole & zero
poles_mod2 = [-(damp_mod2 + M.sqrt(1 - damp_mod2 ** 2) * 1j) * 2 * np.pi * fc_mod2,
                -(damp_mod2 - M.sqrt(1 - damp_mod2 ** 2) * 1j) * 2 * np.pi * fc_mod2]


print("# original pole & zero = ", st_mod2[0].stats.response.response_stages[0]._poles)
# update pole & zero
st_mod2[0].stats.response.response_stages[0]._poles = poles_mod2
#type(st_mod2[0].stats.response.response_stages[0]._poles )
print("# updated pole & zero = ", st_mod2[0].stats.response.response_stages[0]._poles)


print("# original gain value = ",st_mod2[0].stats.response.response_stages[0].stage_gain)
# update sensitivity value
# also x -1 to polarity flip
st_mod2[0].stats.response.response_stages[0].stage_gain =  G_mod2 * -1
print("# updated  gain value = ",st_mod2[0].stats.response.response_stages[0].stage_gain)
# original pole & zero =  [(-35.1858+35.8967j), (-35.1858-35.8967j)]
# updated pole & zero =  [(-9.61327351998477-26.589902758720232j), (-9.61327351998477+26.589902758720232j)]
# original gain value =  81.82
# updated  gain value =  -32

"mod3" uses a set of parameters that provide G2.CHZ waveform fitting better to CN.PGC data

In [17]:
# this set of parameters provide G2.CHZ waveform fitting better to CN.PGC data
# but more analysis required to update the metadata"mod2" uses a set of parameters from the spechseet
G_mod3 = 1.5 # V/m/s
fc_mod3 = 1.0 # corner frequency
damp_mod3 = 0.34 # damping
In [18]:
# pole & zero
poles_mod3 = [-(damp_mod3 + M.sqrt(1 - damp_mod3 ** 2) * 1j) * 2 * np.pi * fc_mod3,
                -(damp_mod3 - M.sqrt(1 - damp_mod3 ** 2) * 1j) * 2 * np.pi * fc_mod3]


print("# original pole & zero = ", st_mod3[0].stats.response.response_stages[0]._poles)
# update pole & zero
st_mod3[0].stats.response.response_stages[0]._poles = poles_mod3
#type(st_mod3[0].stats.response.response_stages[0]._poles )
print("# updated pole & zero = ", st_mod3[0].stats.response.response_stages[0]._poles)


print("# original gain value = ",st_mod3[0].stats.response.response_stages[0].stage_gain)
# update sensitivity value
# also x -1 to polarity flip
st_mod3[0].stats.response.response_stages[0].stage_gain =  G_mod3 * -1
print("# updated  gain value = ",st_mod3[0].stats.response.response_stages[0].stage_gain)
# original pole & zero =  [(-35.1858+35.8967j), (-35.1858-35.8967j)]
# updated pole & zero =  [(-2.1362830044410597-5.908867279715607j), (-2.1362830044410597+5.908867279715607j)]
# original gain value =  81.82
# updated  gain value =  -1.5

Correct instrument response for G2.CHZ data with updated resp files

"mod2" uses a set of parameters from the spechseet

In [19]:
# remove instrument response for st_mod
st_mod2.detrend() # remove liner trend
st_mod2.detrend("demean") # demean
st_mod2.taper(0.05) # cosin taper

# decimate to 100Hz
st_mod2.decimate(factor=5, strict_length=False)
st_mod2 = st_mod2.remove_response( output="VEL" ) # get velocity data (m/s)

_plot = st_mod2.plot()

"mod3" uses a set of parameters that provide G2.CHZ waveform fitting better to CN.PGC data

In [20]:
st_mod3.detrend() # remove liner trend
st_mod3.detrend("demean") # demean
st_mod3.taper(0.05) # cosin taper

# decimate to 100Hz
st_mod3.decimate(factor=5, strict_length=False)
st_mod3 = st_mod3.remove_response( output="VEL" ) # get velocity data (m/s)

_plot = st_mod3.plot()

Update channel name and location code

In [21]:
# updated channel name so that when we plot data we can see the sensor gain value info.
# location code is B2 -> C2 so that tihs treace would be plotted before G2 data. only plotting purpose
st_mod2[0].stats.channel = "CHZ-32V/m/s x -1 fc=4.5 h=0.34"
st_mod2[0].stats.location = "C2"
In [22]:
# updated channel name so that when we plot data we can see the sensor gain value info.
# location code is B2 -> C2 so that tihs treace would be plotted before G2 data. only plotting purpose
st_mod3[0].stats.channel = "CHZ-1.5V/m/s x -1 fc=1.0 h=0.34"
st_mod3[0].stats.location = "C2"

Filtering and plot

In [23]:
st_all = st2.copy() + st_mod2.copy() + st_mod3.copy ()+ st.copy()

st_all.detrend() # remove liner trend
st_all.detrend("demean") # demean
st_all.taper(0.05) # cosin taper

fl = 0.05 # in Hz 
fh = 0.10 # in Hz

fl = 0.10 # in Hz 
fh = 0.50 # in Hz


st_all.filter(type='bandpass', freqmin=fl, freqmax=fh, corners=6, zerophase=True)


### plotting the first arrival part ###
t1 = UTCDateTime("2017-09-08T04:56:00.00")
t2 = UTCDateTime("2017-09-08T05:00:00.00")
#_plot = st_all.plot(size=(800, 800))
_plot = st_all.plot(size=(1000, 700), starttime=t1, endtime=t2)
### plotting the first arrival part ###

### plotting surace wave(?)  part ###
t1 = UTCDateTime("2017-09-08T05:10:00.00")
t2 = UTCDateTime("2017-09-08T05:25:00.00")
#_plot = st_all.plot(size=(800, 800))
_plot = st_all.plot(size=(1000, 700), starttime=t1, endtime=t2)
### plotting surace wave(?) part ###
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: