Problem with computing ray paths

[tianzeliu]

Hi There, I found that solver.trace_ray() fails to output ray paths in some parts of my velocity model. In those cases, it only gives the end points of the rays. Do you know why is that? Thanks!

Tianze

[malcolmw]

Hi, @tianzeliu, Can you post a minimum working example (with any necessary data files) that reproduces your error?

- Malcolm

[tianzeliu]

Hi Malcom, Below is a minimum working example of my code. It reads a NetCDF-format velocity-model file, which is attached to this message. The model represents a low-velocity fault zone sandwiched between normal crustal rocks. The source is located outside the fault zone at the surface, and the receiver is located inside the fault zone at depth.

Tianze

from netCDF4 import Dataset
import numpy as np
import pykonal

## Inputs
# Velocity model
path_mod = './VerticalFaultZoneLayeredWallRock_Vp5.00_Vs2.89_HomoCrust.nc'

# Source location
x_src = 25.699493367599036
y_src = 9.441364007690407

# Receiver location
x_rec = -30.
y_rec = -2.
dep_rec = 5.

## Read the input velocity model
print('Reading the 3D velocity model...')
rootgrp = Dataset(path_mod)

vect_x = rootgrp.variables['x'][:]
vect_y = rootgrp.variables['y'][:]
vect_dep = rootgrp.variables['depth'][:]
array_vp = rootgrp.variables['vp'][:]
array_vs = rootgrp.variables['vs'][:]

vect_z = np.amax(vect_dep)-np.flip(vect_dep)
array_vp = np.flip(array_vp, axis=2)
array_vs = np.flip(array_vs, axis=2)

rootgrp.close()

dx = vect_x[1]-vect_x[0]
dy = vect_y[1]-vect_y[0]
dz = vect_z[1]-vect_z[0]

x_min = np.amin(vect_x)
x_max = np.amax(vect_x)

y_min = np.amin(vect_y)
y_max = np.amax(vect_y)

z_min = np.amin(vect_z)
z_max = np.amax(vect_z)

xnum = len(vect_x)
ynum = len(vect_y)
znum = len(vect_z)

## Compute the P travel time
xind_src = round((x_src-x_min)/dx)
yind_src = round((y_src-y_min)/dy)
zind_src = znum-1 # Source is at the surface

print('Computing the P travel times...')
solver = pykonal.EikonalSolver(coord_sys="cartesian")
solver.velocity.min_coords = [0, 0, 0]
solver.velocity.node_intervals = [dx, dy, dz]
solver.velocity.npts = [xnum, ynum, znum]
solver.velocity.values = array_vp

src_ind = (xind_src, yind_src, zind_src)
solver.traveltime.values[src_ind] = 0
solver.unknown[src_ind] = False
solver.trial.push(*src_ind)

solver.solve()

array_tp = solver.traveltime.values

## Find the ray path
print('Computing the P ray paths...')
z_rec =  z_max-dep_rec
ray = solver.trace_ray(np.array([x_rec, y_rec, z_rec]))
print(ray)

VerticalFaultZoneLayeredWallRock_Vp5.00_Vs2.89_HomoCrust.nc.zip

[malcolmw]

It looks to me like the receiver location is outside of your computational grid. The minimum x coordinate of the grid is 0, whereas the x coordinate of the receiver is -30. You will have to enlarge your computational grid to include the receiver location.

You might also find the PointSourceSolver useful (https://malcolmw.github.io/pykonal-docs/examples/point_source_solver.html). It uses a refined grid in the source region to improve accuracy and it makes specifying the source location easier (the source does not need to be fixed to a grid point). To use it, change

solver = pykonal.EikonalSolver(coord_sys="cartesian")

to

solver = pykonal.solver.PointSourceSolver(coord_sys="cartesian")

and

src_ind = (xind_src, yind_src, zind_src)
solver.traveltime.values[src_ind] = 0
solver.unknown[src_ind] = False
solver.trial.push(*src_ind)

to

solver.src_loc = x_src, y_src, 0

[tianzeliu]

Ah, I see. I forgot to change the lower bounds of the axes. Thanks a lot!

[malcolmw]

No problem! I'm going to close this issue, but let me know if you have any other problems.

- Malcolm