Instaseis Tutorial
Part 3: Plot a Record Section



Advanced Exercise 1: Plot Record Section

Use Instaseis to calculate a record section of your choice.



  1. ObsPy Trace objects provide a function times() to conveniently get the time axis for plotting and the normalize() function to normalize seismograms to a common amplitude.
  2. Use high pass filtering and deep sources if you want to enhance body waves
  3. np.linspace() can help to generate equidistant stations
  4. obspy.taup.TauPyModel can be used to compare with ray theoretical arrivals

Basic lines to set up the notebook and some paths.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os
import obspy'ggplot')
plt.rcParams['figure.figsize'] = (10, 8)

Import Instaseis and open the database:

In [2]:
import instaseis
db = instaseis.open_db("data/database")


In [3]:
from obspy.taup import TauPyModel
from collections import defaultdict

m = TauPyModel(model="ak135")

# some paramters
depth_in_km = 150.
mindist = 0.
maxdist = 180.
numrec = 30
fmin = 0.02
fmax = 0.1
component = "Z"
phases = ["P", "PP", "Pdiff", "S", "SS", "PS"]
colors = ["r", "r", "r", "b", "b", "g"]

# define instaseis source
src = instaseis.Source.from_strike_dip_rake(
    latitude=0., longitude=0.,
    depth_in_m=depth_in_km * 1e3,
    M0=1e+21, strike=32., dip=62., rake=90.)

# storage for traveltimes
distances = defaultdict(list)
ttimes = defaultdict(list)

# loop over distances
for dist in np.linspace(mindist, maxdist, numrec):
    # define receiver
    rec = instaseis.Receiver(latitude=0, longitude=dist)

    # generate seismogram, filter and plot
    tr = db.get_seismograms(source=src, receiver=rec, components=[component])[0]
    tr.filter('highpass', freq=fmin)
    tr.filter('lowpass', freq=fmax)
    plt.plot(tr.times(), * 5 + dist, color="black")
    # get traveltimes
    arrivals = m.get_travel_times(distance_in_degree=dist,
    for arr in arrivals:

# plot traveltimes
for color, phase in zip(colors, phases):
    plt.scatter(ttimes[phase], distances[phase], s=40, color=color)

plt.xlim(tr.times()[0], tr.times()[-1])
plt.ylim(mindist - 3, maxdist + 3)
plt.xlabel('time / s')
plt.ylabel('epicentral distance / degree')
/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/obspy/taup/ UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.