ObsPy Tutorial
Downloading Data

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

Authors:

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

ObsPy has clients to directly fetch data via...

  • FDSN webservices (IRIS, Geofon/GFZ, USGS, NCEDC, SeisComp3 instances, ...)
  • ArcLink (EIDA, ...)
  • Earthworm
  • SeedLink (near-realtime servers)
  • NERIES/NERA/seismicportal.eu
  • NEIC
  • SeisHub (local seismological database)

This introduction shows how to use the FDSN webservice client. The FDSN webservice definition is by now the default web service implemented by many data centers world wide. Clients for other protocols work similar to the FDSN client.

Waveform Data

In [2]:
from obspy import UTCDateTime
from obspy.clients.fdsn import Client

client = Client("IRIS")
t = UTCDateTime("2011-03-11T05:46:23")  # Tohoku
st = client.get_waveforms("II", "PFO", "*", "LHZ",
                          t + 10 * 60, t + 30 * 60)
print(st)
st.plot()
2 Trace(s) in Stream:
II.PFO.00.LHZ | 2011-03-11T05:56:23.069500Z - 2011-03-11T06:16:22.069500Z | 1.0 Hz, 1200 samples
II.PFO.10.LHZ | 2011-03-11T05:56:23.069500Z - 2011-03-11T06:16:22.069500Z | 1.0 Hz, 1200 samples
Out[2]:
  • again, waveform data is returned as a Stream object
  • for all custom processing workflows it does not matter if the data originates from a local file or from a web service

Event Metadata

The FDSN client can also be used to request event metadata:

In [3]:
t = UTCDateTime("2011-03-11T05:46:23")  # Tohoku
catalog = client.get_events(starttime=t - 100, endtime=t + 24 * 3600,
                            minmagnitude=7)
print(catalog)
catalog.plot();
3 Event(s) in Catalog:
2011-03-11T06:25:50.740000Z | +38.051, +144.630 | 7.6 MW
2011-03-11T06:15:37.570000Z | +36.227, +141.088 | 7.9 MW
2011-03-11T05:46:23.200000Z | +38.296, +142.498 | 9.1 MW
/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)

Requests can have a wide range of constraints (see ObsPy Documentation):

  • time range
  • geographical (lonlat-box, circular by distance)
  • depth range
  • magnitude range, type
  • contributing agency

Station Metadata

Finally, the FDSN client can be used to request station metadata. Stations can be looked up using a wide range of constraints (see ObsPy documentation):

  • network/station code
  • time range of operation
  • geographical (lonlat-box, circular by distance)
In [4]:
event = catalog[0]
origin = event.origins[0]

# Münster
lon = 7.63
lat = 51.96

inventory = client.get_stations(longitude=lon, latitude=lat,
                                maxradius=2.5, level="station")
print(inventory)
Inventory created at 2019-11-07T13:52:13.000000Z
	Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.37
		    http://service.iris.edu/fdsnws/station/1/query?level=station&latitu...
	Sending institution: IRIS-DMC (IRIS-DMC)
	Contains:
		Networks (5):
			BE, GE, NL, SG, SY
		Stations (19):
			BE.MEM (Membach, Belgium)
			BE.RCHB (Rochefort, Belgium)
			BE.UCC (Uccle, Brussels, Belgium) (2x)
			GE.FLT1 (Temp GEOFON Station Flechtingen, Germany)
			GE.HLG (UKiel/GEOFON Station Helgoland, Germany) (2x)
			GE.IBBN (RUB/GEOFON Station Ibbenbueren, Germany)
			GE.WLF (GEOFON Station Walferdange, Luxembourg) (2x)
			NL.HGN (HEIMANSGROEVE, NETHERLANDS)
			SG.MEMB (Membach, Belgium)
			SG.RCHS (Rochefort, Belgium)
			SY.FLT1 (FLT1 synthetic)
			SY.HGN (HGN synthetic)
			SY.HLG (HLG synthetic)
			SY.IBBN (IBBN synthetic)
			SY.RCHB (RCHB synthetic)
			SY.WLF (WLF synthetic)
		Channels (0):

The level=... keyword is used to specify the level of detail in the requested inventory

  • "network": only return information on networks matching the criteria
  • "station": return information on all matching stations
  • "channel": return information on available channels in all stations networks matching the criteria
  • "response": include instrument response for all matching channels (large result data size!)
In [5]:
inventory = client.get_stations(network="OE", station="DAVA",
                                level="station")
print(inventory)
Inventory created at 2019-11-07T13:52:13.000000Z
	Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.37
		    http://service.iris.edu/fdsnws/station/1/query?network=OE&station=D...
	Sending institution: IRIS-DMC (IRIS-DMC)
	Contains:
		Networks (1):
			OE
		Stations (1):
			OE.DAVA (Damuels, Vorarlberg, Austria)
		Channels (0):

In [6]:
inventory = client.get_stations(network="OE", station="DAVA",
                                level="channel")
print(inventory)
Inventory created at 2019-11-07T13:52:14.000000Z
	Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.37
		    http://service.iris.edu/fdsnws/station/1/query?network=OE&station=D...
	Sending institution: IRIS-DMC (IRIS-DMC)
	Contains:
		Networks (1):
			OE
		Stations (1):
			OE.DAVA (Damuels, Vorarlberg, Austria)
		Channels (14):
			OE.DAVA..BHZ, OE.DAVA..BHN, OE.DAVA..BHE, OE.DAVA..HHZ, 
			OE.DAVA..HHN, OE.DAVA..HHE, OE.DAVA..LCQ, OE.DAVA..LHZ, 
			OE.DAVA..LHN, OE.DAVA..LHE, OE.DAVA..UFC, OE.DAVA..VHZ, 
			OE.DAVA..VHN, OE.DAVA..VHE

For waveform requests that include instrument correction, the appropriate instrument response information can be attached to waveforms automatically:
(Of course, for work on large datasets, the better choice is to download all station information and avoid the internal repeated webservice requests)

In [7]:
t = UTCDateTime("2011-03-11T05:46:23")  # Tohoku
st = client.get_waveforms("II", "PFO", "*", "LHZ",
                          t + 10 * 60, t + 30 * 60, attach_response=True)
st.plot()

st.remove_response()
st.plot()
Out[7]:

All data requested using the FDSN client can be directly saved to file using the filename="..." option. The data is then stored exactly as it is served by the data center, i.e. without first parsing by ObsPy and outputting by ObsPy.

In [8]:
client.get_events(starttime=t-100, endtime=t+24*3600, minmagnitude=7,
                  filename="/tmp/requested_events.xml")
client.get_stations(network="OE", station="DAVA", level="station",
                    filename="/tmp/requested_stations.xml")
client.get_waveforms("II", "PFO", "*", "LHZ", t + 10 * 60, t + 30 * 60,
                     filename="/tmp/requested_waveforms.mseed")
!ls -lrt /tmp/requested*
-rw-r--r--  1 lion  wheel   4028 Nov  7 14:52 /tmp/requested_events.xml
-rw-r--r--  1 lion  wheel   1423 Nov  7 14:52 /tmp/requested_stations.xml
-rw-r--r--  1 lion  wheel  16384 Nov  7 14:52 /tmp/requested_waveforms.mseed

FDSN Client Exercise

Use the FDSN client to assemble a waveform dataset for on event.

  • search for a large earthquake (e.g. by depth or in a region of your choice, use option limit=5 to keep network traffic down)
In [ ]:
 
In [9]:
from obspy.clients.fdsn import Client

client = Client()
catalog = client.get_events(minmagnitude=7, limit=5, mindepth=400)
print(catalog)
catalog.plot()
event = catalog[0]
print(event)
5 Event(s) in Catalog:
2018-09-06T15:49:14.420000Z | -18.495, +179.356 | 7.9 mww
2018-08-24T09:04:06.780000Z | -11.042,  -70.817 | 7.1 Mww
2018-08-19T00:19:40.670000Z | -18.113, -178.154 | 8.2 mww
2017-01-10T06:13:47.250000Z |  +4.463, +122.575 | 7.3 Mww
2015-11-24T22:50:54.370000Z | -10.060,  -71.018 | 7.6 mww
Event:	2018-09-06T15:49:14.420000Z | -18.495, +179.356 | 7.9 mww

	            resource_id: ResourceIdentifier(id="smi:service.iris.edu/fdsnws/event/1/query?eventid=10944928")
	             event_type: 'earthquake'
	    preferred_origin_id: ResourceIdentifier(id="smi:service.iris.edu/fdsnws/event/1/query?originid=34433936")
	 preferred_magnitude_id: ResourceIdentifier(id="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=189562093")
	                   ---------
	     event_descriptions: 1 Elements
	                origins: 1 Elements
	             magnitudes: 1 Elements
  • search for stations to look at waveforms for the event. stations should..
    • be available at the time of the event
    • have a vertical 1 Hz stream ("LHZ", to not overpower our network..)
    • be in a narrow angular distance around the event (e.g. 90-91 degrees)
    • adjust your search so that only a small number of stations (e.g. 3-6) match your search criteria
In [ ]:
 
In [10]:
origin = event.origins[0]
t = origin.time

inventory = client.get_stations(longitude=origin.longitude, latitude=origin.latitude,
                                minradius=101, maxradius=102,
                                starttime=t, endtime =t+100,
                                channel="LHZ", matchtimeseries=True)
print(inventory)
Inventory created at 2019-11-07T13:52:22.000000Z
	Created by: IRIS WEB SERVICE: fdsnws-station | version: 1.1.37
		    http://service.iris.edu/fdsnws/station/1/query?starttime=2018-09-06...
	Sending institution: IRIS-DMC (IRIS-DMC)
	Contains:
		Networks (6):
			C1, II, IU, N4, NM, US
		Stations (15):
			C1.TA01 (Daracena)
			II.HOPE (Hope Point, South Georgia Island)
			IU.LVC (Limon Verde, Chile)
			IU.OTAV (Otavalo, Ecuador)
			N4.146B (Union, MS, USA)
			N4.L40A (Anamosa, IA, USA)
			N4.N41A (Harden Midland Farm, Stronghurst, IL, ISA)
			N4.SPMN (Marine on St. Croix, MN, USA)
			N4.Y45B (Coffeeville, MS, USA)
			NM.FVM (French Village, MO)
			NM.MPH (Memphis, TN)
			NM.PBMO (Three Rivers Community College, Poplar Bluff, MO)
			NM.PVMO (UM Delta Research Portageville, MO)
			NM.SLM (St. Louis, MO)
			US.OXF (Oxford, Mississippi, USA)
		Channels (0):

  • for each of these stations download data of the event, e.g. a couple of minutes before to half an hour after the event
  • put all data together in one stream (put the get_waveforms() call in a try/except/pass block to silently skip stations that actually have no data available)
  • print stream info, plot the raw data
In [ ]:
 
In [11]:
from obspy import Stream
st = Stream()

for network in inventory:
    for station in network:
        try:
            st += client.get_waveforms(network.code, station.code, "*", "LHZ",
                                       t - 5 * 60, t + 30 * 60, attach_response=True)
        except:
            pass

print(st)
st.plot()
18 Trace(s) in Stream:
C1.TA01..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
II.HOPE.00.LHZ | 2018-09-06T15:44:15.069536Z - 2018-09-06T16:19:14.069536Z | 1.0 Hz, 2100 samples
II.HOPE.10.LHZ | 2018-09-06T15:44:15.069538Z - 2018-09-06T16:19:14.069538Z | 1.0 Hz, 2100 samples
IU.LVC.00.LHZ  | 2018-09-06T15:44:15.069538Z - 2018-09-06T16:19:14.069538Z | 1.0 Hz, 2100 samples
IU.LVC.10.LHZ  | 2018-09-06T15:44:15.069538Z - 2018-09-06T16:19:14.069538Z | 1.0 Hz, 2100 samples
IU.OTAV.00.LHZ | 2018-09-06T15:44:15.069536Z - 2018-09-06T16:19:14.069536Z | 1.0 Hz, 2100 samples
IU.OTAV.10.LHZ | 2018-09-06T15:44:15.069536Z - 2018-09-06T16:19:14.069536Z | 1.0 Hz, 2100 samples
N4.146B..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
N4.L40A..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
N4.N41A..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
N4.SPMN..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
N4.Y45B..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
NM.FVM..LHZ    | 2018-09-06T15:44:14.999998Z - 2018-09-06T16:19:13.999998Z | 1.0 Hz, 2100 samples
NM.MPH..LHZ    | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
NM.PBMO..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
NM.PVMO..LHZ   | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
NM.SLM.00.LHZ  | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
US.OXF.00.LHZ  | 2018-09-06T15:44:15.000000Z - 2018-09-06T16:19:14.000000Z | 1.0 Hz, 2100 samples
Out[11]:
  • correct the instrument response for all stations and plot the corrected data
In [ ]:
 
In [12]:
st.remove_response(water_level=20)
st.plot()
Out[12]:

If you have time, assemble and plot another similar dataset (e.g. like before stations at a certain distance from a big event, or use Transportable Array data for a big event, etc.)