Calculate properties of teleseismic arrivals through a selection of 1D Earth models, using the ObsPy Python software.
To install on Julia versions v0.7 and above:
julia> import Pkg; Pkg.add("https://github.com/anowacki/TauPy.jl")This package uses PyCall.jl package to
access ObsPy. If you have the default PyCall installation, then ObsPy will be
installed automatically via its own Conda environment. If you use your own
Python with PyCall, then you may need to install ObsPy for you installation
via conda, pip, or another means.
(Older versions of TauPy compatible with Julia v0.6 can be installed by doing
Pkg.clone("https://github.com/anowacki/TauPy.jl").)
If you receive and error like ERROR: InitError: Failed to import required Python module geographiclib when you first try using TauPy, then it's likely that
PyCall is set up to use your system python command, but the required packages
aren't installed or available. The easiest way to get things working is:
julia> ENV["PYTHON"] = ""; Pkg.build("PyCall")Restart, and then try to do using TauPy again. Note, however, that PyCall
will from hereon in always use its internal Conda python (which is at
PyCall.python).
TauPy exports three functions:
travel_timepathturning_depth
These take either epicentral distance, or source and receiver coordinates, and
return TauPy.Phase objects containing information about the phase. There are
two types exported by TauPy:
Phase, containing information about a seismic phase calculated using only event depth and epicentral distance; andPhaseGeog, which is the same but for source and receiver locations specified geographically (with longitude and latitude).
The interactive help describes the fields contained by Phases and PhaseGeogs.
To bring this up, type ?Phase or ?PhaseGeog and hit return.
The final positional argument of both travel_time and path is the name of
the seismic phase. This can be a single string, or an array of names. E.g.:
julia> using TauPy
julia> p = travel_time(0, 10, ["P", "PcP"])
2-element Array{Phase{Float64},1}:
Phase{Float64}("ak135", "P", 10.0, 0.0, 144.89570946391675, 13.700630345173362, 45.613198013389635, 45.613198013389635, Float64[], Float64[], Float64[])
Phase{Float64}("ak135", "PcP", 10.0, 0.0, 516.4444277972648, 0.9479529695834205, 2.834193976594543, 2.834193976594543, Float64[], Float64[], Float64[])
By default, all arrivals from a predetermined list are returned, corresponding
to the ‘phase’ "ttall".
With both travel_time and path, specify the Earth model by using the model
keyword argument like so:
julia> arr = travel_time(0, 10, "P", model="sp6")
1-element Array{Phase{Float64},1}:
Phase{Float64}("sp6", "P", 10.0, 0.0, 144.8972605261263, 13.7011317118041, 45.61534012667141, 45.61534012667141, Float64[], Float64[], Float64[])
Available models are listed by calling TauPy.available_models().
Use the travel_time function to quickly calculate the arrival times for
the triplicated arrivals at around 20° epicentral distance:
julia> using TauPy
julia> p = travel_time(110, 20, "P")
5-element Array{TauPy.Phase{Float64},1}:
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 263.806, 10.7956, 34.2707, 52.6707, Float64[], Float64[], Float64[])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 266.524, 11.5422, 37.0166, 58.2286, Float64[], Float64[], Float64[])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 266.525, 11.5214, 36.9391, 58.063, Float64[], Float64[], Float64[])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 267.698, 9.21572, 28.731, 42.7498, Float64[], Float64[], Float64[])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 268.261, 9.5515, 29.8818, 44.7109, Float64[], Float64[], Float64[])
julia> times = getfield.(p, :time)
5-element Array{Float64,1}:
263.80556674138126
266.52428738827155
266.5253559803772
267.6979484052481
268.26087550766334
Good luck in picking all of those!
You can also calculate the ray paths between the event and station:
julia> p = path(110, 20, "P")
5-element Array{TauPy.Phase{Float64},1}:
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 263.806, 10.7956, 34.2707, 52.6707, Float64[], [0.0, 0.120381, 0.142296, 0.1643, 0.208576, 0.298225, 0.472472, 0.652408, 0.667963, 0.683563 … 19.8216, 19.8493, 19.8632, 19.8701, 19.877, 19.9267, 19.9634, 19.9817, 19.9909, 20.0001], [6261.0, 6251.0, 6249.19, 6247.37, 6243.74, 6236.44, 6222.49, 6208.41, 6207.2, 6206.0 … 6343.5, 6347.25, 6349.13, 6350.06, 6351.0, 6359.05, 6365.02, 6368.01, 6369.51, 6371.0])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 266.524, 11.5422, 37.0166, 58.2286, Float64[], [0.0, 0.148344, 0.175382, 0.202552, 0.257294, 0.368423, 0.585572, 0.811497, 0.831113, 0.850797 … 19.8016, 19.8328, 19.8484, 19.8562, 19.864, 19.9189, 19.9596, 19.9799, 19.99, 20.0001], [6261.0, 6251.0, 6249.19, 6247.37, 6243.74, 6236.44, 6222.49, 6208.41, 6207.2, 6206.0 … 6343.5, 6347.25, 6349.13, 6350.06, 6351.0, 6359.05, 6365.02, 6368.01, 6369.51, 6371.0])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 266.525, 11.5214, 36.9391, 58.063, Float64[], [0.0, 0.147386, 0.174248, 0.20124, 0.255622, 0.366009, 0.581664, 0.805968, 0.82544, 0.844979 … 19.7997, 19.8308, 19.8463, 19.854, 19.8618, 19.9165, 19.9571, 19.9773, 19.9875, 19.9976], [6261.0, 6251.0, 6249.19, 6247.37, 6243.74, 6236.44, 6222.49, 6208.41, 6207.2, 6206.0 … 6343.5, 6347.25, 6349.13, 6350.06, 6351.0, 6359.05, 6365.02, 6368.01, 6369.51, 6371.0])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 267.698, 9.21572, 28.731, 42.7498, Float64[], [0.0, 0.0847973, 0.100216, 0.115685, 0.146774, 0.209565, 0.331002, 0.455546, 0.466272, 0.47702 … 19.858, 19.8798, 19.8907, 19.8961, 19.9016, 19.9415, 19.971, 19.9858, 19.9931, 20.0005], [6261.0, 6251.0, 6249.19, 6247.37, 6243.74, 6236.44, 6222.49, 6208.41, 6207.2, 6206.0 … 6343.5, 6347.25, 6349.13, 6350.06, 6351.0, 6359.05, 6365.02, 6368.01, 6369.51, 6371.0])
TauPy.Phase{Float64}("ak135", "P", 20.0, 110.0, 268.261, 9.5515, 29.8818, 44.7109, Float64[], [0.0, 0.0908228, 0.10734, 0.123913, 0.157226, 0.224534, 0.3548, 0.488528, 0.500051, 0.511599 … 19.8505, 19.8734, 19.8848, 19.8905, 19.8963, 19.9381, 19.9691, 19.9845, 19.9923, 20.0], [6261.0, 6251.0, 6249.19, 6247.37, 6243.74, 6236.44, 6222.49, 6208.41, 6207.2, 6206.0 … 6343.5, 6347.25, 6349.13, 6350.06, 6351.0, 6359.05, 6365.02, 6368.01, 6369.51, 6371.0])
If you want to know the arrivals’ turning depths, then turning_depth
is what you want:
julia> turning_depth.(p)
5-element Array{Float64,1}:
465.716
406.874
410.0
665.676
660.0
If you want to know the geographical coordinates of the path for an event and station, then use the source and receiver geographical coordinates:
julia> event_lon, event_lat, sta_lon, sta_lat, dep = 0, 0, 10, 10, 100;
julia> p = path(event_lon, event_lat, dep, sta_lon, sta_lat, "S")
3-element Array{TauPy.PhaseGeog{Float64},1}:
TauPy.PhaseGeog{Float64}("ak135", "S", 0.0, 0.0, 10.0, 10.0, 100.0, 14.106, 350.494, 24.1935, 48.835, 83.5499, Float64[], [0.0, 0.542449, 1.36541, 1.40396, 1.44327, 3.49073, 5.54737, 5.58703, 5.62594, 6.45806 … 9.8127, 9.83169, 9.84117, 9.85065, 9.88817, 9.9256, 9.96294, 9.98158, 9.99089, 10.0002], [0.0, 0.550792, 1.38607, 1.42518, 1.46506, 3.53789, 5.60619, 5.64589, 5.68482, 6.51558 … 9.81819, 9.83664, 9.84585, 9.85506, 9.89149, 9.92782, 9.96405, 9.98213, 9.99116, 10.0002], [6271.0, 6262.21, 6251.8, 6251.4, 6251.0, 6240.69, 6251.0, 6251.4, 6251.8, 6262.21 … 6347.25, 6349.13, 6350.06, 6351.0, 6356.0, 6361.0, 6366.0, 6368.5, 6369.75, 6371.0])
TauPy.PhaseGeog{Float64}("ak135", "S", 0.0, 0.0, 10.0, 10.0, 100.0, 14.106, 364.733, 20.4703, 39.5658, 57.2195, Float64[], [0.0, 0.0878754, 0.192757, 0.19682, 0.200884, 0.315745, 0.431873, 0.52082, 0.666153, 0.762522 … 9.87896, 9.89135, 9.89753, 9.90372, 9.93078, 9.9578, 9.98476, 9.99822, 10.0049, 10.0117], [0.0, 0.0892309, 0.195729, 0.199855, 0.203982, 0.320611, 0.438522, 0.528832, 0.676382, 0.774215 … 9.88255, 9.89457, 9.90058, 9.90658, 9.93285, 9.95906, 9.98522, 9.99827, 10.0048, 10.0113], [6271.0, 6262.21, 6251.8, 6251.4, 6251.0, 6239.72, 6228.43, 6219.86, 6206.0, 6196.9 … 6347.25, 6349.13, 6350.06, 6351.0, 6356.0, 6361.0, 6366.0, 6368.5, 6369.75, 6371.0])
TauPy.PhaseGeog{Float64}("ak135", "S", 0.0, 0.0, 10.0, 10.0, 100.0, 14.106, 364.855, 20.7161, 40.1367, 58.304, Float64[], [0.0, 0.0916586, 0.201103, 0.205344, 0.209587, 0.32951, 0.450821, 0.543782, 0.695759, 0.796594 … 9.86447, 9.87716, 9.8835, 9.88984, 9.91746, 9.94502, 9.97253, 9.98627, 9.99313, 9.99999], [0.0, 0.0930724, 0.204205, 0.208511, 0.212819, 0.334588, 0.457761, 0.552145, 0.706439, 0.808803 … 9.86848, 9.8808, 9.88695, 9.8931, 9.91991, 9.94667, 9.97336, 9.98668, 9.99334, 9.99999], [6271.0, 6262.21, 6251.8, 6251.4, 6251.0, 6239.72, 6228.43, 6219.86, 6206.0, 6196.9 … 6347.25, 6349.13, 6350.06, 6351.0, 6356.0, 6361.0, 6366.0, 6368.5, 6369.75, 6371.0])
Similarly, travel_time also accepts five geographic arguments if you know the
event and station coordinates.
Unfortunately, Obspy's travel time and raypath calculations are somewhat
slow. To speed up repeated calculations of the same times and raypaths,
TauPy implements a cache. To set the size of the cache, use the
TauPy.set_cache_size_mb!(size_mb) function. The cache can be clared
using TauPy.clear_cache!(). These functions are not exported.
To disable the cache for individual calls to path or travel_time,
pass the keyword argument cache=false.