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

Set client (Data Center)

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

In [2]:
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)? This example uses BRK.BK.HHN.00 data

In [3]:
# BRK BHZ data (https://seismo.berkeley.edu/station_book/brk.html; Haviland Hall)
sta = "BRK" # station
com = "HHZ" # component 
com = "HHN" # component 

net = "BK" # network
loc = "00" # location "--" for blank location code

Set time window

This example uses 1-min data for the 2020 M3.1 Alum Rock earthquake

In [4]:
# Alum Rock event
#M 3.1 - 8km NE of Alum Rock, CA
#2020-08-02 12:40:24 (UTC)37.409°N 121.755°W7.9 km depth
start_day = "2020-08-02T12:40:24"
end_day = "2020-08-02T12:41:24"
starttime = UTCDateTime(start_day)
endtime = UTCDateTime(end_day)

Download BRK seismic data

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

In [5]:
st = client.get_waveforms(network=net, station=sta, location=loc, channel=com,
                     starttime=starttime, endtime=endtime, 
                     attach_response=True)
_plot = st.plot()

Correct instrument response

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

In [6]:
st.detrend() # remove liner trend
st.detrend("demean") # demean
st.taper(0.05) # cosin tape
Out[6]:
1 Trace(s) in Stream:
BK.BRK.00.HHN | 2020-08-02T12:40:24.008393Z - 2020-08-02T12:41:23.998393Z | 100.0 Hz, 6000 samples
In [7]:
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)
In [8]:
_plot = st.plot()

Filtering seismic data

first remove liner trend, apply a cosin taper, and then do filtering

In [9]:
st.detrend() # remove liner trend
st.detrend("demean") # demean
st.taper(0.05) # cosin tape

fl = 5 # in Hz 
fh = 10 # in Hz
st.filter(type='bandpass', freqmin=fl, freqmax=fh, corners=6, zerophase=False)
Out[9]:
1 Trace(s) in Stream:
BK.BRK.00.HHN | 2020-08-02T12:40:24.008393Z - 2020-08-02T12:41:23.998393Z | 100.0 Hz, 6000 samples
In [10]:
_plot = st.plot()

Reading DAS data from SEGY format

use _read_regy from obspy

read 4 segy files for the Alum Rock event

Each file has 1885 traces

st_das will have traces from these 4 files. Total traces would be 7540 (1885 * 4)

In [11]:
from obspy.io.segy.core import _read_segy
In [12]:
filename  = "../segy_20200802_124027.876+0000.segy"
filename2 = "../segy_20200802_124037.876+0000.segy"
filename3 = "../segy_20200802_124047.876+0000.segy"
filename4 = "../segy_20200802_124057.876+0000.segy"

st_das = _read_segy(filename)
st_das += _read_segy(filename2)
st_das += _read_segy(filename3)
st_das += _read_segy(filename4)

print(st_das)

7540 Trace(s) in Stream:

Seq. No. in line:    1 | 2020-08-02T12:40:27.000000Z - 2020-08-02T12:40:36.999500Z | 2000.0 Hz, 20000 samples
...
(7538 other traces)
...
Seq. No. in line: 1885 | 2020-08-02T12:40:57.000000Z - 2020-08-02T12:41:06.999500Z | 2000.0 Hz, 20000 samples

[Use "print(Stream.__str__(extended=True))" to print all Traces]

Select Trace

Example uses trace number "1200"

from the 1st segy file -> traces[1200]

from the 2nd segy file -> traces[1200+1885]

st_das_select would be "one trace" from chennal numner "1200"

In [13]:
tr_n = 1200
#tr_n = 1000
#tr_n = 1800
#tr_n = 1100
#tr_n = 800
#tr_n = 500
#tr_n = 1600

das_tr_num = 1885
tr_n2 = tr_n + das_tr_num
tr_n3 = tr_n2 + das_tr_num
tr_n4 = tr_n3 + das_tr_num

print(tr_n)
print(tr_n2)
print(tr_n3)
print(tr_n4)

# adding traces
# it seems that obspy automatically merge all traces becuase the same line number(?)
st_das_select  = st_das[tr_n]
st_das_select += st_das[tr_n2]
st_das_select += st_das[tr_n3]
st_das_select += st_das[tr_n4]

#print(st_das[tr_n])
#print(st_das[tr_n2])
# check 
print(st_das_select)

1200
3085
4970
6855
Seq. No. in line: 1201 | 2020-08-02T12:40:27.000000Z - 2020-08-02T12:41:06.999500Z | 2000.0 Hz, 80000 samples

Add header info.

This is optional. not need for the data analyses below

station code = RFS

network code = DS

component(channel) code = DPX

location code = 1200 (<- DAS trace)

In [14]:
st_das_select.stats['station'] = "RFS"
st_das_select.stats['network'] = "DS"
st_das_select.stats['channel'] = "DPX"
st_das_select.stats['location'] = str(tr_n)

Plot

Checking all data are merged without gaps

In [15]:
_plot = st_das_select.plot()

BRK seismic amplitude x 100000

to mach DAS raw data in counts. This "1000000 factor" for 5-25 Hz data. If you use differnt passbands, please modify it.

In [16]:
st2 = st.copy()
st2[0].data = st2[0].data * 100000

Filtering and plot

Example uses 5-20 Hz bandpass filter

In [17]:
st_all = st2.copy() + st_das_select.copy()

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

# select passband
#st_all.filter("bandpass", freqmin=5, freqmax=25, corners=6, zerophase=True) 
#st_all.filter("bandpass", freqmin=0.1, freqmax=5, corners=6, zerophase=True)
st_all.filter("bandpass", freqmin=5, freqmax=20, corners=6, zerophase=True)
#st_all.filter("bandpass", freqmin=20, freqmax=100, corners=6, zerophase=True)

_plot = st_all.plot(size=(1000, 400))
#_plot = st_all.plot(starttime=starttime, endtime=endtime)


#for tr in st_all:
    #tr.plot(starttime=starttime, endtime=endtime)
    #tr.plot(size=(800, 200))

Spectrogram

In [18]:
st_das_select2 = st_das_select.copy()
st_das_select2.detrend() # remove liner trend
st_das_select2.detrend("demean") # demean
st_das_select2.taper(0.05) # cosin taper

#st_das_select2.filter("bandpass", freqmin=0.1, freqmax=800, corners=6, zerophase=True)
#st_das_select2.filter("bandpass", freqmin=1, freqmax=10, corners=6, zerophase=True)

st_das_select2.spectrogram(log=True, title='DAS trace ' + str(tr_n))
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: