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.1

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=2.0 h=0.70 G=76.2 V/m/s IRIS NRL HS-1-LT
st_mod3 = st.copy() # for geophone  fc=2.0 h=0.61 G=60 V/m/s. fc & h are the same as current metadata. G=60 to match the amplitude

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 G1.CHZ data here to compere the metadata
#st = st.select(id="NV.NSMTC.B2.CNZ")
#st = st.select(id="NV.NSMTC.G2.CHZ")
st = st.select(id="NV.NSMTC.G1.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.G1.CHZ | 2017-09-08T04:49:19.000000Z - 2017-09-08T05:49:19.000000Z | 100.0 Hz, 360001 samples
In [12]:
vars(st[0].stats.response.response_stages[0])
Out[12]:
{'_pz_transfer_function_type': 'LAPLACE (RADIANS/SECOND)',
 '_normalization_frequency': 10.0,
 'normalization_factor': 1.0,
 '_zeros': [0j, 0j],
 '_poles': [(-7.66549+9.95761j), (-7.66549-9.95761j)],
 'stage_sequence_number': 1,
 'input_units': 'M/S',
 'output_units': 'V',
 'input_units_description': 'Velocity in Meters per Second',
 'output_units_description': 'Volts',
 'resource_id': None,
 'resource_id2': None,
 'stage_gain': 78.74,
 'stage_gain_frequency': 10.0,
 'name': None,
 'description': None,
 'decimation_input_sample_rate': None,
 'decimation_factor': None,
 'decimation_offset': None,
 'decimation_delay': None,
 'decimation_correction': None}
In [13]:
# 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-78.74V/m/s fc=2.0 h=0.61"
st[0].stats.channel = "CHZ-78.74V/m/s fc=2.0 h=0.61: current IRIS metadata"

Correct instrument response for CN.PGC

In [14]:
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 G1.CHZ data

In [15]:
# select G1.CHZ data only 
st_mod2 = st_mod2.select(id="NV.NSMTC.G1.CHZ")
print(st_mod2)
st_mod3 = st_mod3.select(id="NV.NSMTC.G1.CHZ")
print(st_mod3)

1 Trace(s) in Stream:
NV.NSMTC.G1.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.G1.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 [16]:
#http://ds.iris.edu/NRL/sensors/oyo_geospace/RESP.XX.NS380..SLZ.HS1LT.3810.115000.2.76
# this from IRIS NRL   
# HS-1-LT 
G_mod2 = 76.2 # V/m/s (Open-Circuit Sensitivity)
fc_mod2 = 2.0 # corner frequency
damp_mod2 = 0.70 # damping
In [17]:
vars(st_mod2[0].stats.response.response_stages[0])
Out[17]:
{'_pz_transfer_function_type': 'LAPLACE (RADIANS/SECOND)',
 '_normalization_frequency': 10.0,
 'normalization_factor': 1.0,
 '_zeros': [0j, 0j],
 '_poles': [(-7.66549+9.95761j), (-7.66549-9.95761j)],
 'stage_sequence_number': 1,
 'input_units': 'M/S',
 'output_units': 'V',
 'input_units_description': 'Velocity in Meters per Second',
 'output_units_description': 'Volts',
 'resource_id': None,
 'resource_id2': None,
 'stage_gain': 78.74,
 'stage_gain_frequency': 10.0,
 'name': None,
 'description': None,
 'decimation_input_sample_rate': None,
 'decimation_factor': None,
 'decimation_offset': None,
 'decimation_delay': None,
 'decimation_correction': None}
In [18]:
# 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 =  [(-7.66549+9.95761j), (-7.66549-9.95761j)]
# updated pole & zero =  [(-8.79645943005142-8.974183634899006j), (-8.79645943005142+8.974183634899006j)]
# original gain value =  78.74
# updated  gain value =  -76.2
In [19]:
vars(st_mod2[0].stats.response.response_stages[0])
Out[19]:
{'_pz_transfer_function_type': 'LAPLACE (RADIANS/SECOND)',
 '_normalization_frequency': 10.0,
 'normalization_factor': 1.0,
 '_zeros': [0j, 0j],
 '_poles': [(-8.79645943005142-8.974183634899006j),
  (-8.79645943005142+8.974183634899006j)],
 'stage_sequence_number': 1,
 'input_units': 'M/S',
 'output_units': 'V',
 'input_units_description': 'Velocity in Meters per Second',
 'output_units_description': 'Volts',
 'resource_id': None,
 'resource_id2': None,
 'stage_gain': -76.2,
 'stage_gain_frequency': 10.0,
 'name': None,
 'description': None,
 'decimation_input_sample_rate': None,
 'decimation_factor': None,
 'decimation_offset': None,
 'decimation_delay': None,
 'decimation_correction': None}
In [20]:
paz_mod2 = {'poles': st_mod2[0].stats.response.response_stages[0]._poles, 'zeros': st_mod2[0].stats.response.response_stages[0]._zeros,'gain': 1, 'sensitivity':1 }
In [21]:
print(paz_mod2)

{'poles': [(-8.79645943005142-8.974183634899006j), (-8.79645943005142+8.974183634899006j)], 'zeros': [0j, 0j], 'gain': 1, 'sensitivity': 1}
In [22]:
from obspy.signal.invsim import paz_to_freq_resp
from obspy.signal.invsim import paz_2_amplitude_value_of_freq_resp
In [23]:
A0_mod2 = 1/(paz_2_amplitude_value_of_freq_resp(paz_mod2, st_mod2[0].stats.response.response_stages[0]._normalization_frequency))
print(A0_mod2)

1.0
In [ ]:

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

In [ ]:
In [24]:
# current IRIS metadaat with -1 polarity flip 
G_mod3 = 78.74 # V/m/s (Open-Circuit Sensitivity)
G_mod3 = 60.0
fc_mod3 = 2.0 # corner frequency
damp_mod3 = 0.61 # damping
In [25]:
# 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 =  [(-7.66549+9.95761j), (-7.66549-9.95761j)]
# updated pole & zero =  [(-7.665486074759095-9.957609836456946j), (-7.665486074759095+9.957609836456946j)]
# original gain value =  78.74
# updated  gain value =  -60.0

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

"mod2" uses a set of parameters from the spechseet

In [26]:
# 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 [27]:
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 [28]:
# updated channel name so that when we plot data we can see the sensor gain value info.
# location code is G2 -> C2 so that tihs treace would be plotted before G2 data. only plotting purpose
#st_mod2[0].stats.channel = "CHZ-"+str(G_mod2)+"V/m/s x -1 fc="+str(fc_mod2)+" h="+str(damp_mod2)
st_mod2[0].stats.channel = "CHZ-"+str(G_mod2)+"V/m/s x -1 fc="+str(fc_mod2)+" h="+str(damp_mod2)+": HS-1-LT IRIS NRL"

st_mod2[0].stats.location = "C2"
In [29]:
# updated channel name so that when we plot data we can see the sensor gain value info.
# location code is G2 -> C1so that tihs treace would be plotted before G2 data. only plotting purpose
st_mod3[0].stats.channel = "CHZ-"+str(G_mod3)+"V/m/s x -1 fc="+str(fc_mod3)+" h="+str(damp_mod3)
#st_mod3[0].stats.channel = "CHZ-"+str(G_mod3)+"V/m/s x -1 fc="+str(fc_mod3)+" h="+str(damp_mod3)+": current IRIS metadata x-1"
st_mod3[0].stats.location = "C1"

Filtering and plot

In [30]:
st_all = st2.copy() + st_mod2.copy() + st_mod3.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 = 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 [31]:
# plot each panel. also check max. amplitude value
for tr in st_all:
    print("# maximum amplitude = ",tr.max())
    #tr.plot(size=(800, 200))
    tr.plot( starttime=t1, endtime=t2, size=(800, 200))

# maximum amplitude =  0.000133466330181

# maximum amplitude =  0.000131618994453

# maximum amplitude =  0.000100577869963

# maximum amplitude =  0.000127609587696

# maximum amplitude =  -9.72387395639e-05
In [ ]:
In [ ]:
In [ ]:
In [ ]: