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 [1]:
%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 [2]:
import obspy
from obspy.clients.fdsn import Client

c_event = Client("USGS")

# Event time.
event_time = obspy.UTCDateTime("2011-03-11T05:46:23.2")

# 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=9)
print(cat)

c = Client("IRIS")
# Download station information at the response level!
inv = c.get_stations(network="II", station="BFO", location="*", channel="BH?",
                     starttime=event_time - 60, endtime=event_time + 3600,
                     level="response")
print(inv)

# Download 3 component waveforms.
st = c.get_waveforms(network="II", station="BFO", location="*",
                     channel="BH?", starttime=event_time - 60,
                     endtime=event_time + 3600)
print(st)
1 Event(s) in Catalog:
2011-03-11T05:46:24.120000Z | +38.297, +142.373 | 9.1 mww | manual
Inventory created at 2019-11-07T13:52:58.000000Z
	Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.37
		    http://service.iris.edu/fdsnws/station/1/query?starttime=2011-03-11...
	Sending institution: IRIS-DMC (IRIS-DMC)
	Contains:
		Networks (1):
			II
		Stations (1):
			II.BFO (Black Forest Observatory, Schiltach, Germany)
		Channels (3):
			II.BFO.00.BHZ, II.BFO.00.BHN, II.BFO.00.BHE
3 Trace(s) in Stream:
II.BFO.00.BHE | 2011-03-11T05:45:23.215600Z - 2011-03-11T06:46:23.165251Z | 20.0 Hz, 73200 samples
II.BFO.00.BHN | 2011-03-11T05:45:23.215100Z - 2011-03-11T06:46:23.164751Z | 20.0 Hz, 73200 samples
II.BFO.00.BHZ | 2011-03-11T05:45:23.215300Z - 2011-03-11T06:46:23.164951Z | 20.0 Hz, 73200 samples

Have a look at the just downloaded data.

In [ ]:
 
In [3]:
inv.plot()
inv.plot_response(0.001)
cat.plot()
st.plot()
/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/obspy/imaging/maps.py:45: UserWarning: basemap/pyproj with proj4 version >= 5 has a bug that results in inverted map axes. Your maps may be wrong. Please use another version of proj4, or use cartopy.
  warnings.warn(msg)
Out[3]:

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 [4]:
coords = inv.get_coordinates("II.BFO.00.BHE")
coords
Out[4]:
{'latitude': 48.330101,
 'longitude': 8.3296,
 'elevation': 589.0,
 'local_depth': 0.0}

Step 2: Determine Coordinates of Event

In [ ]:
 
In [5]:
origin = cat[0].preferred_origin()
print(origin)
Origin
	        resource_id: ResourceIdentifier(id="quakeml:earthquake.usgs.gov/archive/product/origin/official20110311054624120_30/official/1561150942463/product.xml")
	               time: UTCDateTime(2011, 3, 11, 5, 46, 24, 120000)
	          longitude: 142.373
	           latitude: 38.297
	              depth: 29000.0 [uncertainty=0.0]
	            quality: OriginQuality(used_station_count=541, standard_error=1.16, azimuthal_gap=9.5)
	 origin_uncertainty: OriginUncertainty(horizontal_uncertainty=0.0, preferred_description='horizontal uncertainty')
	    evaluation_mode: 'manual'
	      creation_info: CreationInfo(agency_id='US', creation_time=UTCDateTime(2019, 6, 21, 21, 2, 22, 463000))

Step 3: Calculate distance of event and station.

Use obspy.geodetics.locations2degree.

In [ ]:
 
In [6]:
from obspy.geodetics import locations2degrees

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

Step 4: Calculate Theoretical Arrivals

from obspy.taup import TauPyModel
m = TauPyModel(model="ak135")
arrivals = m.get_ray_paths(...)
arrivals.plot()
In [ ]:
 
In [7]:
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();
/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/obspy/taup/tau_branch.py:496: UserWarning: Resizing a TauP array inplace failed due to the existence of other references to the array, creating a new array. See Obspy #2280.
  warnings.warn(msg)
/Users/lion/miniconda3/envs/seismo_live/lib/python3.7/site-packages/ipykernel_launcher.py:9: ObsPyDeprecationWarning: The plot() function is deprecated. Please use arrivals.plot_rays()
  if __name__ == '__main__':

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

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

print(first_arrival)
2011-03-11T05:58:52.936402Z

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

In [ ]:
 
In [9]:
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 [ ]: