Napojení ObsPy na datová centra (EIDA, IRIS)¶
Způsob napojení (klientské rozhraní)¶
- FDSN
- ArcLink
- SeedLink
- EarthWorm
- NEIC
- Syngine
Stahované produkty¶
- seismické signály
- staniční metadata
- zemětřesná metadata
- jiné (např. syntetické signály)
Datová centra s FDSN a ArcLink klienty¶
- IRIS
- EIDA: nódy ODC, GFZ, RESIF, INGV, ETHZ, BGR, NIEP, KOERI, NOA, (LMU)
- a další:
In [1]:
from obspy.clients.fdsn.header import URL_MAPPINGS
for key in sorted(URL_MAPPINGS.keys()):
print("{0:<7} {1}".format(key, URL_MAPPINGS[key]))
In [2]:
%pylab inline
from __future__ import print_function
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
import warnings
warnings.filterwarnings("ignore")
Clienty pro FDSN a Arclink přenosy si pro lepší přehlednost přejmenujeme:
In [3]:
from obspy.clients.fdsn import Client as Client_FDSN
from obspy.clients.arclink import Client as Client_Arclink
from obspy.core import UTCDateTime
Stažení zemětřesných metadat
Způsob stažení: Client_FDSN.get_events()
Datacentra poskytující eventy:
- ISC - obsahuje i finální ISC katalog
- USGS - obsahuje i ~ NEIC PDE katalog
IRISviz https://service.iris.edu/fdsnws/event/1/
- kompilace existujících katalogů
- poskytne eventy z ISC finálního katalogu, když existuje, a z NEIC PDE katalogu v období, kdy ISC ještě není
- ALE: není aktualizovaný a obsahuje díry!
Nadefinujeme si jiný způsob vypsání informací ze zemětřesných metadat než je print(catalog):
In [4]:
def print_cat(catalog):
"""
Jiný výpis než print(catalog)
"""
for event in catalog:
out = ''
origin = None
if event.origins:
origin = event.preferred_origin() or event.origins[0]
origin_time = origin.time.isoformat() + '.00'
out += '{:s} | {:+7.3f}, {:+8.3f}'.format(origin_time[:22], origin.latitude, origin.longitude)
if event.magnitudes:
magnitude = event.preferred_magnitude() or event.magnitudes[0]
out += ' | {:3.1f} {:<3s}'.format(magnitude.mag, magnitude.magnitude_type)
try:
origin_author = origin.creation_info.author or origin.creation_info.agency_id
origin_author = origin_author.upper()
except Exception:
origin_author = 'unk.'
try:
magn_author = magnitude.creation_info.author or magnitude.creation_info.agency_id
magn_author = magn_author.upper()
except Exception:
magn_author = 'unk.'
out += ' | {} ({}/{})'.format(event.event_descriptions[0].text, origin_author, magn_author)
print(out)
return
Můžeme předem filtrovat (podle období, oblasti, hloubky, magnituda), které eventy chceme stáhnout:
In [5]:
# PŘÍKLAD NA RADIUS
starttime = UTCDateTime('2010-01-01')
endtime = UTCDateTime()
# stanice PRA
latitude = 50.0692
longitude = 14.4277
# vzdálenost eventů od PRA (ve stupních)
minradius = 25
maxradius = 30
client_USGS = Client_FDSN('USGS')
catalog = client_USGS.get_events(starttime=starttime, endtime=endtime,
latitude=latitude, longitude=longitude,
minradius=minradius, maxradius=maxradius)
# vykreslení
plt.rcParams['figure.figsize'] = 12, 8
catalog.plot(projection="ortho"),
Out[5]:
Kompilace ISC finálního katalogu a NEIC PDE katalogu (aneb nahrazení stahování z nespolehlivého IRIS).
In [6]:
from obspy.core.event.catalog import Catalog
starttime = UTCDateTime('2016-04-01')
endtime = UTCDateTime('2016-06-01')
minmagnitude = 6
catalog = Catalog()
# ISC reviewed catalog
try:
client = Client_FDSN('ISC')
cat_ISC = client.get_events(starttime=starttime, endtime=endtime,
minmagnitude=minmagnitude, orderby="time-asc",
contributor="ISC") # jen finální ISC
# najdi měsíc, kde už ISC finální katalog není
last_event = cat_ISC[-1]
origin = last_event.preferred_origin() or last_event.origins[0]
year = origin.time.year
month = origin.time.month + 1
if month > 12:
year += 1
month = 1
starttime_US = UTCDateTime(year=year, month=month, day=1)
catalog += cat_ISC
except Exception:
starttime_US = starttime
# US catalog (~ NEIC PDE)
try:
client = Client_FDSN('USGS')
cat_US = client.get_events(starttime=starttime_US, endtime=endtime,
minmagnitude=minmagnitude, orderby="time-asc",
catalog="us") # jen ~ NEIC PDE
catalog += cat_US
except Exception:
pass
print_cat(catalog)
Stažení staničních metadat¶
Arclink¶
Umožňuje routing, umí pouze uložit Dataless SEED
from obspy.clients.arclink import Client
client = Client(user=[email])
client.save_response(route=True)</span>`
FDSN¶
Umí načíst metadata jako Inventory, případně je uložit jako StationXML.
Starý klient neumožňuje routing:
from obspy.clients.fdsn import Client
client = Client([datacentrum nebo nód])
client.get_stations()</span>
Nový klient umožňuje routing (ale neumí ukládat raw metadata do souborů):
from obspy.clients.fdsn import RoutingClient
client = RoutingClient(["eida-routing" | "iris-federator"])
client.get_stations()</span>`
Podíváme se, kde všude najdeme KHC:
In [7]:
from obspy.clients.fdsn import RoutingClient
client = RoutingClient("iris-federator")
inv_KHC = client.get_stations(network="CZ", station="KHC")
print("{}x : {}".format(len(inv_KHC.networks), inv_KHC.sender))
Načteme i odezvy seismometru a porovnáme metadata z dvou různých datacenter:
In [8]:
client_IRIS = Client_FDSN("IRIS")
inv_KHC_fromIRIS = client_IRIS.get_stations(network="CZ", station="KHC", level="response")
client_GFZ = Client_FDSN("GFZ")
inv_KHC_fromGFZ = client_GFZ.get_stations(network="CZ", station="KHC", level="response")
# Vykresli odezvy
fig = plt.figure(figsize=(16,8))
ax1a = plt.subplot2grid((2,2), (0,0))
ax1b = plt.subplot2grid((2,2), (1,0), sharex=ax1a)
ax2a = plt.subplot2grid((2,2), (0,1))
ax2b = plt.subplot2grid((2,2), (1,1), sharex=ax2a)
inv_KHC_fromIRIS.plot_response(0.01, channel="?HZ", axes=(ax1a,ax1b), show=False),
inv_KHC_fromGFZ.plot_response(0.01, channel="?HZ", axes=(ax2a,ax2b), show=False),
ax1a.legend(loc="lower center", ncol=3, fontsize='small')
ax2a.legend(loc="lower center", ncol=3, fontsize='small')
ax1a.set_ylim((1e-4,2e10)); ax2a.set_ylim((1e-4,2e10))
ax1a.set_title('IRIS',loc='left'); ax2a.set_title('GFZ',loc='left')
ax1b.set_xlabel('Frequency [Hz]'); ax2b.set_xlabel('Frequency [Hz]')
ax1a.set_ylabel('Amplitude'); ax1b.set_ylabel('Phase [rad]')
plt.show()
Filtrování staničních metadat:
eida-routing- pouzenetwork,station,location,channel,starttimeaendtimev současné verzi routinguiris-federator- plné fungování, podobně jako u eventů
In [9]:
from obspy.clients.fdsn import RoutingClient
client = RoutingClient("eida-routing", debug=True, timeout=360)
inv_Z3 = client.get_stations(network="CZ", station="KHC", location="", channel="?HZ")
#inv_Z3.plot(projection='local', color_per_network=True),
In [10]:
#from obspy.clients.fdsn import RoutingClient
client = RoutingClient("iris-federator")
inv_CR = client.get_stations(channel="*HZ", minlatitude=48.5, maxlatitude=51.0, minlongitude=12.0, maxlongitude=19)
inv_CR.plot(projection='local', color_per_network=True),
Out[10]:
In [12]:
### ZBYTEČNÉ, POKUD ZROVNA NEFUNGUJE ODC
#inv_CR_Z3 = client.get_stations(network="Z3", channel="*HZ",
# minlatitude=48.5, maxlatitude=51.0, minlongitude=12.0, maxlongitude=19)
#inv_CR_Z3.plot(projection='local', color_per_network=True),
Stažení seismických signálů
| routing | připnutí odezvy | raw mseed | |
|---|---|---|---|
| Arclink | ano | ne | ne |
| FDSN Client | ne | ano | ano |
| FDSN RoutingClient | ano | ne | ne |
Routing: automatické vyhledání příslušného datacentra/nódu
Připnutí odezvy: funguje jako attach_response
RAW mseed: uloží se do souboru v podobě, jak přišel z datacentra
Arclink: některé veřejně nepřístupné sítě umožňují stahování jen přes ArcLink (např. AlpArray)
Možné způsoby stahování¶
get_waveforms- základní filtrování pomocí network, station, location, channel, starttime, endtime (v location a channel fungují?*)get_waveforms_bulk- umožňuje dávkování požadavkůmass_downloader
- ucelený propracovaný modul z doby, kdy nefungoval FDSN RoutingClient
- prohledává všechny FDSN centra
- možnost stahovat signály i metadata
Porovnání signálů stažených pomocí RoutingFDSN a ArcLinkem:
In [13]:
tP = UTCDateTime('2015-08-06T09:34:46')
client = RoutingClient("eida-routing")
st_RFDSN = client.get_waveforms(network="CZ", station="*", channel="*HZ", starttime=tP-60, endtime=tP+60)
st_RFDSN.sort()
print(st_RFDSN[:6])
In [14]:
from obspy.core import Stream
stats = ['DPC','JAVC','KHC','KRLC','KRUC','MORC','NKC','OKC','PRA','PRU','PVCC','TREC','UPC','VRAC']
client = Client_Arclink(user="vecsey@ig.cas.cz")
st_Arclink = Stream()
for stat in stats:
st_Arclink += client.get_waveforms("CZ", stat, "*", "*HZ", starttime=tP-60, endtime=tP+60)
st_Arclink.sort()
print(st_Arclink[:6])
In [15]:
import numpy as np
def plot_section(st_orig, xlim=None, vline=None, title=None):
st = st_orig.copy()
t0 = UTCDateTime()
for i,tr in enumerate(st):
tr.detrend('linear')
tr.stats.distance = float(i)
if tr.stats.starttime < t0:
t0 = tr.stats.starttime
for tr in st:
n = int(round((tr.stats.starttime - t0) * tr.stats.sampling_rate))
pad = np.zeros(n) #; pad.fill(np.nan)
tr.data = np.insert(tr.data, 0, pad)
#tr.data = np.ma.masked_invalid(np.insert(tr.data, 0, pad))
tr.stats.starttime -= n * tr.stats.delta
fig = plt.figure(figsize=(15,9))
st.plot(fig=fig, type='section', orientation='horizontal', scale=2, alpha=0.8)
if vline:
for vl in vline:
fig.gca().axvline(vl-t0, color='r')
if xlim:
fig.gca().set_xlim(xlim)
if title:
fig.gca().set_title(title)
return
In [16]:
plot_section(st_RFDSN, (0,180), (tP-60,tP+60), 'FDSN')
plot_section(st_Arclink, (-34,146), (tP-60,tP+60), 'ArcLink')
ZDE NOTEBOOK UKONČÍME
HINT: Umístění např. AlpArray Z3 stanic v EIDA nódech lze získat i přes webovské rozhraní:
http://www.orfeus-eu.org/eidaws/routing/1/query?net=Z3&start=2014-01-01&format=post
In [17]:
# HINT: které služby poskytuje client:
Client_FDSN("USGS").help('event')
In [19]:
# zastaralost eventů z IRIS
client = Client_FDSN("IRIS")
cat = client.get_events(starttime=UTCDateTime("2014-06-01"), endtime=UTCDateTime("2014-08-01"),
minmagnitude=6, orderby="time-asc")
print_cat(cat)