ObsPy Tutorial
Downloading/Processing Exercise
</div> </div> </div> image: User:Abbaszade656 / Wikimedia Commons / CC-BY-SA-4.0

Workshop for the "Training in Network Management Systems and Analytical Tools for Seismic"

Baku, October 2018

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 2014 South Napa 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.PFO 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 [ ]:
import obspy
from obspy.clients.fdsn import Client

c_event = Client("USGS")

# Event time.
event_time = obspy.UTCDateTime("2014-08-24T10:20:44.0")

# Get the event information. The temporal and magnitude constraints make it unique
cat = c_event.get_events(starttime=event_time - 10, endtime=event_time + 10,
                         minmagnitude=6)
print(cat)

c = Client("IRIS")
# Download station information at the response level!
inv = c.get_stations(network="II", station="PFO", location="00", channel="BHZ",
                     starttime=event_time - 10 * 60, endtime=event_time + 30 * 60,
                     level="response")
print(inv)

# Download 3 component waveforms.
st = c.get_waveforms(network="II", station="PFO", location="00",
                     channel="BHZ", starttime=event_time - 10 * 60,
                     endtime=event_time + 30 * 60)
print(st)

Have a look at the just downloaded data.

In [ ]:
inv.plot()
inv.plot_response(0.001)
cat.plot()
st.plot()

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 [ ]:
coords = inv.get_coordinates("II.PFO.00.BHZ")
coords

Step 2: Determine Coordinates of Event

In [ ]:
origin = cat[0].origins[0]
print(origin)

Step 3: Calculate distance of event and station.

Use obspy.geodetics.locations2degree.

In [ ]:
from obspy.geodetics import locations2degrees

distance = locations2degrees(origin.latitude, origin.longitude,
                             coords["latitude"], coords["longitude"])
print(distance)

Step 4: Calculate Theoretical Arrivals

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

m = TauPyModel(model="ak135")

arrivals = m.get_ray_paths(
    distance_in_degree=distance,
    source_depth_in_km=origin.depth / 1000.0)

arrivals.plot();

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

In [ ]:
first_arrival = origin.time + arrivals[0].time

print(first_arrival)

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

In [ ]:
st.trim(first_arrival - 60, first_arrival + 300)
st.plot();

Step 7: Remove the instrument response

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

taper

In [ ]:
st.remove_response(inventory=inv,
                   pre_filt=(1.0 / 100.0, 1.0 / 50.0, 10.0, 20.0),
                   output="VEL")
st.plot()

Bonus: Interactive IPython widgets

In [ ]:
from ipywidgets 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=(1, 180),
         depth=(0, 700), wavetype=["ttall", "diff"]);