How to interpret index of traveltime.values in sphere coordinates?

[pan3rock]

I thought the index 0 is R, index 1 is 90 - latitude and index 2 is longitude. However, I run the Spherical 2D example and then solver.traveltime.values[0, 1, :] show different values. I can't understand it the source is put at [0, 0, 0] (the North Pole), why are the values different for the same latitude? The values in the figure of that page seems the same. Thanks!

The code is below:

import numpy as np
import pykonal
import matplotlib.pyplot as plt

solver = pykonal.EikonalSolver(coord_sys="spherical")
solver.velocity.min_coords = 6371., 0, 0
solver.velocity.node_intervals = 1, np.pi/32, np.pi/32
solver.velocity.npts = 1, 33, 64
solver.velocity.values = np.ones(solver.velocity.npts)

src_idx = (0, 0, 0)
solver.traveltime.values[src_idx] = 0
solver.unknown[src_idx] = False
solver.trial.push(*src_idx)
solver.solve()

i1 = np.arange(33) * np.degrees(np.pi/32)
i2 = np.arange(64) * np.degrees(np.pi/32)
fig, ax = plt.subplots()
pc = ax.pcolormesh(i2, 90-i1, solver.traveltime.values[0, :, :], cmap='jet')
ax.set_xlabel("longitude")
ax.set_ylabel("latitude")
fig.colorbar(pc, ax=ax)
plt.tight_layout()
plt.show()

solver.traveltime.value shows the same conclusion: different values are obtained for the same longitude. Because the velocities are always 1, the value with longtitude 0 is correct (111.19) but others are wrong.

from pykonal.transformations import geo2sph, sph2xyz

def print_tt(lat, lon):
    loc = np.array([lat, lon, 0])
    loc_sph = geo2sph(loc)
    print(solver.traveltime.value(loc_sph))

print_tt(89, 0)
print_tt(89, 10)

outputs:

111.19492664455889
1248.7380853952743

[malcolmw]

Hi @pan3rock, it looks like you uncovered a bug in the way that pykonal identifies neighbouring grid nodes near the pole... In the example above, there are actually 64 nodes at the north pole (i.e., all nodes with $\rho=\theta=0$) but only those with i_phi in [1, 63] are considered neighbours to that with i_phi = 0. Furthermore, the distance between these nodes is zero, which probably creates a divide by zero error I never handled properly.

If you set the traveltime to zero for all the nodes at the north pole, the example should work as expected.

for i_phi in range(solver.velocity.npts[2]):
    src_idx = (0, 0, i_phi)
    solver.traveltime.values[src_idx] = 0
    solver.unknown[src_idx] = False
    solver.trial.push(*src_idx)
solver.solve()

The above work around is obviously a hack and should be resolved properly within the pykonal package but I don't have time to work on this right now.

[pan3rock]

You're right. After the above corrections, the value of traveltime matches distances along great circle. Thanks a lot.