Wrong ray paths and arrival time using TauP on custom .nd model. Bug in TauPyModel.get_ray_paths and arrivals.plot_rays

[edur409]

Avoid duplicates

• I searched existing issues

Bug Summary

Hi All,

The folowing .nd model was created to do ray tracing on a fruit, therefore the unusual values for P and S velocities.

0.00 200.00 141.42 1.30
298.38 200.00 141.42 1.30
596.75 200.00 141.42 1.30
895.12 200.00 141.42 1.30
1193.50 200.00 141.42 1.30
mantle
1193.50 150.00 106.07 1.30
1491.82 150.00 106.07 1.30
1790.13 150.00 106.07 1.30
2088.45 150.00 106.07 1.30
2386.76 150.00 106.07 1.30
outer-core
2386.76 150.15 106.17 1.30
2386.82 150.15 106.17 1.30
2386.87 150.15 106.17 1.30
2386.92 150.15 106.17 1.30
2386.98 150.15 106.17 1.30
inner-core
2386.98 150.30 106.28 1.30
2386.98 150.30 106.28 1.30
2386.99 150.30 106.28 1.30
2386.99 150.30 106.28 1.30
2387.00 150.30 106.28 1.30

The ray tracing produces this result at epicentral distance 120.0 degrees and a P arrival time of 20.672 seconds is produced by the get_ray_paths function:

The TauP software in Java correctly reports no arrivals for that epicentral distance:

Hope this information is enough to reproduce the results. Thank you so much for keeping TauP as part of the Obspy package.

Regards,
Evert

Code to Reproduce

import numpy as np
import matplotlib.pyplot as plt
from obspy.taup import TauPyModel
from obspy.taup.taup_create import build_taup_model

# Folder where .nd file (velocity model) is stored
filename = ndfolder+'/'+name+'.nd'

# Create the .npz model file from the *.nd file
build_taup_model(filename,output_folder=ndfolder,verbose=True) 

# Load the model into TauP    
model = TauPyModel(model=f"{ndfolder}/{name}.npz") #

sourcedepth = 0.0 # Source depth
angle_test = 120.0 # Epicentral distance to check

phaselist = ['P','PKP','PKIKP']

arrivals = model.get_ray_paths(sourcedepth, angle_test, phase_list=phaselist)

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax = arrivals.plot_rays(plot_type='spherical', phase_list=phaselist,
                   legend=True, fig = fig, ax = ax)
fig.savefig('Obspy_Raypaths_'+name+'_'+phaselist[0]+'.png', dpi=200)
plt.show()

# See all the arrivals.
print(arrivals)

Error Traceback

No response

ObsPy Version?

1.4.0

Operating System?

Ubuntu

Python Version?

3.10.8

Installation Method?

conda

[megies]

@johnrudge maybe you have an idea about this one?

[johnrudge]

I can confirm there's a bug here, but I don't have a quick fix. I'll have a think.

[johnrudge]

I think this is some kind of rounding issue. The 120 degrees epicentral distance is rather special for the given model -- it's the beginning of the P shadow zone. The behaviour in obspy.taup is correct either side of this special value of 120 degrees.

arrivals = model.get_travel_times(0.0, 119.99999, phase_list=['P'])
print(arrivals)

1 arrivals
	P phase arrival at 20.672 seconds

arrivals = model.get_travel_times(0.0, 120.00001, phase_list=['P'])
print(arrivals)

0 arrivals