TauPy

Calculate properties of teleseismic arrivals through a selection of 1D Earth models, using the ObsPy Python software.

Install

To install on Julia versions v1.6 and above: julia 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.

Problems importing geographiclib or obspy.taup

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 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).

Use

TauPy exports three functions:

• travel_time
• path
• turning_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; and
• PhaseGeog, 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.

Specifying the seismic phase

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 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".

Specifying the Earth model

With both travel_time and path, specify the Earth model by using the model keyword argument like so:

```julia 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().

Examples

Use the travel_time function to quickly calculate the arrival times for the triplicated arrivals at around 20° epicentral distance:

```julia 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 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 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 julia> eventlon, eventlat, stalon, stalat, dep = 0, 0, 10, 10, 100;

julia> p = path(eventlon, eventlat, dep, stalon, stalat, "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.

Calculation cache

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 cleared 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.

Related packages

Global 1D traveltimes and raytracing are also available in the following alternative packages:

• TauP.jl wraps Phil Crotwell's eminent TauP Java suite of programs under a LGPL-3 licence and using JavaCall.jl to directly call the TauP library.