ObsPy Tutorial
Handling Waveform Data
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
- read waveform data is returned as a
Streamobject.
In [2]:
from obspy import read
st = read("./data/waveform_PFO.mseed", format="mseed")
print(st)
- UNIX wildcards can be used to read multiple files simultaneously
automatic file format detection, no need to worry about file formats
- currently supported: mseed, sac, segy, seg2, gse1/2, seisan, sh, datamark, css, wav, y, Q (keeps growing...)
- more file formats are included whenever a basic reading routine is provided (or e.g. sufficient documentation on data compression etc.)
- custom user-specific file formats can be added (through plugin) to filetype autodiscovery in local ObsPy installation by user
In [3]:
from obspy import read
st = read("./data/waveform_*")
print(st)
- for MiniSEED files, only reading short portions of files without all of the file getting read into memory is supported (saves time and memory when working on large collections of big files)
In [4]:
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)
The Stream Object¶
- A Stream object is a collection of Trace objects
In [5]:
from obspy import read
st = read("./data/waveform_PFO.mseed")
print(type(st))
In [6]:
print(st)
In [7]:
print(st.traces)
In [8]:
print(st[0])
- More traces can be assembled using
+operator (or using the.append()and.extend()methods)
In [9]:
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)
- convenient (and nicely readable) looping over traces
In [10]:
for tr in st:
print(tr.id)
- Stream is useful for applying the same processing to a larger number of different waveforms or to group Traces for processing (e.g. three components of one station in one Stream)
The Trace Object¶
- a Trace object is a single, contiguous waveform data block (i.e. regularly spaced time series, no gaps)
- a Trace object contains a limited amount of metadata in a dictionary-like object (as
Trace.stats) that fully describes the time series by specifying..- recording location/instrument (network, station, location and channel code)
- start time
- sampling rate
In [11]:
st = read("./data/waveform_PFO.mseed")
tr = st[0] # get the first Trace in the Stream
print(tr)
In [12]:
print(tr.stats)
- For custom applications it is sometimes necessary to directly manipulate the metadata of a Trace.
- The metadata of the Trace will stay consistent, as all values are derived from the starttime, the data and the sampling rate and are updated automatically
In [13]:
print(tr.stats.delta, "|", tr.stats.endtime)
In [14]:
tr.stats.sampling_rate = 5.0
print(tr.stats.delta, "|", tr.stats.endtime)
In [15]:
print(tr.stats.npts)
In [16]:
tr.data = tr.data[:100]
print(tr.stats.npts, "|", tr.stats.endtime)
- convenience methods make basic manipulations simple
In [17]:
tr = read("./data/waveform_PFO.mseed")[0]
tr.plot()
Out[17]:
In [18]:
print(tr)
tr.resample(sampling_rate=100.0)
print(tr)
In [19]:
print(tr)
tr.trim(tr.stats.starttime + 12 * 60, tr.stats.starttime + 14 * 60)
print(tr)
tr.plot()
Out[19]:
In [20]:
tr.detrend("linear")
tr.taper(max_percentage=0.05, type='cosine')
tr.filter("lowpass", freq=0.1)
tr.plot()
Out[20]:
In [21]:
# try tr.<Tab> for other methods defined for Trace
tr.detrend?
- Raw data available as a
numpy.ndarray(asTrace.data)
In [22]:
print(tr.data[:20])
- Data can be directly modified e.g. ..
..by doing arithmetic operations (fast, handled in C by NumPy)
In [23]:
print(tr.data ** 2 + 0.5)
..by using numpy.ndarray builtin methods (also done in C by NumPy)
In [24]:
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
..by using numpy functions (also done in C by NumPy)
In [25]:
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>
..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.
- Trace objects can also be manually generated with data in a
numpy.ndarray(e.g. when needing to parse waveforms from non-standard ascii files)
In [26]:
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()
Out[26]:
- Stream objects can be assembled from manually generated Traces
In [27]:
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()
Out[27]:
Builtin methods defined on Stream / Trace¶
- Most methods that work on a Trace object also work on a Stream object. They are simply executed for every trace. See ObsPy documentation for an overview of available methods (or try
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(), ...
- A
Streamobject can also be exported to many formats, so ObsPy can be used to convert between different file formats.
In [28]:
st = read("./data/waveform_*.sac")
st.write("output_file.mseed", format="MSEED")
!ls -l output_file*
Trace Exercises¶
- Make an
numpy.ndarraywith zeros and (e.g. usenumpy.zeros()) and put an ideal pulse somewhere in it - initialize a
Traceobject with your data array - Fill in some station information (e.g. network, station, ..)
- Print trace summary and plot the trace
- Change the sampling rate to 20 Hz
- Change the starttime of the trace to the start time of this session
- Print the trace summary and plot the trace again
In [ ]:
In [29]:
x = 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()
Out[29]:
- Use
tr.filter(...)and apply a lowpass filter with a corner frequency of 1 Hertz. - Display the preview plot, there are a few seconds of zeros that we can cut off.
In [ ]:
In [30]:
tr.filter("lowpass", freq=1)
tr.plot()
Out[30]:
- Use
tr.trim(...)to remove some of the zeros at start and at the end - show the preview plot again
In [ ]:
In [31]:
tr.trim(tr.stats.starttime + 3, tr.stats.endtime - 5)
tr.plot()
Out[31]:
- Scale up the amplitudes of the trace by a factor of 500
- Add standard normal gaussian noise to the trace (use
np.random.randn()) - Display the preview plot again
In [ ]:
In [32]:
tr.data = tr.data * 500
tr.data = tr.data + np.random.randn(len(tr))
tr.plot()
Out[32]:
Stream Exercises¶
- Read all Tohoku example earthquake data into a stream object ("./data/waveform_*")
- Print the stream summary
In [ ]:
In [33]:
st = read("./data/waveform_*")
print(st)
- Use
st.select()to only keep traces of station BFO in the stream. Show the preview plot.
In [ ]:
In [34]:
st = st.select(station="BFO")
st.plot()
Out[34]:
- trim the data to a 10 minute time window around the first arrival (just roughly looking at the preview plot)
- display the preview plot and spectrograms for the stream (with logarithmic frequency scale, use
wlen=50for the spectrogram plot)
In [ ]:
In [35]:
t1 = UTCDateTime(2011, 3, 11, 5, 55)
st.trim(t1, t1 + 10 * 60)
st.plot()
st.spectrogram(log=True, wlen=50);
- remove the linear trend from the data, apply a tapering and a lowpass at 0.1 Hertz
- show the preview plot again
In [ ]:
In [36]:
st.detrend("linear")
st.filter("lowpass", freq=0.1)
st.plot()
Out[36]: