ObsPy Tutorial
Downloading/Processing Exercise

Seismo-Live: http://seismo-live.org

Authors:

For the this exercise we will download some data from the Tohoku-Oki earthquake, cut out a certain time window around the first arrival and remove the instrument response from the data.

In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8

The first step is to download all the necessary information using the ObsPy FDSN client. Learn to read the documentation!

We need the following things:

  1. Event information about the Tohoku-Oki earthquake. Use the get_events() method of the client. A good provider of event data is the USGS.
  2. Waveform information for a certain station. Choose your favorite one! If you have no preference, use II.BFO which is available for example from IRIS. Use the get_waveforms() method.
  3. Download the associated station/instrument information with the get_stations() method.
In [ ]:
In [ ]:

Have a look at the just downloaded data.

In [ ]:
In [ ]:

Exercise

The goal of this exercise is to cut the data from 1 minute before the first arrival to 5 minutes after it, and then remove the instrument response.

Step 1: Determine Coordinates of Station

In [ ]:
In [ ]:

Step 2: Determine Coordinates of Event

In [ ]:
In [ ]:

Step 3: Calculate distance of event and station.

Use obspy.geodetics.locations2degree.

In [ ]:
In [ ]:

Step 4: Calculate Theoretical Arrivals

from obspy.taup import TauPyModel
m = TauPyModel(model="ak135")
arrivals = m.get_ray_paths(...)
arrivals.plot()
In [ ]:
In [ ]:

Step 5: Calculate absolute time of the first arrivals at the station

In [ ]:
In [ ]:

Step 6: Cut to 1 minute before and 5 minutes after the first arrival

In [ ]:
In [ ]:

Step 7: Remove the instrument response

st.remove_response(inventory=inv, pre_filt=...)

taper

In [ ]:
In [ ]:

Bonus: Interactive IPython widgets

In [ ]:
from IPython.html.widgets import interact
from obspy.taup import TauPyModel

m = TauPyModel("ak135")

def plot_raypaths(distance, depth, wavetype):
    try:
        plt.close()
    except:
        pass
    if wavetype == "ttall":
        phases = ["ttall"]
    elif wavetype == "diff":
        phases = ["Pdiff", "pPdiff"]
    m.get_ray_paths(distance_in_degree=distance,
                    source_depth_in_km=depth,
                    phase_list=phases).plot();
    
interact(plot_raypaths, distance=[0, 180],
         depth=[0, 700], wavetype=["ttall", "diff"]);