malcolmw / pykonal

Travel-time calculator based on the fast-marching method solution to the Eikonal equation.
https://malcolmw.github.io/pykonal-docs/
GNU General Public License v3.0
154 stars 55 forks source link

Using Shear Wave Velocity Model #32

Closed snoltcaraway closed 1 year ago

snoltcaraway commented 1 year ago

I am using this package to try to trace ray paths from a shear wave velocity model, which has shear wave velocity reductions up to 30%, similar to the PyKonal_rays image. When I import the velocity model, I interpolate it based on Vs, and get a matrix that plots how it should look. However, when I use the EikonalSolver (cartesian), the travel times and resulting ray paths do not seem correct. I have tried different source and receiver locations. I am tracing the rays using solver.trace_ray().

Do you have source code or an example on how you got the pykonal_rays figure (figure on the welcome page).

malcolmw commented 1 year ago

Hi @snoltcaraway, try the example below.

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

solver = pykonal.solver.PointSourceSolver(coord_sys='cartesian')
solver.vv.min_coords = 0, 0, 0
solver.vv.node_intervals = 1, 1, 1
solver.vv.npts = 128, 32, 1
solver.vv.values = 4.5 * np.ones(solver.vv.npts)
solver.vv.values[20:44, 8:24] = 2
solver.vv.values[84:108, 8:24] = 9
solver.vv.values += 5*np.random.randn(*solver.vv.npts)
solver.vv.values = scipy.ndimage.gaussian_filter(solver.vv.values, 2)

solver.src_loc = 64, 24, 0
solver.solve()

plt.close('all')
fig, ax = plt.subplots()
qmesh = ax.pcolormesh(
    solver.vv.nodes[:, :, 0, 0], 
    solver.vv.nodes[:, :, 0, 1], 
    solver.vv.values[:, :, 0],
    shading='gouraud',
    cmap=plt.get_cmap('nipy_spectral_r')
)
for rx in np.arange(0, 129, 8):
    ray = solver.trace_ray(np.array([rx, 0.1, 0.]))
    ax.plot(ray[:, 0], ray[:, 1], linewidth=1, color='k')

ax.scatter(solver.src_loc[0], solver.src_loc[1], marker='*', edgecolor='k', facecolor='white', s=1024, zorder=100)

ax.set_aspect(1)
ax.invert_yaxis()
fig.colorbar(qmesh, ax=ax, orientation='horizontal')
snoltcaraway commented 1 year ago

That worked.

I see that you are using a the PointSourceSolver in a heterogenous medium instead of the EikonalSolver. Is there a benefit of choosing one or the other? Eikonal is supposed to computer the shortest travel time correct?

Thank you so much!

malcolmw commented 1 year ago

Great!

The PointSourceSolver class implements a pair of computational grids to improve accuracy, particularly in the near-source region where wavefront curvature is high (see. And yes, both the EikonalSolver and PointSourceSolver compute the first arrival. This is illustrated in White et al. (2020) Figure 10a,c.

malcolmw commented 1 year ago

@snoltcaraway I'm closing this issue, but feel free to reach out if you have any further questions.

Cheers, Malcolm