Avoid duplicates
- [X] I searched existing issues
Bug Summary
Hi Again,
this bug/feature I discussed a few months ago with Philip Crotwell @crotwell , and thought of reporting it here as well.
The issue: No PKIKP phase was calculated at 180 degrees for my model (outer core Vs >0) [see .nd model below].
Philip Crotwell kindly identified that a workaround was setting the Vs for the outer core as zero (liquid). I was interested in knowing the direct transmission time of the P wave through all the layers using a solid outer core (Vs>0); since I am interested only in P waves I set my outer core Vs to zero and the PKIKP phase worked. I leave it here so that a solution is implemented in the future or someone dealing with a similar problem has this workaround in mind.
Model with solid outer core:
0.000 60.000 42.426 1.300
225.000 60.000 42.426 1.300
450.000 60.000 42.426 1.300
675.000 60.000 42.426 1.300
900.000 60.000 42.426 1.300
mantle
900.000 70.001 49.497 1.300
1050.000 70.001 49.497 1.300
1200.000 70.001 49.497 1.300
1350.000 70.001 49.497 1.300
1500.000 70.001 49.497 1.300
outer-core
1500.000 80.001 56.569 1.300
1625.000 80.001 56.569 1.300
1750.000 80.001 56.569 1.300
1875.000 80.001 56.569 1.300
2000.000 80.001 56.569 1.300
inner-core
2000.000 99.999 70.711 1.300
2125.000 99.999 70.711 1.300
2250.000 99.999 70.711 1.300
2375.000 99.999 70.711 1.300
2500.000 99.999 70.711 1.300
Result in Obspy 1.4.0 for phase 'PKIKP' :
model.get_travel_times: 0 arrivals
model.get_ray_paths: 0 arrivals
Result in TauP 2.6.0 for phase 'PKIKP' :
'''Model: 4Layer_model
Distance Depth Phase Travel Ray Param Takeoff Incident Purist Purist
(deg) (km) Name Time (s) p (s/deg) (deg) (deg) Distance Name
-----------------------------------------------------------------------------------'''
Obspy 1.4.0 produces the right answer for phase 'PKP' though:
model.get_travel_times: 1 arrivals
PKP phase arrival at 69.643 seconds
model.get_ray_paths: 1 arrivals
PKP phase arrival at 69.643 seconds
Now, if you set the outer core Vs to zero, the results are as follow:
Result in Obspy 1.4.0 for phase 'PKIKP' (Vs = 0 in outer core) :
model.get_travel_times: 1 arrivals
PKIKP phase arrival at 69.643 seconds
model.get_ray_paths: 1 arrivals
PKIKP phase arrival at 69.643 seconds
And no 'PKP' for this model:
model.get_travel_times: 0 arrivals
model.get_ray_paths: 0 arrivals
Result in TauP 2.6.0 for phase 'PKIKP' (Vs = 0 in outer core) :
'''Model: 4Layer_model_OC_VS0
Distance Depth Phase Travel Ray Param Takeoff Incident Purist Purist
(deg) (km) Name Time (s) p (s/deg) (deg) (deg) Distance Name
-----------------------------------------------------------------------------------'''
180.00 0.0 PKIKP 69.64 0.000 0.00 0.00 180.00 = PKIKP'''
Thank you all for your wonderful efforts to improve these tools for seismology and beyond,
Regards,
Evert
Code to Reproduce
# 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 = 180 # Epicentral distance to check
phaselist = ['PKIKP']
arrivals2 = model.get_travel_times(sourcedepth, angle_test, phase_list=phaselist)
print('model.get_travel_times: ',arrivals2)
arrivals = model.get_ray_paths(sourcedepth, angle_test, phase_list=phaselist)
# See all the arrivals.
print('model.get_ray_paths: ',arrivals)
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()
Error Traceback
Building obspy.taup model for '/home/leo/Documents/Obspy_BUG/nd_models/4Layer_model.nd' ...
filename = /home/leo/Documents/Obspy_BUG/nd_models/4Layer_model.nd
Done reading velocity model.
Radius of model . is 2500.0
Using parameters provided in TauP_config.ini (or defaults if not) to call SlownessModel...
Parameters are:
taup.create.min_delta_p = 0.1 sec / radian
taup.create.max_delta_p = 11.0 sec / radian
taup.create.max_depth_interval = 115.0 kilometers
taup.create.max_range_interval = 0.04363323129985824 degrees
taup.create.max_interp_error = 0.05 seconds
taup.create.allow_inner_core_s = True
Slow model 235 P layers,280 S layers
Done calculating Tau branches.
Done Saving /home/leo/Documents/Obspy_BUG/nd_models/4Layer_model.npz
Method run is done, but not necessarily successful.
model.get_travel_times: 0 arrivals
model.get_ray_paths: 0 arrivals
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[103], line 23
20 print('model.get_ray_paths: ',arrivals)
22 fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
---> 23 ax = arrivals.plot_rays(plot_type='spherical', phase_list=phaselist,
24 legend=True, fig = fig, ax = ax)
25 fig.savefig('Obspy_Raypaths_'+name+'_'+phaselist[0]+'.png', dpi=200)
26 plt.show()
File ~/anaconda3/envs/obspy2/lib/python3.10/site-packages/obspy/taup/tau.py:331, in Arrivals.plot_rays(self, phase_list, plot_type, plot_all, legend, label_arrivals, show, fig, ax, indicate_wave_type)
328 arrivals.append(arrival)
330 if not arrivals:
--> 331 raise ValueError("Can only plot arrivals with calculated ray "
332 "paths.")
334 # get the velocity discontinuities in your model, for plotting:
335 discons = self.model.s_mod.v_mod.get_discontinuity_depths()
ValueError: Can only plot arrivals with calculated ray paths.
ObsPy Version?
1.4.0
Operating System?
Ubuntu
Python Version?
3.10.8
Installation Method?
conda
bug-unconfirmed