%matplotlib inline
from __future__ import print_function
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
Stream
object.from obspy import read
st = read("./data/waveform_PFO.mseed", format="mseed")
print(st)
2 Trace(s) in Stream: II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.10.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.994500Z | 40.0 Hz, 120000 samples
UNIX wildcards can be used to read multiple files simultaneously
automatic file format detection, no need to worry about file formats
from obspy import read
st = read("./data/waveform_*")
print(st)
8 Trace(s) in Stream: GR.BFO..BHE | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples GR.BFO..BHN | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples GR.BFO..BHZ | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.10.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.994500Z | 40.0 Hz, 120000 samples SY.PFO.S3.MXE | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXN | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXZ | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples
from obspy import UTCDateTime
t = UTCDateTime("2011-03-11T05:46:23.015400Z")
st = read("./data/waveform_*", starttime=t + 10 * 60, endtime=t + 12 * 60)
print(st)
8 Trace(s) in Stream: GR.BFO..BHE | 2011-03-11T05:56:23.021088Z - 2011-03-11T05:58:23.021088Z | 20.0 Hz, 2401 samples GR.BFO..BHN | 2011-03-11T05:56:23.021088Z - 2011-03-11T05:58:23.021088Z | 20.0 Hz, 2401 samples GR.BFO..BHZ | 2011-03-11T05:56:23.021088Z - 2011-03-11T05:58:23.021088Z | 20.0 Hz, 2401 samples II.PFO.00.BHZ | 2011-03-11T05:56:23.019500Z - 2011-03-11T05:58:23.019500Z | 20.0 Hz, 2401 samples II.PFO.10.BHZ | 2011-03-11T05:56:23.019500Z - 2011-03-11T05:58:23.019500Z | 40.0 Hz, 4801 samples SY.PFO.S3.MXE | 2011-03-11T05:56:23.084250Z - 2011-03-11T05:58:23.078750Z | 6.2 Hz, 744 samples SY.PFO.S3.MXN | 2011-03-11T05:56:23.084250Z - 2011-03-11T05:58:23.078750Z | 6.2 Hz, 744 samples SY.PFO.S3.MXZ | 2011-03-11T05:56:23.084250Z - 2011-03-11T05:58:23.078750Z | 6.2 Hz, 744 samples
from obspy import read
st = read("./data/waveform_PFO.mseed")
print(type(st))
<class 'obspy.core.stream.Stream'>
print(st)
2 Trace(s) in Stream: II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.10.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.994500Z | 40.0 Hz, 120000 samples
print(st.traces)
[<obspy.core.trace.Trace object at 0x10c4925c0>, <obspy.core.trace.Trace object at 0x10c492be0>]
print(st[0])
II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples
+
operator (or using the .append()
and .extend()
methods)st1 = read("./data/waveform_PFO.mseed")
st2 = read("./data/waveform_PFO_synthetics.mseed")
st = st1 + st2
print(st)
st3 = read("./data/waveform_BFO_BHE.sac")
st += st3
print(st)
5 Trace(s) in Stream: II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.10.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.994500Z | 40.0 Hz, 120000 samples SY.PFO.S3.MXE | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXN | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXZ | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples 6 Trace(s) in Stream: II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.10.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.994500Z | 40.0 Hz, 120000 samples SY.PFO.S3.MXE | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXN | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXZ | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples GR.BFO..BHE | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples
for tr in st:
print(tr.id)
II.PFO.00.BHZ II.PFO.10.BHZ SY.PFO.S3.MXE SY.PFO.S3.MXN SY.PFO.S3.MXZ GR.BFO..BHE
Trace.stats
) that fully describes the time series by specifying..st = read("./data/waveform_PFO.mseed")
tr = st[0] # get the first Trace in the Stream
print(tr)
II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples
print(tr.stats)
network: II station: PFO location: 00 channel: BHZ starttime: 2011-03-11T05:46:23.019500Z endtime: 2011-03-11T06:36:22.969500Z sampling_rate: 20.0 delta: 0.05 npts: 60000 calib: 1.0 _format: MSEED mseed: AttribDict({'dataquality': 'M', 'encoding': 'STEIM1', 'filesize': 356352, 'record_length': 4096, 'byteorder': '>', 'number_of_records': 87})
print(tr.stats.delta, "|", tr.stats.endtime)
0.05 | 2011-03-11T06:36:22.969500Z
tr.stats.sampling_rate = 5.0
print(tr.stats.delta, "|", tr.stats.endtime)
0.2 | 2011-03-11T09:06:22.819500Z
print(tr.stats.npts)
60000
tr.data = tr.data[:100]
print(tr.stats.npts, "|", tr.stats.endtime)
100 | 2011-03-11T05:46:42.819500Z
tr = read("./data/waveform_PFO.mseed")[0]
tr.plot()
print(tr)
tr.resample(sampling_rate=100.0)
print(tr)
II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:23.009500Z | 100.0 Hz, 300000 samples
print(tr)
tr.trim(tr.stats.starttime + 12 * 60, tr.stats.starttime + 14 * 60)
print(tr)
tr.plot()
II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:23.009500Z | 100.0 Hz, 300000 samples II.PFO.00.BHZ | 2011-03-11T05:58:23.019500Z - 2011-03-11T06:00:23.019500Z | 100.0 Hz, 12001 samples
tr.detrend("linear")
tr.taper(max_percentage=0.05, type='cosine')
tr.filter("lowpass", freq=0.1)
tr.plot()
# try tr.<Tab> for other methods defined for Trace
tr.detrend?
numpy.ndarray
(as Trace.data
)print(tr.data[:20])
[ 0.00000000e+00 1.14214131e-11 1.36532184e-10 8.27107696e-10 3.42897869e-09 1.10577664e-08 2.98034791e-08 7.02499050e-08 1.49294808e-07 2.92257039e-07 5.35256796e-07 9.27855448e-07 1.53594149e-06 2.44484947e-06 3.76269890e-06 5.62394072e-06 8.19309900e-06 1.16686963e-05 1.62873519e-05 2.23280418e-05]
..by doing arithmetic operations (fast, handled in C by NumPy)
print(tr.data ** 2 + 0.5)
[ 5.00000000e-01 5.00000000e-01 5.00000000e-01 ..., 2.19703659e+10 2.19251720e+10 2.18793418e+10]
..by using numpy.ndarray
builtin methods (also done in C by NumPy)
print(tr.data.max())
print(tr.data.mean())
print(tr.data.ptp())
# try tr.data.<Tab> for a list of numpy methods defined on ndarray
560511.615943 6903.56237784 1074087.25528
..by using numpy
functions (also done in C by NumPy)
import numpy as np
print(np.abs(tr.data))
# you can try np.<Tab> but there is a lot in there
# try np.a<Tab>
[ 0.00000000e+00 1.14214131e-11 1.36532184e-10 ..., 1.48224039e+05 1.48071510e+05 1.47916672e+05]
..by feeding pointers to existing C/Fortran routines from inside Python!
This is done internally in several places, e.g. for cross correlations, beamforming or in third-party filetype libraries like e.g. libmseed.
numpy.ndarray
(e.g. when needing to parse waveforms from non-standard ascii files)from obspy import Trace
x = np.random.randint(-100, 100, 500)
tr = Trace(data=x)
tr.stats.station = "XYZ"
tr.stats.starttime = UTCDateTime()
tr.plot()
from obspy import Stream
tr2 = Trace(data=np.random.randint(-300, 100, 1000))
tr2.stats.starttime = UTCDateTime()
tr2.stats.sampling_rate = 10.0
st = Stream([tr, tr2])
st.plot()
Stream
/ Trace
¶st.<Tab>
).st.filter()
- Filter all attached traces.st.trim()
- Cut all traces.st.resample()
/ st.decimate()
- Change the sampling rate.st.trigger()
- Run triggering algorithms.st.plot()
/ st.spectrogram()
- Visualize the data.st.attach_response()
/st.remove_response()
, st.simulate()
- Instrument correctionst.merge()
, st.normalize()
, st.detrend()
, st.taper()
, ...Stream
object can also be exported to many formats, so ObsPy can be used to convert between different file formats.st = read("./data/waveform_*.sac")
st.write("output_file.mseed", format="MSEED")
!ls -l output_file*
-rw-r--r-- 1 lion staff 737280 Sep 14 22:51 output_file.mseed
numpy.ndarray
with zeros and (e.g. use numpy.zeros()
) and put an ideal pulse somewhere in itTrace
object with your data arrayx = np.zeros(300)
x[100] = 1.0
tr = Trace(data=x)
tr.stats.station = "ABC"
tr.stats.sampling_rate = 20.0
tr.stats.starttime = UTCDateTime(2014, 2, 24, 15, 0, 0)
print(tr)
tr.plot()
.ABC.. | 2014-02-24T15:00:00.000000Z - 2014-02-24T15:00:14.950000Z | 20.0 Hz, 300 samples
tr.filter(...)
and apply a lowpass filter with a corner frequency of 1 Hertz.tr.filter("lowpass", freq=1)
tr.plot()
tr.trim(...)
to remove some of the zeros at start and at the endtr.trim(tr.stats.starttime + 3, tr.stats.endtime - 5)
tr.plot()
np.random.randn()
)tr.data = tr.data * 500
tr.data = tr.data + np.random.randn(len(tr))
tr.plot()
st = read("./data/waveform_*")
print(st)
8 Trace(s) in Stream: GR.BFO..BHE | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples GR.BFO..BHN | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples GR.BFO..BHZ | 2011-03-11T05:46:23.021088Z - 2011-03-11T06:36:22.971088Z | 20.0 Hz, 60000 samples II.PFO.00.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.969500Z | 20.0 Hz, 60000 samples II.PFO.10.BHZ | 2011-03-11T05:46:23.019500Z - 2011-03-11T06:36:22.994500Z | 40.0 Hz, 120000 samples SY.PFO.S3.MXE | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXN | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples SY.PFO.S3.MXZ | 2011-03-11T05:47:31.587750Z - 2011-03-11T06:36:22.974250Z | 6.2 Hz, 18152 samples
st.select()
to only keep traces of station BFO in the stream. Show the preview plot.st = st.select(station="BFO")
st.plot()
wlen=50
for the spectrogram plot)t1 = UTCDateTime(2011, 3, 11, 5, 55)
st.trim(t1, t1 + 10 * 60)
st.plot()
st.spectrogram(log=True, wlen=50);
st.detrend("linear")
st.filter("lowpass", freq=0.1)
st.plot()