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

# 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

Making SmartSolo instrument response

from IRIS NRL

In [3]:
import obspy
from obspy import read_inventory
from obspy.core.inventory import Inventory, Network, Station, Channel, Site
from obspy.clients.nrl import NRL

Gain 0dB response file

output file is station_SmartSolo_0dB.xml

In [4]:
# lat, lon ele from BK.BKS
st_lat = 37.87622
st_lon = -122.23558
st_ele = 243.9
st_dep = 25.6

# We'll first create all the various objects. These strongly follow the
# hierarchy of StationXML files.
inv = Inventory(
    # We'll add networks later.
    networks=[],
    # The source should be the id whoever create the file.
    source="ObsPy-Tutorial")

net = Network(
    # This is the network code according to the SEED standard.
    code="AA",
    # A list of stations. We'll add one later.
    stations=[],
    description="SmartSolo test",
    # Start-and end dates are optional.
    start_date=obspy.UTCDateTime(2019, 1, 1))

sta1 = Station(
    # This is the station code according to the SEED standard.
    code="B1568",
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    creation_date=obspy.UTCDateTime(2019, 1, 1),
    site=Site(name="Inside the Byerly Vault "))

sta2 = Station(
    # This is the station code according to the SEED standard.
    code="B1482",
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    creation_date=obspy.UTCDateTime(2019, 1, 1),
    site=Site(name="Outside the Byerly Vault "))


cha_z = Channel(
    # This is the channel code according to the SEED standard.
    code="DPZ",
    # This is the location code according to the SEED standard.
    location_code="00",
    # Note that these coordinates can differ from the station coordinates.
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    depth=st_dep,
    azimuth=0.0,
    dip=-90.0,
    sample_rate=500)

cha_n = Channel(
    # This is the channel code according to the SEED standard.
    code="DPN",
    # This is the location code according to the SEED standard.
    location_code="00",
    # Note that these coordinates can differ from the station coordinates.
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    depth=st_dep,
    azimuth=0.0,
    dip=0.0,
    sample_rate=500)
    
    
cha_e = Channel(
    # This is the channel code according to the SEED standard.
    code="DPE",
    # This is the location code according to the SEED standard.
    location_code="00",
    # Note that these coordinates can differ from the station coordinates.
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    depth=st_dep,
    azimuth=90.0,
    dip=0.0,
    sample_rate=500)
    
# By default this accesses the NRL online. Offline copies of the NRL can
# also be used instead
nrl = NRL()
# The contents of the NRL can be explored interactively in a Python prompt,
# see API documentation of NRL submodule:
# http://docs.obspy.org/packages/obspy.clients.nrl.html
# Here we assume that the end point of data logger and sensor are already
# known:

#print(nrl.sensors['DTCC (manuafacturers of SmartSolo)']['5 Hz']['Rc=1850, Rs=430000'])
#print(nrl.dataloggers['DTCC (manufacturers of SmartSolo']['SmartSolo IGU-16']['0 dB (1)']['500']['Minimum Phase']['Off'])
response = nrl.get_response( # doctest: +SKIP
    sensor_keys=['DTCC (manuafacturers of SmartSolo)', '5 Hz', 'Rc=1850, Rs=430000'],
    datalogger_keys=['DTCC (manufacturers of SmartSolo', 'SmartSolo IGU-16', '0 dB (1)', '500', 'Minimum Phase', 'Off'])
    #datalogger_keys=['DTCC (manufacturers of SmartSolo', 'SmartSolo IGU-16', '0 dB (1)', '500', 'Linear Phase', 'Off'])


# Now tie it all together.
cha_z.response = response
cha_n.response = response
cha_e.response = response

sta1.channels.append(cha_z)
sta1.channels.append(cha_n)
sta1.channels.append(cha_e)


sta2.channels.append(cha_z)
sta2.channels.append(cha_n)
sta2.channels.append(cha_e)

net.stations.append(sta1)
net.stations.append(sta2)

inv.networks.append(net)

# And finally write it to a StationXML file. We also force a validation against
# the StationXML schema to ensure it produces a valid StationXML file.
#
# Note that it is also possible to serialize to any of the other inventory
# output formats ObsPy supports.
inv.write("station_SmartSolo_0dB.xml", format="stationxml", validate=True)

Gain 36dB response file

output file is station_SmartSolo_36dB.xml

In [5]:
# lat, lon ele from BK.BKS
st_lat = 37.87622
st_lon = -122.23558
st_ele = 243.9
st_dep = 25.6

# We'll first create all the various objects. These strongly follow the
# hierarchy of StationXML files.
inv = Inventory(
    # We'll add networks later.
    networks=[],
    # The source should be the id whoever create the file.
    source="ObsPy-Tutorial")

net = Network(
    # This is the network code according to the SEED standard.
    code="AA",
    # A list of stations. We'll add one later.
    stations=[],
    description="SmartSolo test",
    # Start-and end dates are optional.
    start_date=obspy.UTCDateTime(2019, 1, 1))

sta1 = Station(
    # This is the station code according to the SEED standard.
    code="B1568",
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    creation_date=obspy.UTCDateTime(2019, 1, 1),
    site=Site(name="Inside the Byerly Vault "))

sta2 = Station(
    # This is the station code according to the SEED standard.
    code="B1482",
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    creation_date=obspy.UTCDateTime(2019, 1, 1),
    site=Site(name="Outside the Byerly Vault "))


cha_z = Channel(
    # This is the channel code according to the SEED standard.
    code="DPZ",
    # This is the location code according to the SEED standard.
    location_code="00",
    # Note that these coordinates can differ from the station coordinates.
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    depth=st_dep,
    azimuth=0.0,
    dip=-90.0,
    sample_rate=500)

cha_n = Channel(
    # This is the channel code according to the SEED standard.
    code="DPN",
    # This is the location code according to the SEED standard.
    location_code="00",
    # Note that these coordinates can differ from the station coordinates.
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    depth=st_dep,
    azimuth=0.0,
    dip=0.0,
    sample_rate=500)
    
    
cha_e = Channel(
    # This is the channel code according to the SEED standard.
    code="DPE",
    # This is the location code according to the SEED standard.
    location_code="00",
    # Note that these coordinates can differ from the station coordinates.
    latitude=st_lat,
    longitude=st_lon,
    elevation=st_ele,
    depth=st_dep,
    azimuth=90.0,
    dip=0.0,
    sample_rate=500)
    
# By default this accesses the NRL online. Offline copies of the NRL can
# also be used instead
nrl = NRL()
# The contents of the NRL can be explored interactively in a Python prompt,
# see API documentation of NRL submodule:
# http://docs.obspy.org/packages/obspy.clients.nrl.html
# Here we assume that the end point of data logger and sensor are already
# known:

#print(nrl.sensors['DTCC (manuafacturers of SmartSolo)']['5 Hz']['Rc=1850, Rs=430000'])
#print(nrl.dataloggers['DTCC (manufacturers of SmartSolo']['SmartSolo IGU-16']['36 dB (64)']['500']['Minimum Phase']['Off'])
response = nrl.get_response( # doctest: +SKIP
    sensor_keys=['DTCC (manuafacturers of SmartSolo)', '5 Hz', 'Rc=1850, Rs=430000'],
    datalogger_keys=['DTCC (manufacturers of SmartSolo', 'SmartSolo IGU-16', '36 dB (64)', '500', 'Minimum Phase', 'Off'])


# Now tie it all together.
cha_z.response = response
cha_n.response = response
cha_e.response = response

sta1.channels.append(cha_z)
sta1.channels.append(cha_n)
sta1.channels.append(cha_e)


sta2.channels.append(cha_z)
sta2.channels.append(cha_n)
sta2.channels.append(cha_e)

net.stations.append(sta1)
net.stations.append(sta2)

inv.networks.append(net)

# And finally write it to a StationXML file. We also force a validation against
# the StationXML schema to ensure it produces a valid StationXML file.
#
# Note that it is also possible to serialize to any of the other inventory
# output formats ObsPy supports.
inv.write("station_SmartSolo_36dB.xml", format="stationxml", validate=True)
In [6]:
inv_0dB =read_inventory("station_SmartSolo_0dB.xml")
print(inv_0dB)

Inventory created at 2020-08-25T12:58:31.896188Z
	Created by: ObsPy 1.2.2
		    https://www.obspy.org
	Sending institution: ObsPy-Tutorial
	Contains:
		Networks (1):
			AA
		Stations (2):
			AA.B1482 (Outside the Byerly Vault )
			AA.B1568 (Inside the Byerly Vault )
		Channels (6):
			AA.B1482.00.DPZ, AA.B1482.00.DPN, AA.B1482.00.DPE, AA.B1568.00.DPZ
			AA.B1568.00.DPN, AA.B1568.00.DPE
In [7]:
inv_36dB =read_inventory("station_SmartSolo_36dB.xml")
print(inv_36dB)

Inventory created at 2020-08-25T12:58:33.095761Z
	Created by: ObsPy 1.2.2
		    https://www.obspy.org
	Sending institution: ObsPy-Tutorial
	Contains:
		Networks (1):
			AA
		Stations (2):
			AA.B1482 (Outside the Byerly Vault )
			AA.B1568 (Inside the Byerly Vault )
		Channels (6):
			AA.B1482.00.DPZ, AA.B1482.00.DPN, AA.B1482.00.DPE, AA.B1568.00.DPZ
			AA.B1568.00.DPN, AA.B1568.00.DPE

Set client (Data Center)

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

In [8]:
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 [9]:
# use Z, N or E
channel = "Z"
#channel = "N"
#channel = "E"
In [10]:
# BKS HNZ (EpiSensor) data as reference
sta = "BKS" # station
net = "BK" # network
#com = "HN"+channel # component
com = "HN"+channel+",HH"+channel # component # read borh broadband (HHZ) and accelerometer (HNZ) data
loc = "00" #  location
In [11]:
# separately added Nodal seismic data 
sta2 = "B1568" # station
net2 = "AA" # network # not used
com2 = "DP"+channel # component 
loc2 = "00" # location # not used
#client2 = Client("IRIS") # data from IRIS

sta3 = "B1482" # station
net3 = "AA" # network # not used
com3 = "DP"+channel # component 
loc3 = "00" # location # not used

Set time window

This example uses 30-sec data from a local M2.5 event

In [12]:
#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 Broadband (HHZ) and EpiSensor (HNZ) data

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

In [13]:
# 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 [14]:
# st2 is  SmartSolo 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+"/"+sta2+"."+com2+"."+mseed_time+"T00:00:00.000.mseed")

# add B1482
st2 += read(dirTMP+"/"+sta3+"."+com3+"."+mseed_time+"T00:00:00.000.mseed")

# trim
st2.trim(starttime, endtime) 
_plot = st2.plot()

Copy st2 data

one is for 0dB response file and the other is for 36dB response file

In [15]:
st2_0dB = st2.copy() # use it woth 0 dB respomse file
st2_36dB = st2.copy() # use if with 36 dB response file

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 with "0 dB" response file

use "decimate" (500sps -> 100sps) and "simulate" from obspy

In [17]:
st2_0dB.detrend() # remove liner trend
st2_0dB.detrend("demean") # demean
st2_0dB.taper(0.05) # cosin taper
# decimate to 100Hz
st2_0dB.decimate(factor=5, strict_length=False)

st2_0dB = st2_0dB.remove_response( output="VEL", inventory=inv_0dB) # 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_0dB.plot()

Add gain value informaion to channel name

In [18]:
st2_0dB[0].stats.channel = st2_36dB[0].stats.channel +" Gain 0 dB"
st2_0dB[1].stats.channel = st2_36dB[1].stats.channel +" Gain 0 dB"

Correct instrument response for SmartSolo seismic data with "36 dB" response file

use "decimate" (500sps -> 100sps) and "simulate" from obspy

In [19]:
st2_36dB.detrend() # remove liner trend
st2_36dB.detrend("demean") # demean
st2_36dB.taper(0.05) # cosin taper
# decimate to 100Hz
st2_36dB.decimate(factor=5, strict_length=False)

st2_36dB = st2_36dB.remove_response( output="VEL", inventory=inv_36dB) # 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_36dB.plot()

Add gain value informaion to channel name

Station name uses "A" instead of B so that these two waveforms are plotted as top two traces

In [20]:
st2_36dB[0].stats.station = "A1482" 
st2_36dB[1].stats.station = "A1568" 
st2_36dB[0].stats.channel = st2_36dB[0].stats.channel +" Gain 36 dB"
st2_36dB[1].stats.channel = st2_36dB[1].stats.channel +" Gain 36 dB"

Filtering and plot

Example uses 1-15 Hz passband

In [21]:
st_all = st2_0dB.copy() + st2_36dB.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 = .5 # in Hz 
fh = 2 # in Hz
fl = .5 # in Hz 
fh = 5 # in Hz

fl = 5
fh = 10

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")
t2 = UTCDateTime("2019-10-15T05:24:20.00")

#_plot = st_all.plot(size=(800, 400))
_plot = st_all.plot(size=(900, 800), starttime=t1, endtime=t2)
In [22]:
for tr in st_all:
    print("# max. amplitude = ",tr.max())
    tr.plot(starttime=t1, endtime=t2)
    #tr.plot(size=(1200,300))
    

# max. amplitude =  -1.35200760412e-05

# max. amplitude =  2.03463752386e-05

# max. amplitude =  -2.11251188144e-07

# max. amplitude =  3.17912113103e-07

# max. amplitude =  -1.28428735189e-05

# max. amplitude =  -1.3117474133e-05
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: