1 Obtaining Seismic Data Using Obspy

1.1 Basic usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime

# Create a data service from 'IRIS' data center.
client = Client('IRIS')

# Setting the start time of waveforms.
t1 = UTCDateTime('2019-09-04T00:00:00')

# End time of waveforms: 3600 seconds shift relative to t1.
t2 = t1 + 3600

# Get waveforms. 
st = client.get_waveforms(network='CO', station='CSB',
						  location='*', channel='LHZ',
						  starttime=t1, endtime=t2,
						  attach_response=False)

You can set attach_response as True if you want to download instrument response. However, the internet service may go wrong when you download data of many stations. Fortunately, python offers us the exception handling to deal with this situation. For instance, the following python codes give you a very simple example to fix it.

1.2 Fix internet service issue

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

for sta in stations:
	try:
		st = client.get_waveforms(network='CO', station=sta,
							  location='*', channel='LHZ',
							  starttime=t1, endtime=t2,
							  attach_response=False)
	except Exception:
		print('....')
	finally:
		pass

1.3 Seismic data management center

Change the data center code if you want to download data from other organizations. These centers are

1
2
3
4

from obspy.clients.fdsn.header import URL_MAPPINGS

for k in URL_MAPPINGS.keys():
    print(k)

Codes of data service centers:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

BGR
ETH
EMSC
GEONET
GFZ
ICGC
INGV
IPGP
IRIS
ISC
KNMI
KOERI
LMU
NCEDC
NIEP
NOA
ODC
ORFEUS
RESIF
RASPISHAKE
SCEDC
TEXNET
UIB-NORSAR
USGS
USP

1.4 Other data center(s)

However, seismic data in Australian seismic network are not included in the above data centers,

and here we show you to retrieve seismic data from AusPass.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

import obspy
from obspy.clients.fdsn import Client as FDSN_Client
from obspy.core import UTCDateTime
from matplotlib import pyplot as plt

auspass = FDSN_Client('http://auspass.edu.au:8080',user_agent='email@address.com')

t1 = UTCDateTime('2018-03-29T00:00:00')
t2 = t1 + 500
st = auspass.get_waveforms(network='S1', station='AUAYR', location='*',
                           channel='BHZ', starttime=t1, endtime=t2)
st.plot()

1.5 Active source data PH5

For controlled (active) source seismic experiment data, you can request from PH5WS. Here we show you how to obtain active source seismic data using Obspy. For more details, you can find them from https://github.com/PIC-IRIS/PH5/wiki/PH5-Web-Services-Shot-Gather.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

import obspy
from obspy.clients.fdsn.client import Client
from obspy.core import UTCDateTime


DATASELECT = 'http://service.iris.edu/ph5ws/dataselect/1'
client = Client(service_mappings={'dataselect': DATASELECT},
               debug=False)

t1 = UTCDateTime('2017-08-25T00:00:00')
t2 = t1 + 2
st = client.get_waveforms(network='YI', station='60', 
                        location='*', channel='DPZ',
                        starttime=t1, endtime=t2)

st.plot()

1.6 The Swiss Seismological Service (SED)

Using standard FDSN web service for the SED, we can directly download seismic data from the data center SED. For example, we download whole-day BHZ seismic data recorded by the station ACB of network CH:

1

wget "http://eida.ethz.ch/fdsnws/dataselect/1/query?net=CH&sta=ACBT&loc=--&cha=BHZ&starttime=2019-01-01T00:00:00&endtime=2019-01-02T00:00:00&nodata=404" -O 2019-01-01_ACB_BHZ.mseed

The wildcard *, of course, can be supported by this style, and here we retrieve whole-day BHZ seismic data recorded by all stations of the network CH:

1

wget "http://eida.ethz.ch/fdsnws/dataselect/1/query?net=CH&sta=*&loc=--&cha=BHZ&starttime=2019-01-01T00:00:00&endtime=2019-01-02T00:00:00&nodata=404" -O 2019-01-01_CH_ALL_BHZ.mseed

2 Remove Instrument Response

stream ( or trace).remove_response(pre_filt=[f1, f2, f3, f4], output='DISP') The parameter output defines the output physical quantity.
DISP: displacement in meters;
VEL: velocity in meters per second;
ACC: acceleration in meters per squared seconds.
You can set the 4 corner frequencies in parameter pre_filt if you want to do pre-filtering on your data, but be careful when you choose frequencies that are not over Nyquist frequency.

3 Get Information about Stations

1
2
3
4

st = client.stations(network='CO', station='CSB',
                           location='*', channel='LHZ',
                           starttime=t1, endtime=t2,
                           level='station')

The parameter level defines what kind of returning information.
network: network code;
station: station code;
channel: channel code;
response: instrument response.
Usually, we set level as response to download instrument response for the later data processing.

4 Solution to Empty sac attribution of SAC

This problem caused by that there is no attribution in trace.stats. You can use AttribDict to add the sac field to trace.stats. Here, we suppose that you have a trace (abbreviated to tr).

1
2
3

from obspy import AttribDict

tr.stats.sac = AttribDict({})

5 Calculating the distance between given two points

Supposing that you’ve installed the package geopy in your python.

1
2
3
4
5
6
7

from geopy.distance import geodesic

# Point format: (latitude, longitude).
p1 = (5, 5)
p2 = (10, 10)
# Calculating the distance in kilometers with 'WGS-84' ellisoid model.
d = geodesic(p1, p2, ellipsoid=''WGS-84').km

Bibliography

http://auspass.edu.au/help/obspy_request.html