Module rgrid

Raytracing on rectilinear grids

This module contains two classes to perform traveltime computation and raytracing on rectilinear grids:

  • Grid2d for 2D media

  • Grid3d for 3D media

Three algorithms are implemented

  • the Shortest-Path Method

  • the Fast-Sweeping Method

  • the Dynamic Shortest-Path Method

Slowness model can be defined in two ways:

  1. slowness constant within the voxels of the grid (the default)

  2. slowness defined at nodes of the grid

This code is part of ttcr ( https://github.com/groupeLIAMG/ttcr )

class ttcrpy.rgrid.Grid2d

class to perform raytracing with 2D rectilinear grids

x

node coordinates along x

Type:

np.ndarray

z

node coordinates along z

Type:

np.ndarray

dx

node separation along x

Type:

float

dz

node separation along z

Type:

float

shape

number of parameters along each dimension

Type:

(int, int)

nparams

total number of parameters for grid

Type:

int

n_threads

number of threads for raytracing

Type:

int

Constructor
Grid2d(x, z, n_threads=1, cell_slowness=1, method='SPM', aniso='iso', eps=1.e-15, maxit=20, weno=1, rotated_template=0, nsnx=10, nsnz=10, n_secondary=3, n_tertiary=3, radius_factor_tertiary=3.0, tt_from_rp=0) → Grid2d
Parameters:

  • x (np.ndarray) – node coordinates along x

  • z (np.ndarray) – node coordinates along z

  • n_threads (int) – number of threads for raytracing (default is 1)

  • cell_slowness (bool) – slowness defined for cells (True) or nodes (False) (default is 1)

  • method (string) –

    raytracing method (default is SPM)

    • ’FSM’ : fast marching method

    • ’SPM’ : shortest path method

    • ’DSPM’ : dynamic shortest path method

  • aniso (string) –

    type of anisotropy (implemented only for the SPM method)

    • ’iso’ : isotropic medium

    • ’elliptical’ : elliptical anisotropy

    • ’tilted_elliptical’ : tilted elliptical anisotropy

    • ’vti_psv’ : vertical transverse isotropy, P and SV waves

    • ’vti_sh’ : vertical transverse isotropy, SH waves

    • ’weakly_anelliptical’ : Weakly-Anelliptical formulation of B. Rommel

  • eps (double) – convergence criterion (FSM) (default is 1e-15)

  • maxit (int) – max number of sweeping iterations (FSM) (default is 20)

  • weno (bool) – use 3rd order weighted essentially non-oscillatory operator (FSM) (default is True)

  • rotated_template (bool) – use rotated templates (FSM)

  • nsnx (int) – number of secondary nodes in x (SPM) (default is 10)

  • nsnz (int) – number of secondary nodes in z (SPM) (default is 10)

  • n_secondary (int) – number of secondary nodes (DSPM) (default is 3)

  • n_tertiary (int) – number of tertiary nodes (DSPM) (default is 3)

  • radius_factor_tertiary (double) – multiplication factor used to compute radius of sphere around source that includes tertiary nodes (DSPM). The radius is the average edge length multiplied by this factor (default is 3)

  • tt_from_rp (bool) – compute traveltime using raypaths (available for FSM and DSPM only) (default is False)

Notes

For raytracing in anisotropic media, the convention for inputting slowness depends on the model. For elliptical anisotropy, the method set_slowness is used to input horizontal slowness, while for weakly anelliptical anisotropy, the method is used to input vertical slowness.

compute_D(coord)

Return matrix of interpolation weights for velocity data points constraint

Parameters:

coord (np.ndarray, shape (npts, 2)) – coordinates of data points

Returns:

D – Matrix of interpolation weights

Return type:

scipy csr_matrix, shape (npts, nparams)

Note

In the current implementation, no check is made to see if the coordinates are on a node, edge, or corner.

compute_K(order=1)

Compute smoothing matrices

Parameters:

order (int) – order of smoothing operator, accept 1 or 2 (1 by default)

Returns:

Kx, Kz – matrices for derivatives along x & z

Return type:

tuple of csr_matrix

static data_kernel_straight_rays(Tx, Rx, grx, grz, aniso=False) → L

Raytracing with straight rays in 2D

Parameters:

  • Tx (np.ndarray) –

    source coordinates, nTx by 2

    • 1st column contains X coordinates,

    • 2nd contains Z coordinates

  • Rx (np.ndarray) –

    receiver coordinates, nTx by 2

    • 1st column contains X coordinates,

    • 2nd contains Z coordinates

  • grx (np.ndarray) – grid node coordinates along x

  • grz (np.ndarray) – grid node coordinates along z

  • aniso (bool) – compute L for elliptically anisotropic medium (True) or isotropic medium (False)

Returns:

L – data kernel matrix (tt = L*slowness)

Return type:

scipy csr_matrix

Note

Tx and Rx should contain the same number of rows, each row corresponding to a source-receiver pair

dx

node separation along x

Type:

float

dz

node separation along x

Type:

float

get_grid_traveltimes(thread_no=0)

Obtain traveltimes computed at primary grid nodes

Parameters:

thread_no (int) – thread used to computed traveltimes (default is 0)

Returns:

tt

Return type:

np ndarray, shape (nx, nz)

get_number_of_cells()
Returns:

number of cells in grid

Return type:

int

get_number_of_nodes()
Returns:

number of nodes in grid

Return type:

int

get_s0(hypo, slowness=None)

Return slowness at source points

Parameters:

  • hypo (np.ndarray with 5 columns) –

    hypo holds source information, i.e.

    • 1st column is event ID number

    • 2nd column is origin time

    • 3rd column is source easting (X)

    • 4th column is source elevation (Z)

  • slowness (np ndarray, shape (nx, nz) (optional)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

get_slowness()

Returns slowness of grid

Returns:

slowness

Return type:

np ndarray, shape (nx, nz)

is_outside(pts)

Check if points are outside grid

Parameters:

pts (np ndarray, shape(npts, 3)) – coordinates of points to check

Returns:

True if at least one point outside grid

Return type:

bool

n_threads

number of threads for raytracing

Type:

int

nparams

total number of parameters for grid

Type:

int

raytrace(source, rcv, slowness=None, xi=None, theta=None, Vp0=None, Vs0=None, delta=None, epsilon=None, gamma=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) → tt, rays, L

Perform raytracing

Parameters:

  • source (2D np.ndarray with 2 or 3 columns) – see notes below

  • rcv (2D np.ndarray with 2 columns) – Columns correspond to x, y and z coordinates

  • slowness (np ndarray, shape (nx, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order) if None, slowness must have been assigned previously

  • xi (np ndarray, shape (nx, nz) (None by default)) – xi at grid cells (only for SPM & cell_slowness=True) xi may also have been flattened (with default ‘C’ order) if None, xi must have been assigned previously

  • theta (np ndarray, shape (nx, nz) (None by default)) – theta at grid cells (only for SPM & cell_slowness=True) theta may also have been flattened (with default ‘C’ order) if None, theta must have been assigned previously

  • Vp0 (np ndarray, shape (nx, nz) (None by default)) – Vp0 at grid cells (only for SPM & cell_slowness=True) Vp0 may also have been flattened (with default ‘C’ order) if None, Vp0 must have been assigned previously

  • Vs0 (np ndarray, shape (nx, nz) (None by default)) – Vs0 at grid cells (only for SPM & cell_slowness=True) Vs0 may also have been flattened (with default ‘C’ order) if None, Vs0 must have been assigned previously

  • delta (np ndarray, shape (nx, nz) (None by default)) – delta at grid cells (only for SPM & cell_slowness=True) delta may also have been flattened (with default ‘C’ order) if None, delta must have been assigned previously

  • epsilon (np ndarray, shape (nx, nz) (None by default)) – epsilon at grid cells (only for SPM & cell_slowness=True) epsilon may also have been flattened (with default ‘C’ order) if None, epsilon must have been assigned previously

  • gamma (np ndarray, shape (nx, nz) (None by default)) – gamma at grid cells (only for SPM & cell_slowness=True) gamma may also have been flattened (with default ‘C’ order) if None, gamma must have been assigned previously

  • thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor

  • aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event

  • compute_L (bool (False by default)) – Compute matrices of partial derivative of travel time w/r to slowness

  • return_rays (bool (False by default)) – Return raypaths

Returns:

  • tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)

  • rays (list of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)

  • L (scipy csr_matrix) – Matrix of partial derivative of travel time w/r to slowness

Notes

If source has 2 columns:

  • Columns correspond to x and z coordinates

  • Origin time (t0) is 0 for all points

If source has 3 columns:

  • 1st column corresponds to origin times

  • 2nd & 3rd columns correspond to x and z coordinates

source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.

set_Vp0(v)

Assign vertical Vp to grid (VTI medium)

Parameters:

v (np ndarray, shape (nx, nz)) – v may also have been flattened (with default ‘C’ order)

set_Vs0(v)

Assign vertical Vs to grid (VTI medium)

Parameters:

v (np ndarray, shape (nx, nz)) – v may also have been flattened (with default ‘C’ order)

set_delta(d)

Assign Thomsen delta parameter to grid (VTI medium, P-SV waves)

Parameters:

d (np ndarray, shape (nx, nz)) – d may also have been flattened (with default ‘C’ order)

set_epsilon(e)

Assign Thomsen epsilon parameter to grid (VTI medium, P-SV waves)

Parameters:

e (np ndarray, shape (nx, nz)) – e may also have been flattened (with default ‘C’ order)

set_gamma(g)

Assign Thomsen gamma parameter to grid (VTI medium, SH waves)

Parameters:

g (np ndarray, shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_s2(g)

Assign weakly anelliptical parameter s2

Parameters:

g (np ndarray, shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_s4(g)

Assign weakly anelliptical parameter s4

Parameters:

g (np ndarray, shape (nx, nz)) – g may also have been flattened (with default ‘C’ order)

set_slowness(slowness)

Assign slowness to grid

Parameters:

slowness (np ndarray, shape (nx, nz)) – slowness may also have been flattened (with default ‘C’ order)

set_tilt_angle(theta)

Assign tilted elliptical anisotropy angle to grid

Parameters:

theta (np ndarray, shape (nx, nz)) – theta may also have been flattened (with default ‘C’ order)

set_traveltime_from_raypath(ttrp)

Set option to compute traveltime using raypath

Parameters:

ttrp (bool) – option value

set_use_thread_pool(use_thread_pool)

Set option to use thread pool instead of parallel loop

Parameters:

use_thread_pool (bool) – option value

set_velocity(velocity)

Assign velocity to grid

Parameters:

velocity (np ndarray, shape (nx, nz)) – velocity may also have been flattened (with default ‘C’ order)

set_xi(xi)

Assign elliptical anisotropy ratio to grid

Parameters:

xi (np ndarray, shape (nx, nz)) – xi may also have been flattened (with default ‘C’ order)

shape

number of parameters along each dimension

Type:

list of int

to_vtk(fields, filename)

Save grid variables and/or raypaths to VTK format

Parameters:

  • fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the grid, or a list of raypath coordinates.

  • filename (str) – Name of file without extension for saving (extension vtr will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.

Notes

VTK files can be visualized with Paraview (https://www.paraview.org)

x

node coordinates along x

Type:

np.ndarray

z

node coordinates along z

Type:

np.ndarray

class ttcrpy.rgrid.Grid3d

class to perform raytracing with 3D rectilinear grids

x

node coordinates along x

Type:

np.ndarray

y

node coordinates along y

Type:

np.ndarray

z

node coordinates along z

Type:

np.ndarray

dx

node separation along x

Type:

float

dy

node separation along y

Type:

float

dz

node separation along z

Type:

float

shape

number of parameters along each dimension

Type:

(int, int, int)

nparams

total number of parameters for grid

Type:

int

n_threads

number of threads for raytracing

Type:

int

Constructor
Grid3d(x, y, z, n_threads=1, cell_slowness=1, method='FSM', tt_from_rp=1, interp_vel=0, eps=1.e-15, maxit=20, weno=1, nsnx=5, nsny=5, nsnz=5, n_secondary=2, n_tertiary=2, radius_factor_tertiary=3.0, translate_grid=False) → Grid3d
Parameters:

  • x (np.ndarray) – node coordinates along x

  • y (np.ndarray) – node coordinates along y

  • z (np.ndarray) – node coordinates along z

  • n_threads (int) – number of threads for raytracing (default is 1)

  • cell_slowness (bool) – slowness defined for cells (True) or nodes (False) (default is 1)

  • method (string) –

    raytracing method (default is FSM)

    • ’FSM’ : fast marching method

    • ’SPM’ : shortest path method

    • ’DSPM’ : dynamic shortest path

  • tt_from_rp (bool) – compute traveltimes from raypaths (FSM or DSPM only) (default is 1)

  • interp_vel (bool) – interpolate velocity instead of slowness at nodes (for cell_slowness == False or FSM) (defauls is False)

  • eps (double) – convergence criterion (FSM) (default is 1e-15)

  • maxit (int) – max number of sweeping iterations (FSM) (default is 20)

  • weno (bool) – use 3rd order weighted essentially non-oscillatory operator (FSM) (default is True)

  • nsnx (int) – number of secondary nodes in x (SPM) (default is 5)

  • nsny (int) – number of secondary nodes in y (SPM) (default is 5)

  • nsnz (int) – number of secondary nodes in z (SPM) (default is 5)

  • n_secondary (int) – number of secondary nodes (DSPM) (default is 2)

  • n_tertiary (int) – number of tertiary nodes (DSPM) (default is 2)

  • radius_factor_tertiary (double) – multiplication factor used to compute radius of sphere around source that includes tertiary nodes (DSPM). The radius is the average edge length multiplied by this factor (default is 3)

  • translate_grid (bool) – Translate the grid such that origin is (0, 0, 0) to perform computations, which may increase accuracy when large values, e.g. UTM coordinates, are used. When raytracing, src and rcv should be given in the original system, and output raypath coordinates are also given in the original system (default if False)

static builder(filename, n_threads=1, method='FSM', tt_from_rp=1, interp_vel=0, eps=1.e-15, maxit=20, weno=1, nsnx=5, nsny=5, nsnz=5, n_secondary=2, n_tertiary=2, radius_factor_tertiary=3.0, translate_grid=0)

Build instance of Grid3d from VTK file

Parameters:

  • filename (str) – Name of file holding a vtkRectilinearGrid. The grid must have point or cell attribute named either ‘Slowness’, ‘slowness’, ‘Velocity’, ‘velocity’, or ‘P-wave velocity’

  • Constructor (Other parameters are defined in) –

Returns:

grid – grid instance

Return type:

Grid3d

compute_D(coord)

Return matrix of interpolation weights for velocity data points constraint

Parameters:

coord (np.ndarray, shape (npts, 3)) – coordinates of data points

Returns:

D – Matrix of interpolation weights

Return type:

scipy csr_matrix, shape (npts, nparams)

Note

In the current implementation, no check is made to see if the coordinates are on a node, edge, face, or corner.

compute_K()

Compute smoothing matrices (2nd order derivative)

Returns:

Kx, Ky, Kz – matrices for derivatives along x, y, & z

Return type:

tuple of csr_matrix

static data_kernel_straight_rays(Tx, Rx, grx, gry, grz, centers) -> L, (xc, yc, zc)

Raytracing with straight rays in 3D

Parameters:

  • Tx (np.ndarray) –

    source coordinates, nTx by 3

    • 1st column contains X coordinates,

    • 2nd contains Y coordinates

    • 3rd contains Z coordinates

  • Rx (np.ndarray) –

    receiver coordinates, nTx by 3

    • 1st column contains X coordinates,

    • 2nd contains Y coordinates

    • 3rd contains Z coordinates

  • grx (np.ndarray) – grid node coordinates along x

  • gry (np.ndarray) – grid node coordinates along y

  • grz (np.ndarray) – grid node coordinates along z

  • centers (bool) – return coordinates of center of cells (False by default)

Returns:

  • L (scipy csr_matrix) – data kernel matrix (tt = L*slowness)

  • (xc, yc, zc) (tuple of np.ndarray) – vectors of coordinates of center of cells

Note

Tx and Rx should contain the same number of rows, each row corresponding to a source-receiver pair

dx

node separation along x

Type:

float

dy

node separation along y

Type:

float

dz

node separation along z

Type:

float

get_grid_traveltimes(thread_no=0)

Obtain traveltimes computed at primary grid nodes

Parameters:

thread_no (int) – thread used to computed traveltimes (default is 0)

Returns:

tt – traveltimes

Return type:

np ndarray, shape (nx, ny, nz)

get_number_of_cells()
Returns:

number of cells in grid

Return type:

int

get_number_of_nodes()
Returns:

number of nodes in grid

Return type:

int

get_s0(hypo, slowness=None)

Return slowness at source points

Parameters:

  • hypo (np.ndarray with 5 columns) –

    hypo holds source information, i.e.

    • 1st column is event ID number

    • 2nd column is origin time

    • 3rd column is source easting

    • 4th column is source northing

    • 5th column is source elevation

  • slowness (np ndarray, shape (nx, ny, nz) (optional)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order)

Returns:

s0 – slowness at source points

Return type:

np.ndarray

get_slowness()

Returns slowness of grid

Returns:

slowness

Return type:

np ndarray, shape (nx, ny, nz)

ind(i, j, k)

Return node index

Parameters:

  • i (int) – index of node along x

  • j (int) – index of node along y

  • k (int) – index of node along z

Returns:

node index for a “flattened” grid

Return type:

int

indc(i, j, k)

return cell index

Parameters:

  • i (int) – index of cell along x

  • j (int) – index of cell along y

  • k (int) – index of cell along z

Returns:

cell index for a “flattened” grid

Return type:

int

is_outside(pts)

Check if points are outside grid

Parameters:

pts (np ndarray, shape(npts, 3)) – coordinates of points to check

Returns:

True if at least one point outside grid

Return type:

bool

n_threads

number of threads for raytracing

Type:

int

nparams

total number of parameters for grid

Type:

int

raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, compute_M=False, return_rays=False)
raytrace(source, rcv, slowness=None, thread_no=None,

aggregate_src=False, compute_L=False, compute_M=False, return_rays=False) -> tt, rays, M, L

Perform raytracing

Parameters:

  • source (2D np.ndarray with 3, 4 or 5 columns) – see notes below

  • rcv (2D np.ndarray with 3 columns) – Columns correspond to x, y and z coordinates

  • slowness (np ndarray, shape (nx, ny, nz) (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) slowness may also have been flattened (with default ‘C’ order) if None, slowness must have been assigned previously

  • thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor

  • aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event

  • compute_L (bool (False by default)) –

    Compute matrices of partial derivative of travel time w/r to slowness (implemeted for the SPM & DSPM with slowness

    defined at cells).

  • compute_M (bool (False by default)) – Compute matrices of partial derivative of travel time w/r to velocity Note : compute_M and compute_L are mutually exclusive

  • return_rays (bool (False by default)) – Return raypaths

Returns:

  • tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)

  • rays (list of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)

  • M (list of csr_matrix) – matrices of partial derivative of travel time w/r to velocity. the number of matrices is equal to the number of sources

  • L (scipy csr_matrix) – Matrix of partial derivative of travel time w/r to slowness. if input argument source has 5 columns, L is a list of matrices and the number of matrices is equal to the number of sources otherwise, L is a single csr_matrix

Notes

If source has 3 columns:

  • Columns correspond to x, y and z coordinates

  • Origin time (t0) is 0 for all points

If source has 4 columns:

  • 1st column corresponds to origin times

  • 2nd, 3rd & 4th columns correspond to x, y and z coordinates

If source has 5 columns:

  • 1st column corresponds to event ID

  • 2nd column corresponds to origin times

  • 3rd, 4th & 5th columns correspond to x, y and z coordinates

For the latter case (5 columns), source and rcv should contain the same number of rows, each row corresponding to a source-receiver pair. For the 2 other cases, source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.

set_slowness(slowness)

Assign slowness to grid

Parameters:

slowness (np ndarray, shape (nx, ny, nz)) – slowness may also have been flattened (with default ‘C’ order)

set_traveltime_from_raypath(ttrp)

Set option to compute traveltime using raypath

Parameters:

ttrp (bool) – option value

set_use_thread_pool(use_thread_pool)

Set option to use thread pool instead of parallel loop

Parameters:

use_thread_pool (bool) – option value

set_velocity(velocity)

Assign velocity to grid

Parameters:

velocity (np ndarray, shape (nx, ny, nz)) – velocity may also have been flattened (with default ‘C’ order)

shape

number of parameters along each dimension

Type:

list of int

to_vtk(fields, filename)

Save grid variables and/or raypaths to VTK format

Parameters:

  • fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the grid, or a list of raypath coordinates.

  • filename (str) – Name of file without extension for saving (extension vtr will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.

Notes

VTK files can be visualized with Paraview (https://www.paraview.org)

x

node coordinates along x

Type:

np.ndarray

y

node coordinates along y

Type:

np.ndarray

z

node coordinates along z

Type:

np.ndarray

ttcrpy.rgrid.set_verbose(v)

Set verbosity level for C++ code

Parameters:

v (int) – verbosity level