Commit 952ed581 authored by Mijian Xu's avatar Mijian Xu 😷
Browse files

Add taup

parent 0f6c6264
...@@ -18,3 +18,5 @@ doc/.texpadtmp/bqmail.out ...@@ -18,3 +18,5 @@ doc/.texpadtmp/bqmail.out
doc/.texpadtmp/bqmail.synctex.gz doc/.texpadtmp/bqmail.synctex.gz
doc/.texpadtmp/bqmail.toc doc/.texpadtmp/bqmail.toc
bin/bqmail_conti
...@@ -8,11 +8,12 @@ ...@@ -8,11 +8,12 @@
# 2015/02/11 # 2015/02/11
# 2015/04/29 # 2015/04/29
# 2015/05/01 # 2015/05/01
# 2015/09/26
# #
def Usage(): def Usage():
print('Usage:') print('Usage:')
print('python bqmail.py -Nnetwork -Sstation -Yyear1/month1/day1/year2/month2/day2 -Bsec_begin/sec_end -Cchannel -Llocation -cdatetimefile -Fformat head.cfg') print('python bqmail.py -Nnetwork -Sstation -Yyear1/month1/day1/year2/month2/day2 -Bsec_begin/sec_end [-Cchannel] [-Plat/lon/phase] [-Llocation] [-cdatetimefile] [-Fformat] head.cfg')
print('-N -- Network.') print('-N -- Network.')
print('-S -- Station.') print('-S -- Station.')
print('-Y -- Date range.') print('-Y -- Date range.')
...@@ -29,8 +30,9 @@ def Usage(): ...@@ -29,8 +30,9 @@ def Usage():
import datetime import datetime
import os, re import os, re
import sys, getopt import sys, getopt
from smtplib import SMTP
import time import time
import taup
import distaz
try: try:
import configparser import configparser
config = configparser.ConfigParser() config = configparser.ConfigParser()
...@@ -40,7 +42,7 @@ except: ...@@ -40,7 +42,7 @@ except:
try: try:
opts,args = getopt.getopt(sys.argv[1:], "hN:S:C:Y:B:L:c:F:") opts,args = getopt.getopt(sys.argv[1:], "hN:S:C:Y:B:L:c:F:P:")
except: except:
print('Arguments are not found!') print('Arguments are not found!')
Usage() Usage()
...@@ -49,7 +51,7 @@ if opts == []: ...@@ -49,7 +51,7 @@ if opts == []:
Usage() Usage()
sys.exit(1) sys.exit(1)
isph = 0
iscustom = 0 iscustom = 0
isyrange = 0 isyrange = 0
chan = "BH?" chan = "BH?"
...@@ -74,6 +76,12 @@ for op, value in opts: ...@@ -74,6 +76,12 @@ for op, value in opts:
loca = value loca = value
elif op == "-F": elif op == "-F":
fformat = value fformat = value
elif op == "-P":
stla = float(value.split('/')[0])
stlo = float(value.split('/')[1])
phase = value.split('/')[2]
isph = 1
mod=taup.TauPyModel(model='iasp91')
elif op == "-h": elif op == "-h":
Usage() Usage()
sys.exit(1) sys.exit(1)
...@@ -110,9 +118,6 @@ INST = config.get("info","INST") ...@@ -110,9 +118,6 @@ INST = config.get("info","INST")
EMAIL = config.get("info","EMAIL") EMAIL = config.get("info","EMAIL")
MEDIA = config.get("info","MEDIA") MEDIA = config.get("info","MEDIA")
ALTERNATEMEDIA = MEDIA ALTERNATEMEDIA = MEDIA
hosts = config.get("smtp","hosts")
port = config.get("smtp","port")
passwd = config.get("smtp","passwd")
if fformat.lower() == 'seed': if fformat.lower() == 'seed':
recipient = 'breq_fast@iris.washington.edu' recipient = 'breq_fast@iris.washington.edu'
elif fformat.lower() == 'miniseed': elif fformat.lower() == 'miniseed':
...@@ -150,14 +155,20 @@ else: ...@@ -150,14 +155,20 @@ else:
lon=float(evenum_split[8]) lon=float(evenum_split[8])
dep=float(evenum_split[9]) dep=float(evenum_split[9])
mw=float(evenum_split[10]) mw=float(evenum_split[10])
date = datetime.datetime(year,mon,day,hour,min,sec) + datetime.timedelta(seconds=btime) evt_time = datetime.datetime(year,mon,day,hour,min,sec)
dateend = datetime.datetime(year,mon,day,hour,min,sec) + datetime.timedelta(seconds=etime) if datemin <= evt_time <= datemax:
if datemin <= date <= datemax: if isph == 1:
event.append([date.strftime('%Y'),date.strftime('%m'),date.strftime('%d'),date.strftime('%H'),date.strftime('%M'),date.strftime('%S'),dateend.strftime('%Y'),dateend.strftime('%m'),dateend.strftime('%d'),dateend.strftime('%H'),dateend.strftime('%M'),dateend.strftime('%S')]) dis = distaz.distaz(stla, stlo, lat, lon).delta
arr = mod.get_travel_times(source_depth_in_km=dep, distance_in_degree=dis, phase_list=[phase])
if len(arr) != 0:
head = ("From: %s\r\nTo: %s\r\n\r\n" % (EMAIL, recipient)) arr_time = evt_time + datetime.timedelta(seconds=arr[0].time)
msg = head date = arr_time - datetime.timedelta(seconds=btime)
dateend = arr_time + datetime.timedelta(seconds=etime)
else:
date = evt_time + datetime.timedelta(seconds=btime)
dateend = evt_time + datetime.timedelta(seconds=etime)
event.append([date.strftime('%Y %m %d %H %M %S'), dateend.strftime('%Y %m %d %H %M %S')])
msg = ''
msg += '.NAME '+NAME+'\n' msg += '.NAME '+NAME+'\n'
msg += '.INST '+INST+'\n' msg += '.INST '+INST+'\n'
msg += '.MAIL\n' msg += '.MAIL\n'
...@@ -170,13 +181,11 @@ msg += '.ALTERNATE MEDIA '+ALTERNATEMEDIA+'\n' ...@@ -170,13 +181,11 @@ msg += '.ALTERNATE MEDIA '+ALTERNATEMEDIA+'\n'
msg += '.LABEL '+LABEL+'\n' msg += '.LABEL '+LABEL+'\n'
msg += '.END\n' msg += '.END\n'
for row in event: for row in event:
msg += station+' '+network+' '+row[0]+' '+row[1]+' '+row[2]+' '+row[3]+' '+row[4]+' '+row[5]+' '+row[6]+' '+row[7]+' '+row[8]+' '+row[9]+' '+row[10]+' '+row[11]+' 1 '+chan+' '+loca+'\n' msg += station+' '+network+' '+row[0]+' '+row[1]+' 1 '+chan+' '+loca+'\n'
with open('tmp.bq','w') as fid_msg:
smtp = SMTP(host=hosts, port=port) fid_msg.write(msg)
smtp.set_debuglevel(0) os.system('mail '+recipient+'<tmp.bq')
smtp.login(EMAIL, passwd)
smtp.sendmail(EMAIL, recipient, msg)
smtp.quit()
time.sleep(5)
print("Successful sending the mail of "+network+"."+station+" to IRIS DMC!!!") print("Successful sending the mail of "+network+"."+station+" to IRIS DMC!!!")
os.system('rm tmp.bq')
time.sleep(4)
...@@ -16,13 +16,12 @@ def Usage(): ...@@ -16,13 +16,12 @@ def Usage():
print('-H -- Request continuous wave by hour.') print('-H -- Request continuous wave by hour.')
print('-F -- File format (SEED or miniseed). Default: SEED') print('-F -- File format (SEED or miniseed). Default: SEED')
print('head.cfg -- Config file.') print('head.cfg -- Config file.')
print('Example: ./bqmail_conti.py -Iex_sta.lst -Y2003/12/3/2003/12/3 -H24 head.cfg') print('Example: ./bqmail_conti.py -Iex_sta.lst -Y2003/12/3/2003/12/4 -H24 head.cfg')
import datetime import datetime
import os, re import os, re
import sys, getopt import sys, getopt
from smtplib import SMTP
import time import time
try: try:
import configparser import configparser
...@@ -91,9 +90,6 @@ INST = config.get("info","INST") ...@@ -91,9 +90,6 @@ INST = config.get("info","INST")
EMAIL = config.get("info","EMAIL") EMAIL = config.get("info","EMAIL")
MEDIA = config.get("info","MEDIA") MEDIA = config.get("info","MEDIA")
ALTERNATEMEDIA = MEDIA ALTERNATEMEDIA = MEDIA
hosts = config.get("smtp","hosts")
port = config.get("smtp","port")
passwd = config.get("smtp","passwd")
if fformat.lower() == 'seed': if fformat.lower() == 'seed':
recipient = 'breq_fast@iris.washington.edu' recipient = 'breq_fast@iris.washington.edu'
elif fformat.lower() == 'miniseed': elif fformat.lower() == 'miniseed':
...@@ -112,18 +108,13 @@ for stainfo in fid.readlines(): ...@@ -112,18 +108,13 @@ for stainfo in fid.readlines():
else: else:
sta.append([stainfo_sp[0], stainfo_sp[1], '']) sta.append([stainfo_sp[0], stainfo_sp[1], ''])
smtp = SMTP(host=hosts, port=port)
smtp.set_debuglevel(0)
smtp.login(EMAIL, passwd)
nowtime = datemin nowtime = datemin
while 1: while 1:
if nowtime >= datemax: if nowtime >= datemax:
break break
endtime = nowtime + datetime.timedelta(hours=timeval) endtime = nowtime + datetime.timedelta(hours=timeval)
LABEL = 'IRIS_'+nowtime.strftime('%Y')+'.'+nowtime.strftime('%m')+'.'+nowtime.strftime('%d')+'.'+nowtime.strftime('%H') LABEL = 'IRIS_'+nowtime.strftime('%Y')+'.'+nowtime.strftime('%m')+'.'+nowtime.strftime('%d')+'.'+nowtime.strftime('%H')
head = ("From: %s\r\nTo: %s\r\n\r\n" % (EMAIL, recipient)) msg = ''
msg = head
msg += '.NAME '+NAME+'\n' msg += '.NAME '+NAME+'\n'
msg += '.INST '+INST+'\n' msg += '.INST '+INST+'\n'
msg += '.MAIL\n' msg += '.MAIL\n'
...@@ -136,9 +127,11 @@ while 1: ...@@ -136,9 +127,11 @@ while 1:
msg += '.LABEL '+LABEL+'\n' msg += '.LABEL '+LABEL+'\n'
msg += '.END\n' msg += '.END\n'
for sta_row in sta: for sta_row in sta:
msg += sta_row[1]+' '+sta_row[0]+' '+nowtime.strftime('%Y')+' '+nowtime.strftime('%m')+' '+nowtime.strftime('%d')+' '+nowtime.strftime('%H')+' '+nowtime.strftime('%M')+' 00.0 '+endtime.strftime('%Y')+' '+endtime.strftime('%m')+' '+endtime.strftime('%d')+' '+endtime.strftime('%H')+' '+endtime.strftime('%M')+' 00.0 1 '+chan+' '+sta_row[2]+'\n' msg += sta_row[1]+' '+sta_row[0]+' '+nowtime.strftime('%Y %m %d %H %M %S')+' '+endtime.strftime('%Y %m %d %H %M %S')+' 1 '+chan+' '+sta_row[2]+'\n'
smtp.sendmail(EMAIL, recipient, msg) with open('tmp.bq','w') as fid_msg:
time.sleep(5) fid_msg.write(msg)
os.system('mail '+recipient+'<tmp.bq')
print("Successful sending the mail between "+nowtime.strftime('%Y')+'.'+nowtime.strftime('%m')+'.'+nowtime.strftime('%d')+'.'+nowtime.strftime('%H')+" and "+endtime.strftime('%Y')+'.'+endtime.strftime('%m')+'.'+endtime.strftime('%d')+'.'+endtime.strftime('%H')+"!!!") print("Successful sending the mail between "+nowtime.strftime('%Y')+'.'+nowtime.strftime('%m')+'.'+nowtime.strftime('%d')+'.'+nowtime.strftime('%H')+" and "+endtime.strftime('%Y')+'.'+endtime.strftime('%m')+'.'+endtime.strftime('%d')+'.'+endtime.strftime('%H')+"!!!")
nowtime = nowtime + datetime.timedelta(hours=timeval) nowtime = nowtime + datetime.timedelta(hours=timeval)
smtp.quit() time.sleep(4.5)
os.system('rm tmp.bq')
...@@ -3,7 +3,7 @@ mkdir ./bin ...@@ -3,7 +3,7 @@ mkdir ./bin
ln -s `pwd`/bqmail.py `pwd`/bin/bqmail ln -s `pwd`/bqmail.py `pwd`/bin/bqmail
ln -s `pwd`/searchDMC.py `pwd`/bin/searchDMC ln -s `pwd`/searchDMC.py `pwd`/bin/searchDMC
ln -s `pwd`/updateCatalog.py `pwd`/bin/updateCatalog ln -s `pwd`/updateCatalog.py `pwd`/bin/updateCatalog
ln -s `pwd`/bqmail_continue.py `pwd`/bin/bqmail_continue ln -s `pwd`/bqmail_conti.py `pwd`/bin/bqmail_conti
if [ `uname` == "Darwin" ]; then if [ `uname` == "Darwin" ]; then
echo "export PATH=`pwd`/bin:\$PATH" >> ~/.bash_profile echo "export PATH=`pwd`/bin:\$PATH" >> ~/.bash_profile
else else
......
...@@ -8,7 +8,8 @@ ...@@ -8,7 +8,8 @@
# 2015/02/11 # 2015/02/11
# 2015/04/29 # 2015/04/29
# 2015/07/19 # 2015/07/19
# 2015/09/26
#
import distaz, math import distaz, math
try: try:
...@@ -21,7 +22,7 @@ import sys, getopt ...@@ -21,7 +22,7 @@ import sys, getopt
import glob import glob
def Usage(): def Usage():
print('Usage: python searchDMC.py -NNetwork -Sstation -Rlon1/lon2/lat1/lat2 -Dlon/lat/dis1/dis2 -Yyear1/mon1/day1/year2/mon2/day2 -Cchannel -K') print('Usage: python searchDMC.py -NNetwork -Sstation -Rlon1/lon2/lat1/lat2 -Dlat/lon/dis1/dis2 -Yyear1/mon1/day1/year2/mon2/day2 -Cchannel -K -G')
print('-N -- Network.') print('-N -- Network.')
print('-S -- Station.') print('-S -- Station.')
print('-R -- Search range.') print('-R -- Search range.')
...@@ -32,7 +33,7 @@ def Usage(): ...@@ -32,7 +33,7 @@ def Usage():
try: try:
opts,args = getopt.getopt(sys.argv[1:], "hR:D:KY:C:N:S:") opts,args = getopt.getopt(sys.argv[1:], "hR:D:KY:C:N:S:G")
except: except:
print('Arguments are not found!') print('Arguments are not found!')
Usage() Usage()
...@@ -44,6 +45,7 @@ if opts == []: ...@@ -44,6 +45,7 @@ if opts == []:
iskml = 0 iskml = 0
islalo = 0 islalo = 0
isyrange = 0 isyrange = 0
isgmt = 0
lat_lon = '' lat_lon = ''
yrange = '' yrange = ''
chan = '' chan = ''
...@@ -67,6 +69,8 @@ for op, value in opts: ...@@ -67,6 +69,8 @@ for op, value in opts:
isyrange = 1 isyrange = 1
elif op == "-C": elif op == "-C":
chan = 'chan='+value+'&' chan = 'chan='+value+'&'
elif op == "-G":
isgmt = 1
elif op == "-h": elif op == "-h":
Usage() Usage()
sys.exit(1) sys.exit(1)
...@@ -84,8 +88,8 @@ if lat_lon != '': ...@@ -84,8 +88,8 @@ if lat_lon != '':
lalo = lon1+'_'+lon2+'_'+lat1+'_'+lat2 lalo = lon1+'_'+lon2+'_'+lat1+'_'+lat2
lalo_label = 'minlat='+lat1+'&maxlat='+lat2+'&minlon='+lon1+'&maxlon='+lon2+'&' lalo_label = 'minlat='+lat1+'&maxlat='+lat2+'&minlon='+lon1+'&maxlon='+lon2+'&'
else: else:
lon = lat_lon_split[0] lon = lat_lon_split[1]
lat = lat_lon_split[1] lat = lat_lon_split[0]
dis1 = float(lat_lon_split[2]) dis1 = float(lat_lon_split[2])
dis2 = float(lat_lon_split[3]) dis2 = float(lat_lon_split[3])
lon1 = str(0) lon1 = str(0)
...@@ -106,12 +110,12 @@ if isyrange: ...@@ -106,12 +110,12 @@ if isyrange:
day2 = yrange_sp[5] day2 = yrange_sp[5]
yrange = 'timewindow='+year1+'/'+mon1+'/'+day1+'-'+year2+'/'+mon2+'/'+day2+'&' yrange = 'timewindow='+year1+'/'+mon1+'/'+day1+'-'+year2+'/'+mon2+'/'+day2+'&'
url += network+station+lalo_label+yrange+chan url += network+station+lalo_label+yrange+chan+'archive=on'
url = url[0:-1]
response = rq.urlopen(url) response = rq.urlopen(url)
html = str(response.read()) html = str(response.read())
find_re = re.compile(r'<station\s.+?"/>',re.DOTALL) find_re = re.compile(r'<station\s.+?"/>',re.DOTALL)
stations = []
for info in find_re.findall(html): for info in find_re.findall(html):
sta_info = re.split('\w+="|"\s+?\w+="|"\s/>',info) sta_info = re.split('\w+="|"\s+?\w+="|"\s/>',info)
if sta_info == []: if sta_info == []:
...@@ -127,9 +131,33 @@ for info in find_re.findall(html): ...@@ -127,9 +131,33 @@ for info in find_re.findall(html):
if not islalo and lat_lon != '': if not islalo and lat_lon != '':
delta = distaz.distaz(float(lat),float(lon),float(stlat),float(stlon)) delta = distaz.distaz(float(lat),float(lon),float(stlat),float(stlon))
if dis1 < delta.delta < dis2: if dis1 < delta.delta < dis2:
print(netname+' '+staname+' '+stlat+' '+stlon+' '+yrange1+' '+yrange2) stations.append([netname, staname, float(stlat), float(stlon), yrange1, yrange2])
else: else:
print(netname+' '+staname+' '+stlat+' '+stlon+' '+yrange1+' '+yrange2) stations.append([netname, staname, float(stlat), float(stlon), yrange1, yrange2])
for station in stations:
print('%s %s %5.2f %5.2f %s %s' % (station[0], station[1], station[2], station[3], station[4], station[5]))
if isgmt:
with open('Stations.gmt', 'w+') as gmt:
if not islalo and lat_lon != '':
center = [float(lat), float(lon)]
elif islalo:
center = [(float(lat1)+float(lat2))/2, (float(lon1)+float(lon2))/2]
else:
netlats = []
netlons = []
for sta in stations:
netlats.append(sta[2])
netlons.append(sta[3])
center = [(max(netlats)+min(netlats))/2, (max(netlons)+min(netlons))/2]
gmt.write('#!/bin/sh\n')
gmt.write('ps=stations.ps\n\n')
gmt.write('gmt pscoast -Rg -JE%5.2f/%5.2f/5i -Ba30 -Dc -A10000 -Glightgray -Wthinnest -K > $ps\n' % (center[1], center[0]))
gmt.write('gmt psxy -R -K -O -J -St0.03i -Gred3 -W0.3p >> $ps << eof\n')
for sta in stations:
gmt.write('%5.2f %5.2f\n' % (sta[3], sta[2]))
gmt.write('eof\n')
if iskml: if iskml:
if network != '': if network != '':
......
package obspy.taup
==================
Copyright
---------
GLPv3
Copyright (c) 2011-2015 by:
* The ObsPy Dev Team
Overview
--------
Travel time calculation tool for ObsPy.
ObsPy is an open-source project dedicated to provide a Python framework for
processing seismological data. It provides parsers for common file formats and
seismological signal processing routines which allow the manipulation of
seismological time series (see Beyreuther et al. 2010, Megies et al. 2011).
The goal of the ObsPy project is to facilitate rapid application development
for seismology.
For more information visit http://www.obspy.org.
# -*- coding: utf-8 -*-
"""
obspy.taup - Ray Theoretical Travel Times and Paths
===================================================
:copyright:
The ObsPy Development Team (devs@obspy.org)
:license:
GNU Lesser General Public License, Version 3
(http://www.gnu.org/copyleft/lesser.html)
This package started out as port of the Java TauP Toolkit by [Crotwell1999]_ so
please look there for more details about the algorithms used and further
information. It can be used to calculate theoretical arrival times for
arbitrary seismic phases in a 1D spherically symmetric background model.
Furthermore it can output ray paths for all phases and derive pierce points of
rays with model discontinuities.
Basic Usage
-----------
Let's start by initializing a :class:`~obspy.taup.tau.TauPyModel` instance.
Models can be initialized by specifying the name of a model provided by ObsPy.
Names of available builtin models (in ``obspy/taup/data`` folder) are provided
by :const:`~obspy.taup.BUILTIN_TAUP_MODELS`.
>>> from obspy.taup import TauPyModel
>>> model = TauPyModel(model="iasp91")
Model initialization is a fairly expensive operation so make sure to do it only
if necessary. Custom built models can be initialized by specifying an absolute
path to a model in ObsPy's ``.npz`` model format instead of just a model name.
See below for how to build a ``.npz`` model file.
Travel Times
^^^^^^^^^^^^
The models' main method is the
:meth:`~obspy.taup.tau.TauPyModel.get_travel_times` method; as the name
suggests it returns travel times for the chosen phases, distance, source depth,
and model. Per default it returns arrivals for a number of phases.
>>> arrivals = model.get_travel_times(source_depth_in_km=55,
... distance_in_degree=67)
>>> print(arrivals) # doctest: +NORMALIZE_WHITESPACE
28 arrivals
P phase arrival at 647.036 seconds
pP phase arrival at 662.230 seconds
sP phase arrival at 668.702 seconds
PcP phase arrival at 674.868 seconds
PP phase arrival at 794.975 seconds
PKiKP phase arrival at 1034.106 seconds
pPKiKP phase arrival at 1050.535 seconds
sPKiKP phase arrival at 1056.727 seconds
S phase arrival at 1176.947 seconds
pS phase arrival at 1195.500 seconds
SP phase arrival at 1196.827 seconds
sS phase arrival at 1203.128 seconds
PS phase arrival at 1205.418 seconds
SKS phase arrival at 1239.088 seconds
SKKS phase arrival at 1239.107 seconds
ScS phase arrival at 1239.515 seconds
SKiKP phase arrival at 1242.400 seconds
pSKS phase arrival at 1260.313 seconds
sSKS phase arrival at 1266.919 seconds
SS phase arrival at 1437.417 seconds
PKIKKIKP phase arrival at 1855.260 seconds
SKIKKIKP phase arrival at 2063.556 seconds
PKIKKIKS phase arrival at 2069.749 seconds
SKIKKIKS phase arrival at 2277.833 seconds
PKIKPPKIKP phase arrival at 2353.930 seconds
PKPPKP phase arrival at 2356.420 seconds
PKPPKP phase arrival at 2358.925 seconds
SKIKSSKIKS phase arrival at 3208.154 seconds
If you know which phases you are interested in, you can also specify them
directly which speeds up the calculation as unnecessary phases are not
calculated. Please note that it is possible to construct *any* phases that
adhere to the naming scheme which is detailed later.
>>> arrivals = model.get_travel_times(source_depth_in_km=100,
... distance_in_degree=45,
... phase_list=["P", "PSPSPS"])
>>> print(arrivals) # doctest: +NORMALIZE_WHITESPACE
3 arrivals
P phase arrival at 485.204 seconds
PSPSPS phase arrival at 4983.023 seconds
PSPSPS phase arrival at 5799.225 seconds
Each arrival is represented by an :class:`~obspy.taup.helper_classes.Arrival`
object which can be queried for various attributes.
>>> arr = arrivals[0]
>>> arr.ray_param, arr.time, arr.incident_angle # doctest: +ELLIPSIS
(453.7188..., 485.2041..., 24.3968...)
Ray Paths
^^^^^^^^^
To also calculate the paths travelled by the rays to the receiver, use the
:meth:`~obspy.taup.tau.TauPyModel.get_ray_paths` method.
>>> arrivals = model.get_ray_paths(500, 130)
>>> arrival = arrivals[0]
The result is a NumPy record array containing ray parameter, time, distance
and depth to use however you see fit.
>>> arrival.path.dtype
dtype([('p', '<f8'), ('time', '<f8'), ('dist', '<f8'), ('depth', '<f8')])
Pierce Points
^^^^^^^^^^^^^
If you only need the pierce points of ray paths with model discontinuities,
use the :meth:`~obspy.taup.tau.TauPyModel.get_pierce_points` method which
results in pierce points being stored as a record array on the arrival object.
>>> arrivals = model.get_pierce_points(500, 130)
>>> arrivals[0].pierce.dtype
dtype([('p', '<f8'), ('time', '<f8'), ('dist', '<f8'), ('depth', '<f8')])
Plotting
--------
If ray paths have been calculated, they can be plotted using the
:meth:`~obspy.taup.tau.Arrivals.plot` method:
>>> arrivals = model.get_ray_paths(source_depth_in_km=500,
... distance_in_degree=130)
>>> arrivals.plot() # doctest: +SKIP
.. plot::
:width: 50%
:align: center
from obspy.taup import TauPyModel
TauPyModel().get_ray_paths(500, 130).plot()
Plotting will only show the requested phases:
>>> arrivals = model.get_ray_paths(
... source_depth_in_km=500,
... distance_in_degree=130,
... phase_list=["Pdiff", "Sdiff", "pPdiff", "sSdiff"])
>>> arrivals.plot() # doctest: +SKIP
.. plot::
:width: 50%
:align: center
from obspy.taup import TauPyModel
TauPyModel().get_ray_paths(
500, 130,
phase_list=["Pdiff", "Sdiff", "pPdiff", "sSdiff"]).plot()
Additionally, Cartesian coordinates may be used instead of a polar grid:
>>> arrivals = model.get_ray_paths(source_depth_in_km=500,
... distance_in_degree=130,
... phase_list=["ttbasic"])
>>> arrivals.plot(plot_type="cartesian") # doctest: +SKIP
.. plot::
:width: 75%
:align: center
from obspy.taup import TauPyModel
TauPyModel().get_ray_paths(
500, 130, phase_list=["ttbasic"]).plot(plot_type="cartesian")
More examples of plotting may be found in the :doc:`ObsPy tutorial
</tutorial/code_snippets/travel_time>`.
Phase naming in obspy.taup
--------------------------
.. note::
This section is a modified copy from the Java TauP Toolkit documentation so
all credit goes to the authors of that.
A major feature of ``obspy.taup`` is the implementation of a phase name parser
that allows the user to define essentially arbitrary phases through the Earth.
Thus, ``obspy.taup`` is extremely flexible in this respect since it is not