wiki:ObspyTutorial
Last modified 8 months ago

This page has been moved to  http://tutorial.obspy.org!

ObsPy Tutorial (Old)

Table of Content


Python Introduction for Seismologists

Here we want to give a small, incomplete intro to the Python programming language, with links to useful packages and further resources. The key features are explained via the following Python script: 

#!/usr/bin/env python
import glob
from obspy.core import read

for file in glob.glob("*.z"):
    st = read(file)
    tr = st[0]
    msg = "%s %s %f %f" % (tr.stats.station, str(tr.stats.starttime), tr.data.mean(), tr.data.std())
    print msg
  • Line 1: Specifies the location of the python interpreter (not necessary on windows).
  • Lines 2-3: Import modules/libraries/packages in the current namespace. The glob module, which allows wildcard matching on filenames, is imported here. All functions inside this module can be accessed via glob.function() afterwards. Furthermore, a single function read() from the obspy.core module is imported, which is used to read various different seismogram file formats.
  • Line 4: Starts a for-loop using the glob function of the module glob, glob.glob(), on all files ending with z. Note: The length of all loops in Python is determined by the indentation level. (Do not mix spaces and tabs in your program code for indentation, this produces bugs that are not easy to identify).
  • Line 5: Uses the read function from the obspy.core module to read in the seismogram
  • Line 6: Assigns the first Trace object of the list-like Stream object to the variable tr
  • Line 7: Python has no printf-like function, but a Python counterpart for the sprintf function is the % operator acting on a format string. Here we print the header attributes station and the string representation of starttime as well as the return value of the methods mean() and std() acting on the data subobjects of the Trace (which are of type numpy.ndarray)
  • Line 8: Prints string of line 7

As Python is an interpreter language, we recommend to use a shell for rapid development and trying things out. We recommend the IPython shell, which supports tab completion, history expansion and various other features. E.g. type help(glob.glob) or glob.glob? to see the help of the glob() function (the module must be imported beforehand).

Further Resources:

Reading Seismograms

Seismograms of the formats SAC, MINISEED, GSE2 (CM6 compressed) and Q can be read into Python variables. The ObsPy Stream object resembles a data-only SEED (MiniSEED) file. It can contain several Traces, i.e. gap-less continuous sample chunks. Stream objects can be created via the read method.

The following example shows how a single GSE2-formatted seismogram is read into a Stream object. There is only one Trace in the seismogram: 

>>> from obspy.core import read
>>> st = read('http://examples.obspy.org/RJOB_061005_072159.ehz.new')
>>> print st
1 Trace(s) in Stream:
.RJOB..Z | 2005-10-06T07:21:59.849998Z - 2005-10-06T07:24:59.844998Z | 200.0 Hz, 36000 samples

Seismogram meta data are accessed via the stats keyword on each trace: 

>>> print st[0].stats
         network:
         station: RJOB
        location:
         channel: Z
       starttime: 2005-10-06T07:21:59.849998Z
         endtime: 2005-10-06T07:24:59.844998Z
   sampling_rate: 200.0
           delta: 0.005
            npts: 36000
           calib: 0.0948999971151
         _format: GSE2
            gse2: AttribDict({'instype': '      ', 'datatype': 'CM6', 'hang': -1.0, 'auxid': 'RJOB', 'vang': -1.0, 'calper': 1.0})
>>> st[0].stats.station
'RJOB'
>>> st[0].stats.gse2.datatype
'CM6'

The actual waveform data may be retrieved via the data keyword on each trace: 

>>> st[0].data
[-38  12  -4 ..., -14  -3  -9]
>>> st[0].data[0:3]
array([-38,  12,  -4])
>>> len(st[0])
36000

The Stream object offers a plotting method for fast previews (requires the obspy.imaging module): 

>>> st.plot(color='k')

See WaveformPlotting for more information about plotting seismograms.

Filtering Seismograms

The following script shows how to filter a seismogram. The example uses a zero-phase-shift low-pass filter with a corner frequency of 1 Hz using 2 corners. This is done in two runs forward and backward, so we end up with 4 corners de facto.

The available filters are:

  • "bandpass"
  • "bandstop"
  • "lowpass"
  • "highpass" 
import numpy as np
import matplotlib.pyplot as plt
from obspy.core import read

# Read the seismogram
st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")

# There is only one trace in the Stream object, let's work on that trace...
tr = st[0]

# Filtering with a lowpass on a copy of the original Trace
tr_filt = tr.copy()
tr_filt.filter('lowpass', freq=1.0, corners=2, zerophase=True)

# Now let's plot the raw and filtered data...
t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta)
plt.subplot(211)
plt.plot(t, tr.data, 'k')
plt.ylabel('Raw Data')
plt.subplot(212)
plt.plot(t, tr_filt.data, 'k')
plt.ylabel('Lowpassed Data')
plt.xlabel('Time [s]')
plt.suptitle(tr.stats.starttime)
plt.show()

Downsampling Seismograms

The following script shows how to downsample a seismogram. Currently, a simple integer decimation is supported. If not explicitely disabled, a low-pass filter is applied prior to decimation in order to prevent aliasing. For comparison, the non-decimated but filtered data is plotted as well. Applied processing steps are documented in trace.stats.processing of every single Trace. Note the shift that is introduced because by default the applied filters are not of zero-phase type. This can be avoided by manually applying a zero-phase filter and deactivating automatic filtering during downsampling (no_filter=True). 

import numpy as np
import matplotlib.pyplot as plt
from obspy.core import read

# Read the seismogram
st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")

# There is only one trace in the Stream object, let's work on that trace...
tr = st[0]

# Downsample the 200 Hz data by a factor of 4 to 50 Hz
# Note that this automatically includes a lowpass filtering with corner frequency 20 Hz
# We work on a copy of the original data just to demonstrate the effects of downsampling
tr_new = tr.copy()
tr_new.downsample(decimation_factor=4, strict_length=False)

# For comparison also only filter the original data (same filter options as in automatically
# applied filtering during downsampling, corner frequency 0.4 * new sampling rate)
tr_filt = tr.copy()
tr_filt.filter('lowpass', freq=0.4*tr.stats.sampling_rate/4.0)

# Now let's plot the raw and filtered data...
t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta)
t_new = np.arange(0, tr_new.stats.npts / tr_new.stats.sampling_rate, tr_new.stats.delta)

plt.plot(t, tr.data, 'k', label='Raw', alpha=0.3)
plt.plot(t, tr_filt.data, 'b', label='Lowpassed', alpha=0.7)
plt.plot(t_new, tr_new.data, 'r', label='Lowpassed/Downsampled', alpha=0.7)
plt.xlabel('Time [s]')
plt.xlim(82, 83.5)
plt.suptitle(tr.stats.starttime)
plt.legend()
plt.show()

Seismogram Envelopes

The following script shows how to filter a seismogram and plot it together with its envelope. This example uses a zero-phase-shift bandpass to filter the data with corner frequencies 1 and 3 Hz, using 2 corners (two runs due to zero-phase option, thus 4 corners overall). Then we calculate the envelope and plot it together with the Trace. Data can be found here Download. 

import numpy as np
import matplotlib.pyplot as plt
from obspy.core import read
import obspy.signal

st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
data = st[0].data
npts = st[0].stats.npts
samprate = st[0].stats.sampling_rate

# Filtering the Stream object
st_filt = st.copy()
st_filt.filter('bandpass', freqmin=1, freqmax=3, corners=2, zerophase=True)

# Envelope of filtered data
data_envelope = obspy.signal.filter.envelope(st_filt[0].data)

# The plotting, plain matplotlib
t = np.arange(0, npts/samprate, 1/samprate)
plt.plot(t, st_filt[0].data, 'k')
plt.plot(t, data_envelope, 'k:')
plt.title(st[0].stats.starttime)
plt.ylabel('Filtered Data w/ Envelope')
plt.xlabel('Time [s]')
plt.xlim(80, 90)
plt.show()

Seismometer Correction/Simulation

The following script shows how to simulate a 1Hz seismometer from a STS-2 seismometer with the given poles and zeros. Poles, zeros, gain (A0 normalization factor) and sensitivity (overall sensitivity) are specified as keys of a dictionary. 

from obspy.core import read
from obspy.signal import cornFreq2Paz

paz_sts2 = {'poles': [-0.037004+0.037016j, -0.037004-0.037016j,
                      -251.33+0j,
                      -131.04-467.29j, -131.04+467.29j],
            'zeros': [0j, 0j],
            'gain': 60077000.0,
            'sensitivity': 2516778400.0}
paz_1hz = cornFreq2Paz(1.0, damp=0.707) # 1Hz instrument
paz_1hz['sensitivity'] = 1.0

st = read()
# make a copy to keep our original data
st_orig = st.copy()

# Simulate instrument given poles, zeros and gain of
# the original and desired instrument
st.simulate(paz_remove=paz_sts2, paz_simulate=paz_1hz)

# plot original and simulated data
st_orig.plot()
st.plot()

For more customized plotting we could also work with matplotlib manually from here: 

import numpy as np
import matplotlib.pyplot as plt

tr = st[0]
tr_orig = st_orig[0]

t = np.arange(tr.stats.npts) / tr.stats.sampling_rate

plt.subplot(211)
plt.plot(t, tr_orig.data, 'k')
plt.ylabel('STS-2 [counts]')
plt.subplot(212)
plt.plot(t, tr.data, 'k')
plt.ylabel('1Hz Instrument [m/s]')
plt.xlabel('Time [s]')
plt.show()

Plotting Spectrogram

The following lines of code demonstrate how to make a spectrogram plot of an Obspy Stream object.

Lots of options can be customized, see  spectrogram(). For example, the colormap of the plot can easily be adjusted by importing a predefined colormap from matplotlib.cm, nice overviews of available matplotlib colormaps are given at:
 http://www.astro.lsa.umich.edu/~msshin/science/code/matplotlib_cm/
 http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps 

from obspy.core import read

st = read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
st.spectrogram(log=True)

Poles and Zeros, Frequency Response

The following lines show how to calculate and visualize the frequency response of a LE-3D/1s seismometer with sampling interval 0.005s and 16384 points of fft. Two things have to be taken into account for the phase (actually for the imaginary part of the response):

  • the fft that is used is defined as exp(-i*phi), but this minus sign is missing for the visualization, so we have to add it again
  • we want the phase to go from 0 to 2*pi, instead of the output from atan2 that goes from -pi to pi 
import numpy as np
import matplotlib.pyplot as plt
from obspy.signal import pazToFreqResp

paz = {}
paz['poles'] = [-4.440+4.440j, -4.440-4.440j, -1.083+0.0j]
paz['zeros'] = [0.0 +0.0j, 0.0 +0.0j, 0.0 +0.0j]
paz['gain'] = 0.4

h,f = pazToFreqResp(paz['poles'], paz['zeros'], paz['gain'], 0.005, 16384, freq=True)

plt.figure()
plt.subplot(121)
plt.loglog(f, abs(h))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')

plt.subplot(122)
phase = np.unwrap(np.arctan2(-h.imag, h.real)) #take negative of imaginary part
plt.semilogx(f, phase)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [radian]')
plt.suptitle('Frequency Response of LE-3D/1s Seismometer') #title, centered above both subplots
plt.subplots_adjust(wspace=0.3) #make more room inbetween subplots for the ylabel of right plot
plt.show()

The following lines show how to retrieve waveforms and poles and zeros from ArcLink. The poles and zeros are then used to correct for the instrument response and to simulate a 1Hz instrument with damping 0.707. Plots are provided. 

import numpy as np
import matplotlib.pyplot as plt
from obspy.core import UTCDateTime
from obspy.arclink import Client
from obspy.signal import cornFreq2Paz, seisSim

# Retrieve waveforms via ArcLink
client = Client(host="webdc.eu", port=18001)
t = UTCDateTime("2009-08-24 00:20:03")
st = client.getWaveform('BW', 'RJOB', '', 'EHZ', t, t+30)
paz = client.getPAZ('BW', 'RJOB', '', 'EHZ', t, t+30)
paz = paz.values()[0]

# 1Hz instrument
one_hertz = cornFreq2Paz(1.0)
# Correct for frequency response of the instrument
res = seisSim(st[0].data.astype('float32'), st[0].stats.sampling_rate, paz, inst_sim=one_hertz)
# Correct for overall sensitivity
res = res / paz['sensitivity']

# Plot the seismograms
sec = np.arange(len(res)) / st[0].stats.sampling_rate
plt.subplot(211)
plt.plot(sec,st[0].data, 'k')
plt.title("%s %s" % (st[0].stats.station, t))
plt.ylabel('STS-2')
plt.subplot(212)
plt.plot(sec, res, 'k')
plt.xlabel('Time [s]')
plt.ylabel('1Hz CornerFrequency')
plt.show()

Beachball plot

The following lines show how to create a graphical representation of a focal mechanism. 

from obspy.imaging.beachball import Beachball

mt = [0.91, -0.89, -0.02, 1.78, -1.55, 0.47]
Beachball(mt, size=200, linewidth=2, facecolor='b')

Basemap plot with beachballs

The following example shows how to plot beachballs into a basemap plot together with some stations. The example requires the  basemap package ( download site) to be installed, the srtm file used can be downloaded  here. 

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from obspy.imaging.beachball import Beach
import gzip

# read in topo data (on a regular lat/lon grid)
# (srtm data from: http://srtm.csi.cgiar.org/)
srtm = np.loadtxt(gzip.open("srtm_1240-1300E_4740-4750N.asc.gz"), skiprows=8)

# origin of data grid as stated in srtm data file header
# create arrays with all lon/lat values from min to max and
lats = np.linspace(47.8333, 47.6666, srtm.shape[0])
lons = np.linspace(12.6666, 13.0000, srtm.shape[1])

# create Basemap instance with Mercator projection
# we want a slightly smaller region than covered by our srtm data
m = Basemap(projection='merc', lon_0=13, lat_0=48, resolution="h",
            llcrnrlon=12.75, llcrnrlat=47.69, urcrnrlon=12.95, urcrnrlat=47.81)

# create grids and compute map projection coordinates for lon/lat grid
x, y = m(*np.meshgrid(lons,lats))

# Make contour plot
cs = m.contour(x, y, srtm, 40, colors="k", lw=0.5, alpha=0.3)
m.drawcountries(color="red", linewidth=1)

# Draw a lon/lat grid (20 lines for an interval of one degree)
m.drawparallels(np.linspace(47, 48, 21), labels=[1,1,0,0], fmt="%.2f", dashes=[2,2])
m.drawmeridians(np.linspace(12, 13, 21), labels=[0,0,1,1], fmt="%.2f", dashes=[2,2])

# Plot station positions and names into the map
# again we have to compute the projection of our lon/lat values
lats = [47.761659, 47.7405, 47.755100, 47.737167]
lons = [12.864466, 12.8671, 12.849660, 12.795714]
names = [" RMOA", " RNON", " RTSH", " RJOB"]
x, y = m(lons, lats)
m.scatter(x, y, 200, color="r", marker="v", edgecolor="k", zorder=3)
for i in range(len(names)):
    plt.text(x[i], y[i], names[i], va="top", family="monospace", weight="bold")

# Add beachballs for two events
lats = [47.751602, 47.75577]
lons = [12.866492, 12.893850]
x, y = m(lons, lats)
# Two focal mechanisms for the beachball routine, specified as [strike, dip, rake]
focmecs = [[80, 50, 80], [85, 30, 90]]
ax = plt.gca()
for i in range(len(focmecs)):
    b = Beach(focmecs[i], xy=(x[i], y[i]), width=1000, linewidth=1)
    b.set_zorder(10)
    ax.add_collection(b)

plt.show()

The first lines of our SRTM data file (from  CGIAR) look like this: 

ncols         400
nrows         200
xllcorner     12°40'E
yllcorner     47°40'N
xurcorner     13°00'E
yurcorner     47°50'N
cellsize      0.00083333333333333
NODATA_value  -9999
682 681 685 690 691 689 678 670 675 680 681 679 675 671 674 680 679 679 675 671 668 664 659 660 656 655 662 666 660 659 659 658 ....

Some notes:

  • The Python package  GDAL allows you to directly read a GeoTiff into numpy array (geo = gdal.Open("file.geotiff"); x = geo.ReadAsArray())
  • GeoTiff elevation data is available e.g. from  ASTER
  • Shading/Illumination? can be added. See the basemap example  plotmap_shaded.py for more info.

Anything to MiniSEED

The following lines show how you can convert anything to MiniSEED format. In the example, a few lines of a weather station output are written to a MiniSEED file. The correct meta information starttime, the sampling_rate, station name and so forth are also encoded (Note: Only the ones given are allowed by the MiniSEED standard). Converting arbitrary ASCII to MiniSEED is extremely helpful if you want to send log messages, output of meteorologic stations or anything else via the seedlink protocol. 

import numpy as np
from obspy.core import read, Trace, Stream, UTCDateTime

weather = """
00.0000 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000000
00.0002 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000001
00.0005 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000002
00.0008 0.0 ??? 4.7 97.7 1015.4 0.0 010308 000003
00.0011 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000004
00.0013 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000005
00.0016 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000006
00.0019 0.0 ??? 4.7 97.7 1015.0 0.0 010308 000007
"""

# Convert to numpy character array
data = np.fromstring(weather, dtype='|S1')

# Fill header attributes
stats = {'network': 'BW', 'station': 'RJOB', 'location': '',
         'channel': 'WLZ', 'npts': len(data), 'sampling_rate': 0.1,
         'mseed' : {'dataquality' : 'D'}}
stats['starttime'] = UTCDateTime() #the current time
st = Stream([Trace(data=data, header=stats)])
st.write("weather.mseed", format='MSEED', encoding=0, reclen=256) #encoding 0==ASCII

# Show that it worked, convert numpy character array back to string
st1 = read("weather.mseed")
print st1[0].data.tostring()

Checking Data availability for data files

Often, you have a bunch of data and want to know which station is available at what time. For this purpose, ObsPy ships the  obspy-scan script (automatically available after installation), which automatically detects the format (MiniSEED, SAC, SACXY, GSE2, SH-ASC, SH-Q, SEISAN) from the header of the data files. Gaps are plotted as vertical red lines, start times of available data are plotted as crosses --- the data itself are plotted as horizontal lines.

The script can be used to scan through 1000s of files (already used with 30000 files, execution time ca. 45min), month/year ranges are plotted automatically. It opens an interactive plot in which you can zoom in...

Execute something like following line from the command prompt (NOT from inside python or ipython), use e.g. wildcards to match the files: 

obspy-scan /bay_mobil/mobil/20090622/1081019/*_1.*

Saving Python objects as MATLAB MAT files

The following example shows how to read in a file (MSEED, SAC, SEISAN, GSE2, SH) with Python and save each Trace in the resulting Stream object to one MATLAB MAT file. The data can the be loaded from within MATLAB with the load function. 

from obspy.core import read
from scipy.io import savemat

st = read("test/BW.BGLD..EHE.D.2008.001")
for i, tr in enumerate(st):
    mdict = dict([[j, str(k)] for j,k in tr.stats.iteritems()])
    mdict['data'] = tr.data
    savemat("data-%d.mat" % i, mdict)

Interfacing R-cran from Python

The  rpy2 package allows to interface R from Python. The following example shows how to convert data (numpy.ndarray) to an R matrix and execute the R command summary on it. 

>>> from obspy.core import read
>>> import rpy2.robjects as RO
>>> import rpy2.robjects.numpy2ri
>>> r = RO.r
>>> 
>>> st = read("test/BW.BGLD..EHE.D.2008.001")
>>> M = RO.RMatrix(st[0].data)
>>> print r.summary(M)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-1056.0  -409.0  -393.0  -393.7  -378.0   233.0

Seismogram to ASCII

In the following, a small Python script is shown which converts each Trace of a seismogram file (supported formats: MiniSEED, SAC, Seisan, GSE2, SH_ASC, Q) to an ASCII file and multiplies the data by a calibration factor. 

"""
USAGE: seis2asc.py mseed_file out_file calib_factor
"""
from obspy.core import read
import numpy as np
import sys

try:
    in_file = sys.argv[1]
    of_file = sys.argv[2]
    calib = float(sys.argv[3])
except:
    print __doc__
    raise

st = read(in_file)
for i, tr in enumerate(st):
    f = open("%s_%d" % (of_file, i), "w");
    f.write("# STATION %s\n" % (tr.stats.station))
    f.write("# CHANNEL %s\n" % (tr.stats.channel))
    f.write("# START_TIME %s\n" % (str(tr.stats.starttime)))
    f.write("# SAMP_FREQ %f\n" % (tr.stats.sampling_rate))
    f.write("# NDAT %d\n" % (tr.stats.npts))
    np.savetxt(f, tr.data * calib, fmt="%f")
    f.close()

Merging Seismograms

The following example shows how to merge and plot three seismograms with overlaps, the longest one is taken to be the right one. Please also refer to the documentation of the  merge function. 

#!/usr/bin/env python

from obspy.core import read
import numpy as np, matplotlib.pyplot as plt

# Read in all files starting with dis.
st  = read("http://examples.obspy.org/dis.G.SCZ.__.BHE")
st += read("http://examples.obspy.org/dis.G.SCZ.__.BHE.1")
st += read("http://examples.obspy.org/dis.G.SCZ.__.BHE.2")

# Go through the stream object, determine time range in julian seconds
# and plot the data with a shared x axis
ax = plt.subplot(4,1,1) #dummy for tying axis
for i in range(3):
    plt.subplot(4,1,i+1, sharex=ax)
    t = np.linspace(st[i].stats.starttime.timestamp, 
                    st[i].stats.endtime.timestamp,
                    st[i].stats.npts)
    plt.plot(t, st[i].data)

# Merge the data together and show plot in a similar way
st.merge(method=1)
plt.subplot(4,1,4, sharex=ax)
t = np.linspace(st[0].stats.starttime.timestamp, 
                st[0].stats.endtime.timestamp,
                st[0].stats.npts)
plt.plot(t, st[0].data, 'r')
plt.show()

Beamforming --- FK Analysis

The following code shows how to do an FK Analysis with ObsPy. The data are from the blasting of the AGFA skyscraper in Munich (pickle.dumped stream object with paz attached to each Trace, can also be found  here). We execute sonic using the following settings: The slowness grid is set to corner values of -3.0 to 3.0 s/km with a step fraction of sl_s = 0.03. The window length is 1.0 s, using a step fraction of 0.05s. The data is bandpass filtered, using corners at 1.0 and 8.0 Hz, prewhitening is disabled. semb_thres and vel_thres are set to infinitesimally small numbers and must not be changed. The timestamp will be written in mlabhours, which can be read directly by our plotting routine. stime and etime have to be given in the UTCDateTime format. The output will be stored in 'out'. 

import pickle, urllib
from obspy.core import UTCDateTime
from obspy.signal.array_analysis import sonic
from obspy.signal import cornFreq2Paz

#
# Load data
st = pickle.load(urllib.urlopen("http://examples.obspy.org/agfa.dump"))

#
# Instrument correction to 1Hz corner frequency
paz1hz = cornFreq2Paz(1.0, damp=0.707)
st.simulate(paz_remove='self', paz_simulate=paz1hz)

#
# Execute sonic
kwargs = dict(
    # slowness grid: X min, X max, Y min, Y max, Slow Step
    sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03,
    # sliding window propertieds
    win_len=1.0, win_frac=0.05,
    # frequency properties
    frqlow=1.0, frqhigh=8.0, prewhiten=0,
    # restrict output
    semb_thres=-1e9, vel_thres=-1e9, verbose=True, timestamp='mlabhour',
    stime=UTCDateTime("20080217110515"), etime=UTCDateTime("20080217110545")
)
out = sonic(st, **kwargs)

The output can be plotted by the following code lines. We read the output produced by sonic ('out'), which is: timestamp, rel. power, abs. power, backazimut and slowness. The colorbar corresponds to rel. power. 

import matplotlib.pyplot as plt

labels = 'rel.power abs.power baz slow'.split()

fig = plt.figure()
for i, lab in enumerate(labels):
    ax = fig.add_subplot(4,1,i+1)
    ax.scatter(out[:,0], out[:,i+1], c=out[:,1], alpha=0.6, edgecolors='none')
    ax.set_ylabel(lab)
    ax.xaxis_date()

fig.autofmt_xdate()
fig.subplots_adjust(top=0.95, right=0.95, bottom=0.2, hspace=0)

or by a polar plot. The plot sums the relative power in gridded bins, each defined by backazimuth and slowness of the analysed signal part. The backazimuth is counted clockwise from north, the slowness limits can be set by hand. 

import numpy as np
import matplotlib.cm as cm
import matplotlib.pyplot as plt
from matplotlib.colorbar import ColorbarBase
from matplotlib.colors import Normalize

cmap = cm.hot_r
pi = np.pi

#
# make output human readable, adjust backazimuth to values between 0 and 360
t, rel_power, abs_power, baz, slow = out.T
baz[baz < 0.0] += 360

# choose number of fractions in plot (desirably 360°/N is an integer!) 
N = 30
abins = np.arange(N+1)*360./N 
sbins = np.linspace(0, 3, N+1)

# sum rel power in bins given by abins and sbins
hist, baz_edges, sl_edges = np.histogram2d(baz, slow,
        bins=[abins, sbins], weights=rel_power)

# transfrom to gradient
baz_edges = baz_edges/180*np.pi

# add polar and colorbar axes
fig = plt.figure(figsize=(8,8))
cax = fig.add_axes([0.85, 0.2, 0.05, 0.5])
ax  = fig.add_axes([0.10, 0.1, 0.70, 0.7], polar=True)

dh = abs(sl_edges[1] - sl_edges[0])
dw = abs(baz_edges[1] - baz_edges[0])
# circle through backazimut
for i, row in enumerate(hist):
    bars = ax.bar(left=(pi/2 - (i+1)*dw)*np.ones(N),
                  height=dh*np.ones(N),
                  width=dw, bottom=dh*np.arange(N),
                  color=cmap(row/hist.max()))

ax.set_xticks([pi/2, 0, 3./2*pi, pi])
ax.set_xticklabels(['N', 'E' , 'S', 'W'])
# set slowness limits
ax.set_ylim(0, 3)
ColorbarBase(cax, cmap=cmap,
        norm=Normalize(vmin=hist.min(), vmax=hist.max()))

plt.show()

better polar plot as requested by Tobi

Array Response Function

The following code block shows how to plot the array transfer function for beam forming as a function of wavenumber using the ObsPy function obspy.signal.array_analysis.array_transff_wavenumber(). 

import matplotlib.pyplot as plt 
import numpy as np
from obspy.signal.array_analysis import array_transff_wavenumber

# generate array coordinates
coords = np.array([[  10.,   60., 0.],
                   [ 200.,   50., 0.],
                   [-120.,  170., 0.],
                   [-100., -150., 0.],
                   [  30., -220., 0.]])

# coordinates in km
coords /= 1000.

# set limits for wavenumber differences to analyze
klim = 40. 
kxmin = -klim
kxmax = klim
kymin = -klim
kymax = klim
kstep = klim/100.

# compute transfer function as a function of wavenumber difference
transff = array_transff_wavenumber(coords, klim, kstep, coordsys='xy')

# plot
plt.pcolor(np.arange(kxmin, kxmax + kstep * 1.1, kstep) - kstep / 2., 
           np.arange(kymin, kymax + kstep * 1.1, kstep) - kstep / 2., transff.T)

plt.colorbar()
plt.clim(vmin=0., vmax=1.)
plt.xlim(kxmin, kxmax)
plt.ylim(kymin, kymax)
plt.show()

Correcting arraytransff.png, thanks to Martin van Driel

Continuous Wavelet Transform

Small script doing the continuous wavelet transform using the  mlpy package for infrasound data recorded at Yasur in 2008. Further details on wavelets can be found e.g  http://en.wikipedia.org/wiki/Morlet_wavelet (in the wikipedia article the omega0 factor is denoted as sigma) [really sloppy and possibly incorrect: the omega0 factor tells you how often the wavelet fits into the time window, dj defines the spacing in the scale domain] 

#!/usr/bin/env python 

import matplotlib.pyplot as plt
from obspy.core import read
import numpy as np
import mlpy

tr = read("http://examples.obspy.org/a02i.2008.240.mseed")[0]

omega0 = 8
spec, scale = mlpy.cwt(tr.data, dt=tr.stats.delta, dj=0.05, wf='morlet', p=omega0, 
                       extmethod='none', extlength='powerof2')
# approximate scales through frequencies
freq = (omega0 + np.sqrt(2.0 + omega0**2))/(4*np.pi * scale[1:])

fig = plt.figure()
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2])
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60])
ax3 = fig.add_axes([0.83,0.1,0.03,0.6])
t = np.arange(tr.stats.npts)/tr.stats.sampling_rate
ax1.plot(t, tr.data, 'k')
img = ax2.imshow(np.abs(spec), extent=[t[0], t[-1], freq[-1], freq[0]], aspect='auto') 
fig.colorbar(img, cax=ax3)
plt.show()

Coordinate Conversions

Coordinate conversions can be done conveniently using  pyproj. After looking up the  EPSG codes of source and target coordinate system, the conversion can be done in just a few lines of code. The following example converts the station coordinates of two German stations to the regionally used Gauß-Krüger system: 

import pyproj
lat = [49.6919, 48.1629]
lon = [11.2217, 11.2752]
proj_wgs84 = pyproj.Proj(init="epsg:4326")
proj_gk4 = pyproj.Proj(init="epsg:31468")
x, y = pyproj.transform(proj_wgs84, proj_gk4, lon, lat)

>>> print x, y
[4443947.179, 4446185.667] [5506428.401, 5336354.055]

Hierarchical Clustering

An implementation of hierarchical clustering is provided in the package  hcluster. Among other things, it allows to build clusters from similarity matrices (e.g. computed using the cross-correlation routines in obspy.signal) and make dendrogram plots.
The following example shows how to do this for an already computed similarity matrix. The similarity data are computed from events in an area with induced seismicity in southern Germany and can be fetched from ObsPy's examples page:: 

import pickle, urllib
import matplotlib.pyplot as plt
import hcluster

dissimilarity = pickle.load(urllib.urlopen("http://examples.obspy.org/dissimilarities.pkl"))

plt.subplot(121)
plt.imshow(1 - dissimilarity, interpolation="nearest")

dissimilarity = hcluster.squareform(dissimilarity)
threshold = 0.3
linkage = hcluster.linkage(dissimilarity, method="single")
clusters = hcluster.fcluster(linkage, 0.3, criterion="distance")
plt.xlabel("Event number")
plt.ylabel("Event number")

plt.subplot(122)
hcluster.dendrogram(linkage, color_threshold=0.3)
plt.xlabel("Event number")
plt.ylabel("Dissimilarity")
plt.show()

Building a new dataless-SEED file

A small script how to clone an existing dataless-SEED file (dataless.seed.BW_RNON) and use this as a template to build up a dataless-SEED file for a new station. [Note that the FIR coefficients are not set, in case you require correct FIR coefficients, clone from an existing dataless file with the same seismometer type] 

from urllib import urlopen
from obspy.core import UTCDateTime
from obspy.xseed import Parser

url = "http://examples.obspy.org/dataless.seed.BW_RNON"
p = Parser(urlopen(url).read())
blk = p.blockettes

blk[50][0].network_code = 'BW'
blk[50][0].station_call_letters = 'RMOA'
blk[50][0].site_name = "Moar Alm, Bavaria, BW-Net"
blk[50][0].latitude = 47.761658
blk[50][0].longitude = 12.864466
blk[50][0].elevation = 815.0
blk[50][0].start_effective_date = UTCDateTime("2006-07-18T00:00:00.000000Z")
blk[50][0].end_effective_date = ""
blk[33][1].abbreviation_description = "Lennartz LE-3D/1 seismometer"

mult = len(blk[58])/3
for i, cha in enumerate(['Z', 'N', 'E']):
    blk[52][i].channel_identifier = 'EH%s' % cha
    blk[52][i].location_identifier = ''
    blk[52][i].latitude = blk[50][0].latitude
    blk[52][i].longitude = blk[50][0].longitude
    blk[52][i].elevation = blk[50][0].elevation
    blk[52][i].start_date = blk[50][0].start_effective_date
    blk[52][i].end_date = blk[50][0].end_effective_date
    blk[53][i].number_of_complex_poles = 3
    blk[53][i].real_pole = [-4.444, -4.444, -1.083]
    blk[53][i].imaginary_pole = [+4.444, -4.444, +0.0]
    blk[53][i].real_pole_error = [0, 0, 0]
    blk[53][i].imaginary_pole_error = [0, 0, 0]
    blk[53][i].number_of_complex_zeros = 3
    blk[53][i].real_zero = [0.0, 0.0, 0.0]
    blk[53][i].imaginary_zero = [0.0, 0.0, 0.0]
    blk[53][i].real_zero_error = [0, 0, 0]
    blk[53][i].imaginary_zero_error = [0, 0, 0]
    blk[53][i].A0_normalization_factor = 1.0
    blk[53][i].normalization_frequency = 3.0
    # stage sequence number 1, seismometer gain
    blk[58][i*mult].sensitivity_gain = 400.0
    # stage sequence number 2, digitizer gain
    blk[58][i*mult+1].sensitivity_gain = 1677850.0
    # stage sequence number 0, overall sensitivity
    blk[58][(i+1)*mult-1].sensitivity_gain = 671140000.0

p.writeSEED("dataless.seed.BW_RMOA")

Visualizing Probabilistic Power Spectral Densities

The following code example shows how to use the PPSD class defined in  obspy.signal. The routine is useful for interpretation of e.g. noise measurements for site quality control checks. For more information on the topic see  McNamara and Buland 2004. 

from obspy.core import read
from obspy.xseed import Parser
from obspy.signal import PPSD

st = read("http://examples.obspy.org/BW.KW1..EHZ.D.2011.037")
# select a trace with the desired station/channel combination
tr = st.select(id="BW.KW1..EHZ")[0]
# acquire poles and zeros information, e.g. from a dataless SEED file
parser = Parser("http://examples.obspy.org/dataless.seed.BW_KW1")
paz = parser.getPAZ(tr.id)
# initialize a new PPSD instance. the ppsd object will then make sure that
# only appropriate data go into the probabilistic psd statistics.
ppsd = PPSD(tr.stats, paz)

# now we can add data (either trace or stream objects) to the ppsd estimate
# this step may take a while..
ppsd.add(st)

# we can check what time ranges are represented in the ppsd estimate..
# ppsd.times contains a sorted list of starttimes of the one hour long slices
# that the psds are computed from
print ppsd.times
# the count of 1 hour long segments in the statistics can be checked via..
print "number of psd segments: ", len(ppsd.times)

# adding the same stream again will do nothing, the ppsd object makes sure
# that no overlapping data segments go into the ppsd estimate
ppsd.add(st)
print "number of psd segments: ", len(ppsd.times)

# additional information from other files/sources
# can be added step by step (like client.getWaveform(...))
st = read("http://examples.obspy.org/BW.KW1..EHZ.D.2011.038")
ppsd.add(st)

# the graphical representation of the ppsd can be saved to a file like..
ppsd.plot("/tmp/ppsd.png")
ppsd.plot("/tmp/ppsd.pdf")
# ... or displayed in a matplotlib window
ppsd.plot()

Below the actual PPSD (for a detailed discussion see  McNamara and Buland 2004) is a visualization of the data basis for the PPSD (can also be switched off during plotting). The top row shows data fed into the PPSD, green patches represent available data, red patches represent gaps in streams that were added to the PPSD. The bottom row in blue shows the single psd measurements that go into the histogram. The default processing method fills gaps with zeros, these data segments then show up as single outlying psd lines.

Note: Providing metadata from e.g. a Dataless SEED volume is safer than specifying static poles and zeros information (see  API)

PPSD example

Attachments