# Source code for pylops.utils.seismicevents

import numpy as np
import scipy.signal as filt

def _filterdata(d, nt, wav, wcenter):
r"""Apply filtering to data with wavelet wav"""
dwav = filt.lfilter(wav, 1, d, axis=-1)
dwav = dwav[..., wcenter:]
d = d[..., : (nt - wcenter)]
return d, dwav

[docs]def makeaxis(par):
r"""Create axes t, x, and y axes

Create space and time axes from dictionary containing initial values ot, ox, oy,
sampling steps dt, dx, dy and number of elements nt, nx, ny
for each axis.

Parameters
----------
par : :obj:dict
Dictionary containing initial values, sampling steps, and number of elements

Returns
-------
t : :obj:numpy.ndarray
Time axis
t2 : :obj:numpy.ndarray
Symmetric time axis
x : :obj:numpy.ndarray
x axis
y : :obj:numpy.ndarray
y axis (None, if oy, dy or ny are not provided)

Examples
--------
>>> par = {'ox':0, 'dx':2, 'nx':60,
>>>        'oy':0, 'dy':2, 'ny':100,
>>>        'ot':0, 'dt':4, 'nt':400}
>>> # Create axis
>>> t, t2, x, y = makeaxis(par)
"""
x = par["ox"] + np.arange(par["nx"]) * par["dx"]
t = par["ot"] + np.arange(par["nt"]) * par["dt"]
t2 = np.arange(-par["nt"] + 1, par["nt"]) * par["dt"]

if "oy" in par.keys():
y = par["oy"] + np.arange(par["ny"]) * par["dy"]
else:
y = None
return t, t2, x, y

[docs]def linear2d(x, t, v, t0, theta, amp, wav):
r"""Linear 2D events

Create 2d linear events given propagation velocity, intercept time, angle,
and amplitude of each event

Parameters
----------
x : :obj:numpy.ndarray
space axis
t : :obj:numpy.ndarray
time axis
v : :obj:float
propagation velocity
t0 : :obj:tuple or :obj:float
intercept time at :math:x=0 of each linear event
theta : :obj:tuple or :obj:float
angle (in degrees) of each linear event
amp : :obj:tuple or :obj:float
amplitude of each linear event
wav : :obj:numpy.ndarray
wavelet to be applied to data

Returns
-------
d : :obj:numpy.ndarray
data without wavelet  of size
:math:[n_x \times n_t]
dwav : :obj:numpy.ndarray
data with wavelet  of size
:math:[n_x \times n_t]

Notes
-----
Each event is created using the following relation:

.. math::
t_i(x) = t_{0,i} + p_{x,i} x

where :math:p_{x,i}=\sin( \theta_i)/v

"""
if isinstance(t0, (float, int)):
t0 = (t0,)
if isinstance(theta, (float, int)):
theta = (theta,)
if isinstance(amp, (float, int)):
amp = (amp,)

# identify dimensions
dt = t[1] - t[0]
wcenter = int(len(wav) / 2)
nx = np.size(x)
nt = np.size(t) + wcenter
nevents = np.size(t0)

# create events
d = np.zeros((nx, nt))
for ievent in range(nevents):
px = np.sin(np.deg2rad(theta[ievent])) / v
tevent = t0[ievent] + px * x
tevent = (tevent - t[0]) / dt
itevent = tevent.astype(int)
dtevent = tevent - itevent
for ix in range(nx):
if itevent[ix] < nt - 1 and itevent[ix] >= 0:
d[ix, itevent[ix]] += amp[ievent] * (1 - dtevent[ix])
d[ix, itevent[ix] + 1] += amp[ievent] * dtevent[ix]

# filter events with certain wavelet
d, dwav = _filterdata(d, nt, wav, wcenter)
return d, dwav

[docs]def parabolic2d(x, t, t0, px, pxx, amp, wav):
r"""Parabolic 2D events

Create 2d parabolic events given intercept time,
slowness, curvature, and amplitude of each event

Parameters
----------
x : :obj:numpy.ndarray
space axis
t : :obj:numpy.ndarray
time axis
t0 : :obj:tuple or :obj:float
intercept time at :math:x=0 of each parabolic event
px : :obj:tuple or :obj:float
slowness of each parabolic event
pxx : :obj:tuple or :obj:float
curvature of each parabolic event
amp : :obj:tuple or :obj:float
amplitude of each parabolic event
wav : :obj:numpy.ndarray
wavelet to be applied to data

Returns
-------
d : :obj:numpy.ndarray
data without wavelet of size
:math:[n_x \times n_t]
dwav : :obj:numpy.ndarray
data with wavelet of size
:math:[n_x \times n_t]

Notes
-----
Each event is created using the following relation:

.. math::
t_i(x) = t_{0,i} + p_{x,i} x + p_{xx,i} x^2

"""
if isinstance(t0, (float, int)):
t0 = (t0,)
if isinstance(px, (float, int)):
px = (px,)
if isinstance(pxx, (float, int)):
pxx = (pxx,)
if isinstance(amp, (float, int)):
amp = (amp,)

# identify dimensions
dt = t[1] - t[0]
wcenter = int(len(wav) / 2)
nx = np.size(x)
nt = np.size(t) + wcenter
nevents = np.size(t0)

# create events
d = np.zeros((nx, nt))
for ievent in range(nevents):
tevent = t0[ievent] + px[ievent] * x + pxx[ievent] * x ** 2
tevent = (tevent - t[0]) / dt
itevent = tevent.astype(int)
dtevent = tevent - itevent
for ix in range(nx):
if itevent[ix] < nt - 1 and itevent[ix] >= 0:
d[ix, itevent[ix]] += amp[ievent] * (1 - dtevent[ix])
d[ix, itevent[ix] + 1] += amp[ievent] * dtevent[ix]

# filter events with certain wavelet
d, dwav = _filterdata(d, nt, wav, wcenter)
return d, dwav

[docs]def hyperbolic2d(x, t, t0, vrms, amp, wav):
r"""Hyperbolic 2D events

Create 2d hyperbolic events given intercept time, root-mean-square
velocity, and amplitude of each event

Parameters
----------
x : :obj:numpy.ndarray
space axis
t : :obj:numpy.ndarray
time axis
t0 : :obj:tuple or :obj:float
intercept time at :math:x=0 of each of hyperbolic event
vrms : :obj:tuple or :obj:float
root-mean-square velocity of each hyperbolic event
amp : :obj:tuple or :obj:float
amplitude of each hyperbolic event
wav : :obj:numpy.ndarray
wavelet to be applied to data

Returns
-------
d : :obj:numpy.ndarray
data without wavelet of size :math:[n_x \times n_t]
dwav : :obj:numpy.ndarray
data with wavelet of size :math:[n_x \times n_t]

Notes
-----
Each event is created using the following relation:

.. math::
t_i(x) = \sqrt{t_{0,i}^2 + \frac{x^2}{v_{\text{rms},i}^2}}

"""
if isinstance(t0, (float, int)):
t0 = (t0,)
if isinstance(vrms, (float, int)):
vrms = (vrms,)
if isinstance(amp, (float, int)):
amp = (amp,)

# identify dimensions
dt = t[1] - t[0]
wcenter = int(len(wav) / 2)
nx = np.size(x)
nt = np.size(t) + wcenter
nevents = np.size(t0)

# create events
d = np.zeros((nx, nt))
for ievent in range(nevents):
tevent = np.sqrt(t0[ievent] ** 2 + x ** 2 / vrms[ievent] ** 2)
tevent = (tevent - t[0]) / dt
itevent = tevent.astype(int)
dtevent = tevent - itevent
for ix in range(nx):
if itevent[ix] < nt - 1 and itevent[ix] >= 0:
d[ix, itevent[ix]] += amp[ievent] * (1 - dtevent[ix])
d[ix, itevent[ix] + 1] += amp[ievent] * dtevent[ix]

# filter events with certain wavelet
d, dwav = _filterdata(d, nt, wav, wcenter)
return d, dwav

[docs]def linear3d(x, y, t, v, t0, theta, phi, amp, wav):
r"""Linear 3D events

Create 3d linear events given propagation velocity, intercept time, angles,
and amplitude of each event.

Parameters
----------
x : :obj:numpy.ndarray
space axis in x direction
y : :obj:numpy.ndarray
space axis in y direction
t : :obj:numpy.ndarray
time axis
v : :obj:float
propagation velocity
t0 : :obj:tuple or :obj:float
intercept time at :math:x=0 of each linear event
theta : :obj:tuple or :obj:float
angle in x direction (in degrees) of each linear event
phi : :obj:tuple or :obj:float
angle in y direction (in degrees) of each linear event
amp : :obj:tuple or :obj:float
amplitude of each linear event
wav : :obj:numpy.ndarray
wavelet to be applied to data

Returns
-------
d : :obj:numpy.ndarray
data without wavelet of size
:math:[n_y \times n_x \times n_t]
dwav : :obj:numpy.ndarray
data with wavelet of size
:math:[n_y \times n_x \times n_t]

Notes
-----
Each event is created using the following relation:

.. math::
t_i(x, y) = t_{0,i} + p_{x,i} x + p_{y,i} y

where :math:p_{x,i}=\frac{1}{v} \sin( \theta_i)\cos( \phi_i)
and :math:p_{x,i}=\frac{1}{v} \sin( \theta_i)\sin( \phi_i).

"""
if isinstance(t0, (float, int)):
t0 = (t0,)
if isinstance(theta, (float, int)):
theta = (theta,)
if isinstance(phi, (float, int)):
phi = (phi,)
if isinstance(amp, (float, int)):
amp = (amp,)

# identify dimensions
dt = t[1] - t[0]
wcenter = int(len(wav) / 2)
nx = np.size(x)
ny = np.size(y)
nt = np.size(t) + wcenter
nevents = np.size(t0)

# create events
d = np.zeros((ny, nx, nt))
for ievent in range(nevents):
px = np.sin(np.deg2rad(theta[ievent])) * np.cos(np.deg2rad(phi[ievent])) / v
py = np.sin(np.deg2rad(theta[ievent])) * np.sin(np.deg2rad(phi[ievent])) / v
for iy in range(ny):
tevent = t0[ievent] + px * x + py * y[iy]
tevent = (tevent - t[0]) / dt
itevent = tevent.astype(int)
dtevent = tevent - itevent
for ix in range(nx):
if itevent[ix] < nt - 1 and itevent[ix] >= 0:
d[iy, ix, itevent[ix]] += amp[ievent] * (1 - dtevent[ix])
d[iy, ix, itevent[ix] + 1] += amp[ievent] * dtevent[ix]

# filter events with certain wavelet
d, dwav = _filterdata(d, nt, wav, wcenter)
return d, dwav

[docs]def hyperbolic3d(x, y, t, t0, vrms_x, vrms_y, amp, wav):
r"""Hyperbolic 3D events

Create 3d hyperbolic events given intercept time, root-mean-square
velocities, and amplitude of each event

Parameters
----------
x : :obj:numpy.ndarray
space axis in x direction
y : :obj:numpy.ndarray
space axis in y direction
t : :obj:numpy.ndarray
time axis
t0 : :obj:tuple or :obj:float
intercept time at :math:x=0 of each of hyperbolic event
vrms_x : :obj:tuple or :obj:float
root-mean-square velocity in x direction for each hyperbolic event
vrms_y : :obj:tuple or :obj:float
root-mean-square velocity in y direction for each hyperbolic event
amp : :obj:tuple or :obj:float
amplitude of each hyperbolic event
wav : :obj:numpy.ndarray
wavelet to be applied to data

Returns
-------
d : :obj:numpy.ndarray
data without wavelet of size :math:[n_y \times n_x \times n_t]
dwav : :obj:numpy.ndarray
data with wavelet of size :math:[n_y \times n_x \times n_t]

Notes
-----
Each event is created using the following relation:

.. math::
t_i(x, y) = \sqrt{t_{0,i}^2 + \frac{x^2}{v_{\text{rms}_x, i}^2} +
\frac{y^2}{v_{\text{rms}_y, i}^2}}

Note that velocities do not have a physical meaning here (compared to the
corresponding :func:pylops.utils.seismicevents.hyperbolic2d), they rather
simply control the curvature of the hyperboloid along the spatial axes.

"""
if isinstance(t0, (float, int)):
t0 = (t0,)
if isinstance(vrms_x, (float, int)):
vrms_x = (vrms_x,)
if isinstance(vrms_y, (float, int)):
vrms_y = (vrms_y,)
if isinstance(amp, (float, int)):
amp = (amp,)

# identify dimensions
dt = t[1] - t[0]
wcenter = int(len(wav) / 2)
nx = np.size(x)
ny = np.size(y)
nt = np.size(t) + wcenter
nevents = np.size(t0)

# create events
d = np.zeros((ny, nx, nt))
for ievent in range(nevents):
for iy in range(ny):
tevent = np.sqrt(
t0[ievent] ** 2
+ x ** 2 / vrms_x[ievent] ** 2
+ y[iy] ** 2 / vrms_y[ievent] ** 2
)
tevent = (tevent - t[0]) / dt
itevent = tevent.astype(int)
dtevent = tevent - itevent
for ix in range(nx):
if itevent[ix] < nt - 1 and itevent[ix] >= 0:
d[iy, ix, itevent[ix]] += amp[ievent] * (1 - dtevent[ix])
d[iy, ix, itevent[ix] + 1] += amp[ievent] * dtevent[ix]

# filter events with certain wavelet
d, dwav = _filterdata(d, nt, wav, wcenter)
return d, dwav