In [2]:

    

%matplotlib ipympl

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


    

In [424]:

    

st0 = obspy.read('/home/malcolmw/data/wfs/070/ALCY.2013.070.mseed')
st0.filter('bandpass', freqmin=1, freqmax=20)


    

    
Out[424]:





3 Trace(s) in Stream:
YN.ALCY..EHE | 2013-03-11T00:00:00.004500Z - 2013-03-11T23:59:59.984500Z | 50.0 Hz, 4320000 samples
YN.ALCY..EHN | 2013-03-11T00:00:00.004500Z - 2013-03-11T23:59:59.984500Z | 50.0 Hz, 4320000 samples
YN.ALCY..EHZ | 2013-03-11T00:00:00.004500Z - 2013-03-11T23:59:59.984500Z | 50.0 Hz, 4320000 samples

In [429]:

    

st = st0.slice(
    starttime=obspy.UTCDateTime('2013-03-11T18:25:30'),
    endtime=obspy.UTCDateTime('2013-03-11T18:26:00')
)
# st.normalize()
st.plot(handle=True);


    

In [430]:

    

sample_rate = st[0].stats.sampling_rate

trz = st.select(channel='??Z*')[0].data
trn = st.select(channel='??[N1]*')[0].data
tre = st.select(channel='??[E2]*')[0].data

window_length_in_seconds = 1
nsamp = int(window_length_in_seconds * sample_rate)
shape = (trz.size - nsamp + 1)
xx = np.lib.stride_tricks.as_strided(trz, shape=(shape, nsamp), strides=(trz.itemsize, trz.itemsize))
yy = np.lib.stride_tricks.as_strided(trn, shape=(shape, nsamp), strides=(trz.itemsize, trz.itemsize))
zz = np.lib.stride_tricks.as_strided(tre, shape=(shape, nsamp), strides=(trz.itemsize, trz.itemsize))
xyz = np.stack([xx, yy, zz], axis=1)
eigh = [np.linalg.eigh(np.cov(window)) for window in xyz]
eigval = np.stack([np.flip(eigh[idx][0]) for idx in range(len(eigh))])
eigvec = np.stack([np.fliplr(eigh[idx][1]) for idx in range(len(eigh))])
rect = 1 - (eigval[:, 1] + eigval[:, 2]) / (2 * eigval[:, 0])
plan = 1 - (2 * eigval[:, 2]) / (eigval[:, 0] + eigval[:, 1])
azim = np.arctan2((eigvec[:, 0, 1]), (eigvec[:, 0, 2]))
inc = np.arccos(np.abs(eigvec[:, 0, 0]))


    

In [441]:

    

plt.close('all')

axes = []
fig = plt.figure(figsize=(5, 5))
s_filt = rect * (1 - np.cos(inc))
p_filt = rect * np.cos(inc)
filt = np.diff(p_filt)
ax = fig.add_subplot(4, 1, 1)
axes.append(ax)
# ax.plot(filt, 'k')
ax.plot(rect, label='Rectilinearity')
ax.plot(plan, label='Planarity')
ax.plot(np.cos(inc), label='cos(incidence)')
ax.set_ylim(0, 1)
ax.legend(loc=7)

row = 1
for tr in (trz, trn, tre):
    row += 1
    ax = fig.add_subplot(4, 1, row)
    axes.append(ax)
    tr = tr[len(tr)-len(filt):]
    tr = tr / np.max(np.abs(tr))
    ax.plot(tr, 'r')
    ax.plot(tr * filt, 'k')

for ax in axes:
    ax.set_xlim(0, len(filt))
for ax in axes[1:]:
    ax.set_ylim(-1, 1)
for ax in axes[:-1]:
    ax.set_xticks([])
axes[-1].set_xticklabels(axes[-1].get_xticks() / sample_rate)

plt.subplots_adjust(hspace=0)


    

In [448]:

    

plt.close('all')

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
filt = rect * np.cos(inc)

tr = trz[len(trz)-len(filt):]
ax.plot(tr / np.max(np.abs(tr)), 'k')

param = rect * np.cos(inc)
ax.plot(param)
ax.plot(np.diff(param))


    

    






 
 










    
Out[448]:





[<matplotlib.lines.Line2D at 0x7f15b0d91780>]

In [436]:

    

np.diff?


    

    






Signature: np.diff(a, n=1, axis=-1)
Docstring:
Calculate the n-th discrete difference along the given axis.

The first difference is given by ``out[n] = a[n+1] - a[n]`` along
the given axis, higher differences are calculated by using `diff`
recursively.

Parameters
----------
a : array_like
    Input array
n : int, optional
    The number of times values are differenced. If zero, the input
    is returned as-is.
axis : int, optional
    The axis along which the difference is taken, default is the
    last axis.

Returns
-------
diff : ndarray
    The n-th differences. The shape of the output is the same as `a`
    except along `axis` where the dimension is smaller by `n`. The
    type of the output is the same as the type of the difference
    between any two elements of `a`. This is the same as the type of
    `a` in most cases. A notable exception is `datetime64`, which
    results in a `timedelta64` output array.

See Also
--------
gradient, ediff1d, cumsum

Notes
-----
Type is preserved for boolean arrays, so the result will contain
`False` when consecutive elements are the same and `True` when they
differ.

For unsigned integer arrays, the results will also be unsigned. This
should not be surprising, as the result is consistent with
calculating the difference directly:

>>> u8_arr = np.array([1, 0], dtype=np.uint8)
>>> np.diff(u8_arr)
array([255], dtype=uint8)
>>> u8_arr[1,...] - u8_arr[0,...]
array(255, np.uint8)

If this is not desirable, then the array should be cast to a larger
integer type first:

>>> i16_arr = u8_arr.astype(np.int16)
>>> np.diff(i16_arr)
array([-1], dtype=int16)

Examples
--------
>>> x = np.array([1, 2, 4, 7, 0])
>>> np.diff(x)
array([ 1,  2,  3, -7])
>>> np.diff(x, n=2)
array([  1,   1, -10])

>>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
>>> np.diff(x)
array([[2, 3, 4],
       [5, 1, 2]])
>>> np.diff(x, axis=0)
array([[-1,  2,  0, -2]])

>>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
>>> np.diff(x)
array([1, 1], dtype='timedelta64[D]')
File:      ~/local/anaconda3/envs/py37/lib/python3.7/site-packages/numpy/lib/function_base.py
Type:      function

In [421]:

    

eigvec.shape


    

    
Out[421]:





(702, 3, 3)

In [ ]: