# Copyright 2019 Pascal Audet & Helen Janiszewski
#
# This file is part of OBStools.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
import sys
from scipy.signal import stft, detrend
from scipy.linalg import norm
import matplotlib.pyplot as plt
import numpy as np
import pickle
from obspy.core import Stream, Trace, read
from obstools.atacr import utils, plotting
from pkg_resources import resource_filename
from pathlib import Path
np.seterr(all='ignore')
# np.set_printoptions(threshold=sys.maxsize)
[docs]
class Power(object):
"""
Container for power spectra for each component, with any shape
Attributes
----------
c11 : :class:`~numpy.ndarray`
Power spectral density for component 1 (any shape)
c22 : :class:`~numpy.ndarray`
Power spectral density for component 2 (any shape)
cZZ : :class:`~numpy.ndarray`
Power spectral density for component Z (any shape)
cPP : :class:`~numpy.ndarray`
Power spectral density for component P (any shape)
"""
def __init__(self, c11=None, c22=None, cZZ=None, cPP=None):
self.c11 = c11
self.c22 = c22
self.cZZ = cZZ
self.cPP = cPP
[docs]
class Cross(object):
"""
Container for cross-power spectra for each component pairs, with any shape
Attributes
----------
c12 : :class:`~numpy.ndarray`
Cross-power spectral density for components 1 and 2 (any shape)
c1Z : :class:`~numpy.ndarray`
Cross-power spectral density for components 1 and Z (any shape)
c1P : :class:`~numpy.ndarray`
Cross-power spectral density for components 1 and P (any shape)
c2Z : :class:`~numpy.ndarray`
Cross-power spectral density for components 2 and Z (any shape)
c2P : :class:`~numpy.ndarray`
Cross-power spectral density for components 2 and P (any shape)
cZP : :class:`~numpy.ndarray`
Cross-power spectral density for components Z and P (any shape)
"""
def __init__(self, c12=None, c1Z=None, c1P=None, c2Z=None, c2P=None,
cZP=None):
self.c12 = c12
self.c1Z = c1Z
self.c1P = c1P
self.c2Z = c2Z
self.c2P = c2P
self.cZP = cZP
[docs]
class Rotation(object):
"""
Container for rotated spectra, with any shape
Attributes
----------
cHH : :class:`~numpy.ndarray`
Power spectral density for rotated horizontal component H (any shape)
cHZ : :class:`~numpy.ndarray`
Cross-power spectral density for components H and Z (any shape)
cHP : :class:`~numpy.ndarray`
Cross-power spectral density for components H and P (any shape)
coh : :class:`~numpy.ndarray`
Coherence between horizontal-vertical components
ph : :class:`~numpy.ndarray`
Phase of cross-power spectrum between horizontal-vertical components
ad : :class:`~numpy.ndarray`
Admittance between horizontal-vertical components
tilt_dir : float
Azimuth of the tilt axis
tilt_ang : float
Angle of the tilt axis
coh_value : float
Maximum coherence
phase_value : float
Phase at maximum coherence
admit_value : float
Admittance value at maximum coherence
phi : :class:`~numpy.ndarray`
Directions for which the coherence is calculated
"""
def __init__(self, cHH=None, cHZ=None, cHP=None, coh=None, ph=None,
ad=None, tilt_dir=None, tilt_ang=None, coh_value=None,
phase_value=None, admit_value=None, phi=None):
self.cHH = cHH
self.cHZ = cHZ
self.cHP = cHP
self.coh = coh
self.ph = ph
self.ad = ad
self.tilt_dir = tilt_dir
self.tilt_ang = tilt_ang
self.coh_value = coh_value
self.phase_value = phase_value
self.admit_value = admit_value
self.phi = phi
[docs]
class DayNoise(object):
r"""
A DayNoise object contains attributes that associate
three-component raw (or deconvolved) traces, metadata information
and window parameters. The available methods carry out the quality
control steps and the average daily spectra for windows flagged as
"good".
Note
----
The object is initialized with :class:`~obspy.core.Trace` objects for
H1, H2, HZ and P components. Traces can be empty if data are not
available. Upon saving, those traces are discarded to save disk space.
Attributes
----------
tr1, tr2, trZ, trP : :class:`~obspy.core.Trace` object
Corresponding trace objects for components H1, H2, HZ and HP.
Traces can be empty (i.e., ``Trace()``) for missing components.
window : float
Length of time window in seconds
overlap : float
Fraction of overlap between adjacent windows
key : str
Station key for current object
dt : float
Sampling distance in seconds. Obtained from ``trZ`` object
npts : int
Number of points in time series. Obtained from ``trZ`` object
fs : float
Sampling frequency (in Hz). Obtained from ``trZ`` object
year : str
Year for current object (obtained from UTCDateTime). Obtained from
``trZ`` object
julday : str
Julian day for current object (obtained from UTCDateTime). Obtained
from ``trZ`` object
ncomp : int
Number of available components (either 2, 3 or 4). Obtained from
non-empty ``Trace`` objects
tf_list : Dict
Dictionary of possible transfer functions given the available
components.
Examples
--------
Get demo noise data as DayNoise object
>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
Now check its main attributes
>>> print(*[daynoise.tr1, daynoise.tr2, daynoise.trZ, daynoise.trP], sep="\n")
7D.M08A..BH1 | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
7D.M08A..BH2 | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
7D.M08A..BHZ | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
7D.M08A..BDH | 2012-03-04T00:00:00.000000Z - 2012-03-04T23:59:59.800000Z | 5.0 Hz, 432000 samples
>>> daynoise.window
7200.0
>>> daynoise.overlap
0.3
>>> daynoise.key
'7D.M08A'
>>> daynoise.ncomp
4
>>> daynoise.tf_list
{'ZP': True, 'Z1': True, 'Z2-1': True, 'ZP-21': True, 'ZH': True, 'ZP-H': True}
"""
def __init__(self, tr1=None, tr2=None, trZ=None, trP=None, window=7200.,
overlap=0.3, key=''):
# Load example data if initializing empty object
if tr1 == 'demo' or tr1 == 'Demo':
print("Uploading demo data - March 04, 2012, station 7D.M08A")
exmpl_path = Path(resource_filename('obstools', 'examples'))
fn = exmpl_path / 'data' / '2012.064*.SAC'
st = read(str(fn))
tr1 = st.select(component='1')[0]
tr2 = st.select(component='2')[0]
trZ = st.select(component='Z')[0]
trP = st.select(component='H')[0]
window = 7200.
overlap = 0.3
key = '7D.M08A'
# Check that all traces are valid Trace objects
for tr in [tr1, tr2, trZ, trP]:
if not isinstance(tr, Trace):
raise Exception("Error initializing DayNoise object - "
+ str(tr)+" is not a Trace object")
# Unpack everything
self.tr1 = tr1
self.tr2 = tr2
self.trZ = trZ
self.trP = trP
self.window = window
self.overlap = overlap
self.key = key
# Get trace attributes
zstats = self.trZ.stats
self.dt = zstats.delta
self.npts = zstats.npts
self.fs = zstats.sampling_rate
self.year = zstats.starttime.year
self.julday = zstats.starttime.julday
self.tkey = str(self.year) + '.' + str(self.julday)
# Get number of components for the available, non-empty traces
ncomp = np.sum(
[1 for tr in
Stream(traces=[tr1, tr2, trZ, trP]) if np.any(tr.data)])
self.ncomp = ncomp
# Build list of available transfer functions based on the number of
# components
if self.ncomp == 2:
self.tf_list = {'ZP': True, 'Z1': False, 'Z2-1': False,
'ZP-21': False, 'ZH': False, 'ZP-H': False}
elif self.ncomp == 3:
self.tf_list = {'ZP': False, 'Z1': True, 'Z2-1': True,
'ZP-21': False, 'ZH': True, 'ZP-H': False}
else:
self.tf_list = {'ZP': True, 'Z1': True, 'Z2-1': True,
'ZP-21': True, 'ZH': True, 'ZP-H': True}
self.QC = False
self.av = False
[docs]
def QC_daily_spectra(self, pd=[0.004, 0.2], tol=1.5, alpha=0.05,
smooth=True, fig_QC=False, debug=False, save=None,
form='png'):
"""
Method to determine daily time windows for which the spectra are
anomalous and should be discarded in the calculation of the
transfer functions.
Parameters
----------
pd : list
Frequency corners of passband for calculating the spectra
tol : float
Tolerance threshold. If spectrum > std*tol, window is flagged as
bad
alpha : float
Confidence interval for f-test
smooth : boolean
Determines if the smoothed (True) or raw (False) spectra are used
fig_QC : boolean
Whether or not to produce a figure showing the results of the
quality control
debug : boolean
Whether or not to plot intermediate steps in the QC procedure
for debugging
save : :class:`~pathlib.Path` object
Relative path to figures folder
form : str
File format (e.g., 'png', 'jpg', 'eps')
Attributes
----------
ftX : :class:`~numpy.ndarray`
Windowed Fourier transform for the `X` component (can be either
1, 2, Z or P)
f : :class:`~numpy.ndarray`
Full frequency axis (Hz)
goodwins : list
List of booleans representing whether a window is good (True)
or not (False)
Examples
--------
Perform QC on DayNoise object using default values and plot final
figure
>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra(fig_QC=True)
.. figure:: ../obstools/examples/figures/Figure_3a.png
:align: center
Print out new attribute of DayNoise object
>>> daynoise.goodwins
array([False, True, True, True, True, True, True, True, False,
False, True, True, True, True, True, True], dtype=bool)
"""
# Points in window
ws = int(self.window/self.dt)
# Number of points to overlap
ss = int(self.window*self.overlap/self.dt)
# hanning window
hanning = np.hanning(2*ss)
wind = np.ones(ws)
wind[0:ss] = hanning[0:ss]
wind[-ss:ws] = hanning[ss:ws]
# Get windowed Fourier transforms
ft1 = self.ft1 = None
ft2 = self.ft2 = None
ftZ = self.ftZ = None
ftP = self.ftP = None
# Calculate windowed FFTs and store as transpose
f, t, ftZ = stft(
self.trZ.data, self.fs, return_onesided=False, boundary=None,
padded=False, window=wind, nperseg=ws, noverlap=ss,
detrend='constant')
ftZ = ftZ * ws
self.ftZ = ftZ.T
if self.ncomp == 2 or self.ncomp == 4:
_f, _t, ftP = stft(
self.trP.data, self.fs, return_onesided=False, boundary=None,
padded=False, window=wind, nperseg=ws, noverlap=ss,
detrend='constant')
ftP = ftP * ws
self.ftP = ftP.T
if self.ncomp == 3 or self.ncomp == 4:
_f, _t, ft1 = stft(
self.tr1.data, self.fs, return_onesided=False, boundary=None,
padded=False, window=wind, nperseg=ws, noverlap=ss,
detrend='constant')
_f, _t, ft2 = stft(
self.tr2.data, self.fs, return_onesided=False, boundary=None,
padded=False, window=wind, nperseg=ws, noverlap=ss,
detrend='constant')
ft1 = ft1 * ws
ft2 = ft2 * ws
self.ft1 = ft1.T
self.ft2 = ft2.T
# Store frequency axis
self.f = f
# Get spectrograms for single day-long keys
psd1 = None
psd2 = None
psdZ = None
psdP = None
# Positive frequencies for PSD plots
faxis = int(len(f)/2)
f = f[0:faxis]
psdZ = np.abs(ftZ)**2*2./self.dt
psdZ = psdZ[0:faxis, :]
if self.ncomp == 2 or self.ncomp == 4:
psdP = np.abs(ftP)**2*2/self.dt
psdP = psdP[0:faxis, :]
if self.ncomp == 3 or self.ncomp == 4:
psd1 = np.abs(ft1)**2*2/self.dt
psd1 = psd1[0:faxis, :]
psd2 = np.abs(ft2)**2*2/self.dt
psd2 = psd2[0:faxis, :]
if fig_QC:
if self.ncomp == 2:
plt.figure(1)
plt.subplot(2, 1, 1)
plt.pcolormesh(t, f, np.log(psdZ), shading='auto')
plt.title('Z', fontdict={'fontsize': 8})
plt.subplot(2, 1, 2)
plt.pcolormesh(t, f, np.log(psdP), shading='auto')
plt.title('P', fontdict={'fontsize': 8})
plt.xlabel('Seconds')
plt.tight_layout()
if save:
fname = self.key + '.' + self.tkey + \
'.specgram_Z.P.' + form
if isinstance(save, Path):
fname = save / fname
plt.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plt.show()
elif self.ncomp == 3:
plt.figure(1)
plt.subplot(3, 1, 1)
plt.pcolormesh(t, f, np.log(psd1), shading='auto')
plt.title('H1', fontdict={'fontsize': 8})
plt.subplot(3, 1, 2)
plt.pcolormesh(t, f, np.log(psd2), shading='auto')
plt.title('H2', fontdict={'fontsize': 8})
plt.subplot(3, 1, 3)
plt.pcolormesh(t, f, np.log(psdZ), shading='auto')
plt.title('Z', fontdict={'fontsize': 8})
plt.xlabel('Seconds')
plt.tight_layout()
if save:
fname = self.key + '.' + self.tkey + \
'.specgram_H1.H2.Z.' + form
if isinstance(save, Path):
fname = save / fname
plt.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plt.show()
else:
plt.figure(1)
plt.subplot(4, 1, 1)
plt.pcolormesh(t, f, np.log(psd1), shading='auto')
plt.title('H1', fontdict={'fontsize': 8})
plt.subplot(4, 1, 2)
plt.pcolormesh(t, f, np.log(psd2), shading='auto')
plt.title('H2', fontdict={'fontsize': 8})
plt.subplot(4, 1, 3)
plt.pcolormesh(t, f, np.log(psdZ), shading='auto')
plt.title('Z', fontdict={'fontsize': 8})
plt.subplot(4, 1, 4)
plt.pcolormesh(t, f, np.log(psdP), shading='auto')
plt.title('P', fontdict={'fontsize': 8})
plt.xlabel('Seconds')
plt.tight_layout()
if save:
fname = self.key + '.' + self.tkey + \
'.specgram_H1.H2.Z.P.' + form
if isinstance(save, Path):
fname = save / fname
plt.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plt.show()
# Select bandpass frequencies
ff = (f > pd[0]) & (f < pd[1])
if np.sum([(psd == 0.).any() for psd in
[psd1, psd2, psdZ, psdP] if psd is not None]) > 0.:
smooth = True
if smooth:
# Smooth out the log of the PSDs
sl_psd1 = None
sl_psd2 = None
sl_psdZ = None
sl_psdP = None
sl_psdZ = utils.smooth(np.log(psdZ, where=(psdZ > 0.)), 50, axis=0)
if self.ncomp == 2 or self.ncomp == 4:
sl_psdP = utils.smooth(
np.log(psdP, where=(psdP > 0.)), 50, axis=0)
if self.ncomp == 3 or self.ncomp == 4:
sl_psd1 = utils.smooth(
np.log(psd1, where=(psd1 > 0.)), 50, axis=0)
sl_psd2 = utils.smooth(
np.log(psd2, where=(psd2 > 0.)), 50, axis=0)
else:
# Take the log of the PSDs
sl_psd1 = None
sl_psd2 = None
sl_psdZ = None
sl_psdP = None
sl_psdZ = np.log(psdZ)
if self.ncomp == 2 or self.ncomp == 4:
sl_psdP = np.log(psdP)
if self.ncomp == 3 or self.ncomp == 4:
sl_psd1 = np.log(psd1)
sl_psd2 = np.log(psd2)
# Remove mean of the log PSDs
dsl_psdZ = sl_psdZ[ff, :] - np.mean(sl_psdZ[ff, :], axis=0)
if self.ncomp == 2:
dsl_psdP = sl_psdP[ff, :] - np.mean(sl_psdP[ff, :], axis=0)
dsls = [dsl_psdZ, dsl_psdP]
elif self.ncomp == 3:
dsl_psd1 = sl_psd1[ff, :] - np.mean(sl_psd1[ff, :], axis=0)
dsl_psd2 = sl_psd2[ff, :] - np.mean(sl_psd2[ff, :], axis=0)
dsls = [dsl_psd1, dsl_psd2, dsl_psdZ]
else:
dsl_psd1 = sl_psd1[ff, :] - np.mean(sl_psd1[ff, :], axis=0)
dsl_psd2 = sl_psd2[ff, :] - np.mean(sl_psd2[ff, :], axis=0)
dsl_psdP = sl_psdP[ff, :] - np.mean(sl_psdP[ff, :], axis=0)
dsls = [dsl_psd1, dsl_psd2, dsl_psdZ, dsl_psdP]
if debug:
if self.ncomp == 2:
plt.figure(2)
plt.subplot(2, 1, 1)
plt.semilogx(f, sl_psdZ, 'g', lw=0.5)
plt.subplot(2, 1, 2)
plt.semilogx(f, sl_psdP, 'k', lw=0.5)
plt.tight_layout()
elif self.ncomp == 3:
plt.figure(2)
plt.subplot(3, 1, 1)
plt.semilogx(f, sl_psd1, 'r', lw=0.5)
plt.subplot(3, 1, 2)
plt.semilogx(f, sl_psd2, 'b', lw=0.5)
plt.subplot(3, 1, 3)
plt.semilogx(f, sl_psdZ, 'g', lw=0.5)
plt.tight_layout()
else:
plt.figure(2)
plt.subplot(4, 1, 1)
plt.semilogx(f, sl_psd1, 'r', lw=0.5)
plt.subplot(4, 1, 2)
plt.semilogx(f, sl_psd2, 'b', lw=0.5)
plt.subplot(4, 1, 3)
plt.semilogx(f, sl_psdZ, 'g', lw=0.5)
plt.subplot(4, 1, 4)
plt.semilogx(f, sl_psdP, 'k', lw=0.5)
plt.tight_layout()
plt.show()
# Cycle through to kill high-std-norm windows
moveon = False
goodwins = np.repeat([True], len(t))
indwin = np.argwhere(goodwins == True)
while moveon == False:
ubernorm = np.empty((self.ncomp, np.sum(goodwins)))
for ind_u, dsl in enumerate(dsls):
normvar = np.zeros(np.sum(goodwins))
for ii, tmp in enumerate(indwin):
ind = np.copy(indwin)
ind = np.delete(ind, ii)
normvar[ii] = norm(np.std(dsl[:, ind], axis=1), ord=2)
ubernorm[ind_u, :] = np.median(normvar) - normvar
penalty = np.sum(ubernorm, axis=0)
plt.figure(4)
for i in range(self.ncomp):
plt.plot(range(0, np.sum(goodwins)), detrend(
ubernorm, type='constant')[i], 'o-')
if debug:
plt.show()
else:
plt.close('all')
plt.figure(5)
plt.plot(range(0, np.sum(goodwins)),
np.sum(ubernorm, axis=0), 'o-')
if debug:
plt.show()
else:
plt.close('all')
kill = penalty > tol*np.std(penalty)
if np.sum(kill) == 0:
self.goodwins = goodwins
moveon = True
trypenalty = penalty[np.argwhere(kill == False)].T[0]
if utils.ftest(penalty, 1, trypenalty, 1) < alpha:
goodwins[indwin[kill == True]] = False
indwin = np.argwhere(goodwins == True)
moveon = False
else:
moveon = True
self.goodwins = goodwins
if fig_QC:
power = Power(sl_psd1, sl_psd2, sl_psdZ, sl_psdP)
plot = plotting.fig_QC(
f, power, goodwins, self.ncomp, key=self.key)
# Save or show figure
if save:
fname = self.key + '.' + self.tkey + '.' + 'QC.' + form
if isinstance(save, Path):
fname = save / fname
plot.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plot.show()
self.QC = True
[docs]
def average_daily_spectra(self, calc_rotation=True,
tiltfreqs=[0.005, 0.035],
fig_average=False,
fig_tilt=False,
save=None,
form='png'):
"""
Method to average the daily spectra for good windows. By default, the
method will attempt to calculate the tilt orientation from the
maximum coherence between component H1 and the vertical component HZ
and use the rotated horizontals in the transfer function calculations
for tilt noise removal.
Parameters
----------
calc_rotation : boolean
Whether or not to calculate the tilt direction
tiltfreqs : list
Two floats representing the frequency band at which the tilt is
calculated
fig_average : boolean
Whether or not to produce a figure showing the average daily
spectra
fig_tilt : boolean
Whether or not to produce a figure showing the maximum coherence
and phase between H1 and HZ as function of azimuth measured CW
from H1, and the spectra for these components
save : :class:`~pathlib.Path` object
Relative path to figures folder
form : str
File format (e.g., 'png', 'jpg', 'eps')
Attributes
----------
power : :class:`~obstools.atacr.classes.Power`
Container for the Power spectra
cross : :class:`~obstools.atacr.classes.Cross`
Container for the Cross power spectra
rotation : :class:`~obstools.atacr.classes.Rotation`, optional
Container for the Rotated power and cross spectra
Examples
--------
Average spectra for good windows using default values and plot final
figure
>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra(fig_average=True)
.. figure:: ../obstools/examples/figures/Figure_3b.png
:align: center
Print out available attributes of DayNoise object
>>> daynoise.__dict__.keys()
dict_keys(['tr1', 'tr2', 'trZ', 'trP', 'window', 'overlap', 'key',
'dt', 'npts', 'fs', 'year', 'julday', 'tkey', 'ncomp', 'tf_list',
'QC', 'av', 'f', 'goodwins', 'power', 'cross', 'rotation'])
"""
if not self.QC:
print("Warning: processing daynoise object for " +
"QC_daily_spectra using default values")
self.QC_daily_spectra()
# Extract good windows
c11 = None
c22 = None
cZZ = None
cPP = None
cZZ = np.abs(
np.mean(
self.ftZ[self.goodwins, :]*np.conj(self.ftZ[self.goodwins, :]),
axis=0))
if self.ncomp == 2 or self.ncomp == 4:
cPP = np.abs(
np.mean(
self.ftP[self.goodwins, :]*np.conj(
self.ftP[self.goodwins, :]), axis=0))
if self.ncomp == 3 or self.ncomp == 4:
c11 = np.abs(
np.mean(
self.ft1[self.goodwins, :]*np.conj(
self.ft1[self.goodwins, :]), axis=0))
c22 = np.abs(
np.mean(
self.ft2[self.goodwins, :]*np.conj(
self.ft2[self.goodwins, :]), axis=0))
# Extract bad windows
bc11 = None
bc22 = None
bcZZ = None
bcPP = None
if np.sum(~self.goodwins) > 0:
bcZZ = np.abs(
np.mean(
self.ftZ[~self.goodwins, :]*np.conj(
self.ftZ[~self.goodwins, :]), axis=0))
if self.ncomp == 2 or self.ncomp == 4:
bcPP = np.abs(
np.mean(
self.ftP[~self.goodwins, :]*np.conj(
self.ftP[~self.goodwins, :]), axis=0))
if self.ncomp == 3 or self.ncomp == 4:
bc11 = np.abs(
np.mean(
self.ft1[~self.goodwins, :]*np.conj(
self.ft1[~self.goodwins, :]), axis=0))
bc22 = np.abs(
np.mean(
self.ft2[~self.goodwins, :]*np.conj(
self.ft2[~self.goodwins, :]), axis=0))
# Calculate mean of all good windows if component combinations exist
c12 = None
c1Z = None
c2Z = None
c1P = None
c2P = None
cZP = None
if self.ncomp == 3 or self.ncomp == 4:
c12 = np.mean(
self.ft1[self.goodwins, :] *
np.conj(self.ft2[self.goodwins, :]), axis=0)
c1Z = np.mean(
self.ft1[self.goodwins, :] *
np.conj(self.ftZ[self.goodwins, :]), axis=0)
c2Z = np.mean(
self.ft2[self.goodwins, :] *
np.conj(self.ftZ[self.goodwins, :]), axis=0)
if self.ncomp == 4:
c1P = np.mean(
self.ft1[self.goodwins, :] *
np.conj(self.ftP[self.goodwins, :]), axis=0)
c2P = np.mean(
self.ft2[self.goodwins, :] *
np.conj(self.ftP[self.goodwins, :]), axis=0)
if self.ncomp == 2 or self.ncomp == 4:
cZP = np.mean(
self.ftZ[self.goodwins, :] *
np.conj(self.ftP[self.goodwins, :]), axis=0)
# Store as attributes
self.power = Power(c11, c22, cZZ, cPP)
self.cross = Cross(c12, c1Z, c1P, c2Z, c2P, cZP)
bad = Power(bc11, bc22, bcZZ, bcPP)
if fig_average:
plot = plotting.fig_average(self.f, self.power, bad, self.goodwins,
self.ncomp, key=self.key)
if save:
fname = self.key + '.' + self.tkey + '.' + 'average.' + form
if isinstance(save, Path):
fname = save / fname
plot.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plot.show()
if calc_rotation and self.ncomp >= 3:
cHH, cHZ, cHP, coh, ph, ad, phi, tilt_dir, tilt_ang, coh_value, phase_value, admit_value = \
utils.calculate_tilt(
self.ft1, self.ft2, self.ftZ, self.ftP, self.f,
self.goodwins, tiltfreqs)
self.rotation = Rotation(
cHH, cHZ, cHP, coh, ph, ad, tilt_dir, tilt_ang, coh_value, phase_value, admit_value, phi)
if fig_tilt:
Co = utils.coherence(cHZ, cHH, cZZ)
Ad = utils.admittance(cHZ, cHH)
Ph = utils.phase(cHZ)
plot = plotting.fig_tilt_day(
coh, ph, ad, phi,
Co, Ph, Ad, self.f,
tiltfreqs, tilt_dir, tilt_ang)
# Save or show figure
if save:
fname = self.key + '.' + self.tkey + '.' + 'tilt.' + form
if isinstance(save, Path):
fname = save / fname
plot.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plot.show()
else:
self.rotation = Rotation()
self.av = True
[docs]
def save(self, filename):
"""
Method to save the object to file using `~Pickle`.
Parameters
----------
filename : str
File name
Examples
--------
Run demo through all methods
>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
Save object
>>> daynoise.save('daynoise_demo.pkl')
Check that it has been saved
>>> import glob
>>> glob.glob("./daynoise_demo.pkl")
['./daynoise_demo.pkl']
"""
if not self.av:
print("Warning: saving before having calculated the average " +
"spectra")
# Remove original traces to save disk space
del self.tr1
del self.tr2
del self.trZ
del self.trP
file = open(str(filename), 'wb')
pickle.dump(self, file)
file.close()
[docs]
class StaNoise(object):
"""
A StaNoise object contains attributes that associate
three-component raw (or deconvolved) traces, metadata information
and window parameters.
Note
----
The object is initially a container for
:class:`~obstools.atacr.classes.DayNoise` objects. Once the StaNoise
object is initialized (using the method `init()` or by calling the
`QC_sta_spectra()` method), each individual spectral quantity is unpacked
as an object attribute and the original `DayNoise` objects are removed
from memory. **DayNoise objects cannot be added or appended to the
StaNoise object once this is done**.
In addition, all spectral quantities associated with the
original `DayNoise` objects (now stored as attributes) are discarded as
the object is saved to disk and new container objects are defined and
saved.
Attributes
----------
daylist : list
A list of :class:`~obstools.atacr.classes.DayNoise` objects to process
and produce a station average
initialized : bool
Whether or not the object has been initialized - `False` unless one
of the methods have been called. When `True`, the `daylist` attribute
is deleted from memory
Examples
--------
Initialize empty object
>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise()
Initialize with DayNoise object
>>> from obstools.atacr import DayNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> stanoise = StaNoise(daylist=[daynoise])
Add or append DayNoise object to StaNoise
>>> stanoise = StaNoise()
>>> stanoise += daynoise
>>> stanoise = StaNoise()
>>> stanoise.append(daynoise)
Import demo noise data with 4 DayNoise objects
>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.daylist
[<obstools.atacr.classes.DayNoise at 0x11e3ce8d0>,
<obstools.atacr.classes.DayNoise at 0x121c7ae10>,
<obstools.atacr.classes.DayNoise at 0x121ca5940>,
<obstools.atacr.classes.DayNoise at 0x121e7dd30>]
>>> stanoise.initialized
False
"""
def __init__(self, daylist=None):
def _load_dn(day):
exmpl_path = Path(resource_filename('obstools', 'examples'))
fn = '2012.'+day+'*.SAC'
fn = exmpl_path / 'data' / fn
st = read(str(fn))
tr1 = st.select(component='1')[0]
tr2 = st.select(component='2')[0]
trZ = st.select(component='Z')[0]
trP = st.select(component='H')[0]
window = 7200.
overlap = 0.3
key = '7D.M08A'
return DayNoise(tr1, tr2, trZ, trP, window, overlap, key)
self.daylist = []
self.initialized = False
self.QC = False
self.av = False
if isinstance(daylist, DayNoise):
daylist = [daylist]
elif daylist == 'demo' or daylist == 'Demo':
print("Uploading demo data - March 01 to 04, 2012, station " +
"7D.M08A")
self.daylist = [_load_dn('061'), _load_dn(
'062'), _load_dn('063'), _load_dn('064')]
if not daylist == 'demo' and daylist:
self.daylist.extend(daylist)
def __add__(self, other):
if isinstance(other, DayNoise):
other = StaNoise([other])
if not isinstance(other, StaNoise):
raise TypeError
daylist = self.daylist + other.daylist
return self.__class__(daylist=daylist)
def __iter__(self):
return List(self.daylist).__iter__()
def append(self, daynoise):
if isinstance(daynoise, DayNoise):
self.daylist.append(daynoise)
else:
msg = 'Append only supports a single DayNoise object as argument'
raise TypeError(msg)
return self
def extend(self, daynoise_list):
if isinstance(daynoise_list, list):
for _i in daynoise_list:
# Make sure each item in the list is a Grid object.
if not isinstance(_i, DayNoise):
msg = 'Extend only accepts a list of Daynoise objects.'
raise TypeError(msg)
self.daylist.extend(daynoise_list)
elif isinstance(daynoise_list, StaNoise):
self.daylist.extend(daynoise_list.daylist)
else:
msg = 'Extend only supports a list of DayNoise objects as ' +\
'argument.'
raise TypeError(msg)
return self
[docs]
def init(self):
"""
Method to initialize the `StaNoise` object. This method is used to
unpack the spectral quantities from the original
:class:`~obstools.atacr.classes.DayNoise` objects and allow the
methods to proceed. The original
:class:`~obstools.atacr.classes.DayNoise` objects are deleted from
memory during this process.
Note
----
If the original :class:`~obstools.atacr.classes.DayNoise` objects
have not been processed using their QC and averaging methods, these
will be called first before unpacking into the object attributes.
Attributes
----------
f : :class:`~numpy.ndarray`
Frequency axis for corresponding time sampling parameters
nwins : int
Number of good windows from the
:class:`~obstools.atacr.classes.DayNoise` object
key : str
Station key for current object
ncomp : int
Number of available components (either 2, 3 or 4)
tf_list : Dict
Dictionary of possible transfer functions given the available
components.
c11 : `numpy.ndarray`
Power spectra for component `H1`. Other identical attributes
are available for
the power, cross and rotated spectra: [11, 12, 1Z, 1P, 22, 2Z,
2P, ZZ, ZP, PP, HH, HZ, HP]
phi : `numpy.ndarray`
Array of azimuths used in determining the tilt direction
tilt_dir : float
Tilt direction from maximum coherence between rotated `H1` and
`HZ` components
tilt_ang : float
Tilt angle from maximum coherence between rotated `H1` and
`HZ` components
QC : bool
Whether or not the method
:func:`~obstools.atacr.classes.StaNoise.QC_sta_spectra` has
been called.
av : bool
Whether or not the method
:func:`~obstools.atacr.classes.StaNoise.average_sta_spectra` has
been called.
Examples
--------
Initialize demo data
>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.init()
Check that `daylist` attribute has been deleted
>>> stanoise.daylist
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
AttributeError: 'StaNoise' object has no attribute 'daylist'
>>> stanoise.__dict__.keys()
dict_keys(['initialized', 'c11', 'c22', 'cZZ', 'cPP', 'c12', 'c1Z',
'c1P', 'c2Z', 'c2P', 'cZP', 'cHH', 'cHZ', 'cHP', 'phi', 'tilt_dir',
'tilt_ang', 'f', 'nwins', 'ncomp', 'key', 'tf_list', 'QC', 'av'])
"""
# First, check that the StaNoise object contains at least two
# DayNoise objects
if len(self.daylist) < 2:
raise Exception(
"StaNoise requires at least two DayNoise objects to execute " +
"its methods")
for dn in self.daylist:
if not dn.QC:
dn.QC_daily_spectra()
if not dn.av:
dn.average_daily_spectra()
# Then unpack the DayNoise objects
self.c11 = np.array([dn.power.c11 for dn in self.daylist]).T
self.c22 = np.array([dn.power.c22 for dn in self.daylist]).T
self.cZZ = np.array([dn.power.cZZ for dn in self.daylist]).T
self.cPP = np.array([dn.power.cPP for dn in self.daylist]).T
self.c12 = np.array([dn.cross.c12 for dn in self.daylist]).T
self.c1Z = np.array([dn.cross.c1Z for dn in self.daylist]).T
self.c1P = np.array([dn.cross.c1P for dn in self.daylist]).T
self.c2Z = np.array([dn.cross.c2Z for dn in self.daylist]).T
self.c2P = np.array([dn.cross.c2P for dn in self.daylist]).T
self.cZP = np.array([dn.cross.cZP for dn in self.daylist]).T
self.cHH = np.array([dn.rotation.cHH for dn in self.daylist]).T
self.cHZ = np.array([dn.rotation.cHZ for dn in self.daylist]).T
self.cHP = np.array([dn.rotation.cHP for dn in self.daylist]).T
self.phi = self.daylist[0].rotation.phi
self.tilt_dir = self.daylist[0].rotation.tilt_dir
self.tilt_ang = self.daylist[0].rotation.tilt_ang
self.f = self.daylist[0].f
self.nwins = np.array([np.sum(dn.goodwins) for dn in self.daylist])
self.ncomp = np.min([dn.ncomp for dn in self.daylist])
self.key = self.daylist[0].key
# Build list of available transfer functions for future use
if self.ncomp == 2:
self.tf_list = {'ZP': True, 'Z1': False, 'Z2-1': False,
'ZP-21': False, 'ZH': False, 'ZP-H': False}
elif self.ncomp == 3:
self.tf_list = {'ZP': False, 'Z1': True, 'Z2-1': True,
'ZP-21': False, 'ZH': False, 'ZP-H': False}
else:
self.tf_list = {'ZP': True, 'Z1': True, 'Z2-1': True,
'ZP-21': True, 'ZH': False, 'ZP-H': False}
self.initialized = True
self.QC = False
self.av = False
# Remove DayNoise objects from memory
del self.daylist
[docs]
def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05,
fig_QC=False, debug=False, save=None, form='png'):
"""
Method to determine the days (for given time window) for which the
spectra are anomalous and should be discarded in the calculation of
the long-term transfer functions.
Parameters
----------
pd : list
Frequency corners of passband for calculating the spectra
tol : float
Tolerance threshold. If spectrum > std*tol, window is flagged as
bad
alpha : float
Confidence interval for f-test
fig_QC : boolean
Whether or not to produce a figure showing the results of the
quality control
debug : boolean
Whether or not to plot intermediate steps in the QC procedure for
debugging
save : :class:`~pathlib.Path` object
Relative path to figures folder
form : str
File format (e.g., 'png', 'jpg', 'eps')
Attributes
----------
gooddays : list
List of booleans representing whether a day is good (True) or not
(False)
Examples
--------
Import demo data, call method and generate final figure
>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra(fig_QC=True)
>>> stanoise.QC
True
"""
if self.initialized:
raise Exception("Object has been initialized already - " +
"list of DayNoise objects has been lost and " +
"method cannot proceed")
else:
self.init()
# Select bandpass frequencies
ff = (self.f > pd[0]) & (self.f < pd[1])
# Extract only positive frequencies
faxis = self.f > 0
# Smooth out the log of the PSDs
sl_cZZ = None
sl_c11 = None
sl_c22 = None
sl_cPP = None
sl_cZZ = utils.smooth(np.log(self.cZZ), 50, axis=0)
if self.ncomp == 2 or self.ncomp == 4:
sl_cPP = utils.smooth(np.log(self.cPP), 50, axis=0)
if self.ncomp == 3 or self.ncomp == 4:
sl_c11 = utils.smooth(np.log(self.c11), 50, axis=0)
sl_c22 = utils.smooth(np.log(self.c22), 50, axis=0)
# Remove mean of the log PSDs
dsl_cZZ = sl_cZZ[ff, :] - np.mean(sl_cZZ[ff, :], axis=0)
if self.ncomp == 2:
dsl_cPP = sl_cPP[ff, :] - np.mean(sl_cPP[ff, :], axis=0)
dsls = [dsl_cZZ, dsl_cPP]
elif self.ncomp == 3:
dsl_c11 = sl_c11[ff, :] - np.mean(sl_c11[ff, :], axis=0)
dsl_c22 = sl_c22[ff, :] - np.mean(sl_c22[ff, :], axis=0)
dsls = [dsl_c11, dsl_c22, dsl_cZZ]
else:
dsl_c11 = sl_c11[ff, :] - np.mean(sl_c11[ff, :], axis=0)
dsl_c22 = sl_c22[ff, :] - np.mean(sl_c22[ff, :], axis=0)
dsl_cPP = sl_cPP[ff, :] - np.mean(sl_cPP[ff, :], axis=0)
dsls = [dsl_c11, dsl_c22, dsl_cZZ, dsl_cPP]
if debug:
if self.ncomp == 2:
plt.figure(2)
plt.subplot(2, 1, 1)
plt.semilogx(self.f[faxis], sl_cZZ[faxis], 'g', lw=0.5)
plt.subplot(2, 1, 2)
plt.semilogx(self.f[faxis], sl_cPP[faxis], 'k', lw=0.5)
plt.tight_layout()
elif self.ncomp == 3:
plt.figure(2)
plt.subplot(3, 1, 1)
plt.semilogx(self.f[faxis], sl_c11[faxis], 'r', lw=0.5)
plt.subplot(3, 1, 2)
plt.semilogx(self.f[faxis], sl_c22[faxis], 'b', lw=0.5)
plt.subplot(3, 1, 3)
plt.semilogx(self.f[faxis], sl_cZZ[faxis], 'g', lw=0.5)
plt.tight_layout()
else:
plt.figure(2)
plt.subplot(4, 1, 1)
plt.semilogx(self.f[faxis], sl_c11[faxis], 'r', lw=0.5)
plt.subplot(4, 1, 2)
plt.semilogx(self.f[faxis], sl_c22[faxis], 'b', lw=0.5)
plt.subplot(4, 1, 3)
plt.semilogx(self.f[faxis], sl_cZZ[faxis], 'g', lw=0.5)
plt.subplot(4, 1, 4)
plt.semilogx(self.f[faxis], sl_cPP[faxis], 'k', lw=0.5)
plt.tight_layout()
plt.show()
# Cycle through to kill high-std-norm windows
moveon = False
gooddays = np.repeat([True], self.cZZ.shape[1])
indwin = np.argwhere(gooddays == True)
while moveon == False:
ubernorm = np.empty((self.ncomp, np.sum(gooddays)))
for ind_u, dsl in enumerate(dsls):
normvar = np.zeros(np.sum(gooddays))
for ii, tmp in enumerate(indwin):
ind = np.copy(indwin)
ind = np.delete(ind, ii)
normvar[ii] = norm(np.std(dsl[:, ind], axis=1), ord=2)
ubernorm[ind_u, :] = np.median(normvar) - normvar
penalty = np.sum(ubernorm, axis=0)
if debug:
plt.figure(4)
for i in range(self.ncomp):
plt.plot(range(0, np.sum(gooddays)), detrend(
ubernorm, type='constant')[i], 'o-')
plt.show()
plt.figure(5)
plt.plot(range(0, np.sum(gooddays)),
np.sum(ubernorm, axis=0), 'o-')
plt.show()
kill = penalty > tol*np.std(penalty)
if np.sum(kill) == 0:
self.gooddays = gooddays
self.QC = True
moveon = True
trypenalty = penalty[np.argwhere(kill == False)].T[0]
if utils.ftest(penalty, 1, trypenalty, 1) < alpha:
gooddays[indwin[kill == True]] = False
indwin = np.argwhere(gooddays == True)
moveon = False
else:
moveon = True
self.gooddays = gooddays
self.QC = True
if fig_QC:
power = Power(sl_c11, sl_c22, sl_cZZ, sl_cPP)
plot = plotting.fig_QC(self.f, power, gooddays,
self.ncomp, key=self.key)
if save:
fname = self.key + '.' + 'QC.' + form
if isinstance(save, Path):
fname = save / fname
plot.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plot.show()
[docs]
def average_sta_spectra(self, fig_average=False, save=None, form='png'):
r"""
Method to average the daily station spectra for good windows.
Parameters
----------
fig_average : boolean
Whether or not to produce a figure showing the average daily
spectra
save : :class:`~pathlib.Path` object
Relative path to figures folder
form : str
File format (e.g., 'png', 'jpg', 'eps')
Attributes
----------
power : :class:`~obstools.atacr.classes.Power`
Container for the Power spectra
cross : :class:`~obstools.atacr.classes.Cross`
Container for the Cross power spectra
rotation : :class:`~obstools.atacr.classes.Cross`, optional
Container for the Rotated power and cross spectra
Examples
--------
Average daily spectra for good days using default values and produce
final figure
>>> obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
"""
if not self.QC:
print(
"Warning: processing stanoise object for QC_sta_spectra " +
"using default values")
self.QC_sta_spectra()
# Power spectra
c11 = None
c22 = None
cZZ = None
cPP = None
cZZ = np.sum(self.cZZ[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
if self.ncomp == 2 or self.ncomp == 4:
cPP = np.sum(self.cPP[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
if self.ncomp == 3 or self.ncomp == 4:
c11 = np.sum(self.c11[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
c22 = np.sum(self.c22[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
# Bad days - for plotting
bc11 = None
bc22 = None
bcZZ = None
bcPP = None
if np.sum(~self.gooddays) > 0:
bcZZ = np.sum(
self.cZZ[:, ~self.gooddays]*self.nwins[~self.gooddays],
axis=1)/np.sum(self.nwins[~self.gooddays])
if self.ncomp == 2 or self.ncomp == 4:
bcPP = np.sum(
self.cPP[:, ~self.gooddays]*self.nwins[~self.gooddays],
axis=1)/np.sum(self.nwins[~self.gooddays])
if self.ncomp == 3 or self.ncomp == 4:
bc11 = np.sum(
self.c11[:, ~self.gooddays]*self.nwins[~self.gooddays],
axis=1)/np.sum(self.nwins[~self.gooddays])
bc22 = np.sum(
self.c22[:, ~self.gooddays]*self.nwins[~self.gooddays],
axis=1)/np.sum(self.nwins[~self.gooddays])
# Cross spectra
c12 = None
c1Z = None
c2Z = None
c1P = None
c2P = None
cZP = None
if self.ncomp == 3 or self.ncomp == 4:
c12 = np.sum(self.c12[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
c1Z = np.sum(self.c1Z[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
c2Z = np.sum(self.c2Z[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
if self.ncomp == 4:
c1P = np.sum(self.c1P[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
c2P = np.sum(self.c2P[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
if self.ncomp == 2 or self.ncomp == 4:
cZP = np.sum(self.cZP[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
# Rotated component
cHH = None
cHZ = None
cHP = None
if self.ncomp == 3 or self.ncomp == 4:
cHH = np.sum(self.cHH[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
cHZ = np.sum(self.cHZ[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
if self.ncomp == 4:
cHP = np.sum(self.cHP[:, self.gooddays]*self.nwins[self.gooddays],
axis=1)/np.sum(self.nwins[self.gooddays])
self.power = Power(c11, c22, cZZ, cPP)
self.cross = Cross(c12, c1Z, c1P, c2Z, c2P, cZP)
self.rotation = Rotation(cHH, cHZ, cHP)
bad = Power(bc11, bc22, bcZZ, bcPP)
if fig_average:
plot = plotting.fig_average(
self.f, self.power, bad,
self.gooddays, self.ncomp, key=self.key)
if save:
fname = self.key + '.' + 'average.' + form
if isinstance(save, Path):
fname = save / fname
plot.savefig(
str(fname), dpi=300, bbox_inches='tight', format=form)
else:
plot.show()
self.av = True
[docs]
def save(self, filename):
"""
Method to save the object to file using `~Pickle`.
Parameters
----------
filename : str
File name
Examples
--------
Run demo through all methods
>>> from obstools.atacr import StaNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
Save object
>>> stanoise.save('stanoise_demo.pkl')
Check that it has been saved
>>> import glob
>>> glob.glob("./stanoise_demo.pkl")
['./stanoise_demo.pkl']
"""
if not self.av:
print("Warning: saving before having calculated the average " +
"spectra")
# Remove traces to save disk space
del self.c11
del self.c22
del self.cZZ
del self.cPP
del self.c12
del self.c1Z
del self.c1P
del self.c2Z
del self.c2P
del self.cZP
file = open(filename, 'wb')
pickle.dump(self, file)
file.close()
[docs]
class TFNoise(object):
"""
A TFNoise object contains attributes that store the transfer function
information from multiple components (and component combinations).
Note
----
The object is initialized with either a processed
:class:`~obstools.atacr.classes.DayNoise` or
:class:`~obstools.atacr.classes.StaNoise` object. Each individual
spectral quantity is unpacked as an object attribute, but all of them
are discarded as the object is saved to disk and new container objects
are defined and saved.
Attributes
----------
f : :class:`~numpy.ndarray`
Frequency axis for corresponding time sampling parameters
c11 : `numpy.ndarray`
Power spectra for component `H1`. Other identical attributes are
available for the power, cross and rotated spectra:
[11, 12, 1Z, 1P, 22, 2Z, 2P, ZZ, ZP, PP, HH, HZ, HP]
tilt : float
Tilt direction from maximum coherence between rotated `H1` and
`HZ` components
tf_list : Dict
Dictionary of possible transfer functions given the available
components.
Examples
--------
Initialize a TFNoise object with a DayNoise object. The DayNoise
object must be processed for QC and averaging, otherwise the TFNoise
object will not initialize.
>>> from obstools.atacr import DayNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> tfnoise = TFNoise(daynoise)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/Users/pascalaudet/Softwares/Python/Projects/dev/OBStools/obstools/atacr/classes.py", line 1215, in __init__
Exception: Error: Noise object has not been processed (QC and averaging) - aborting
Now re-initialized with a processed DayNoise object
>>> from obstools.atacr import DayNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise = TFNoise(daynoise)
Initialize a TFNoise object with a processed StaNoise object
>>> from obstools.atacr import StaNoise, TFNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise = TFNoise(stanoise)
"""
def __init__(self, objnoise=None):
if (not objnoise and not isinstance(objnoise, DayNoise) and
not isinstance(objnoise, StaNoise)):
msg = "Error: A TFNoise object must be initialized with only " +\
"one of type DayNoise or StaNoise object"
raise TypeError(msg)
if not objnoise.av:
raise Exception(
"Error: Noise object has not been processed (QC and " +
"averaging) - aborting")
self.f = objnoise.f
self.c11 = objnoise.power.c11
self.c22 = objnoise.power.c22
self.cZZ = objnoise.power.cZZ
self.cPP = objnoise.power.cPP
self.cHH = objnoise.rotation.cHH
self.cHZ = objnoise.rotation.cHZ
self.cHP = objnoise.rotation.cHP
self.c12 = objnoise.cross.c12
self.c1Z = objnoise.cross.c1Z
self.c1P = objnoise.cross.c1P
self.c2Z = objnoise.cross.c2Z
self.c2P = objnoise.cross.c2P
self.cZP = objnoise.cross.cZP
self.tilt = objnoise.rotation.tilt
self.tf_list = objnoise.tf_list
[docs]
class TfDict(dict):
def __init__(self):
self = dict()
def add(self, key, value):
self[key] = value
[docs]
def transfer_func(self):
"""
Method to calculate transfer functions between multiple
components (and component combinations) from the averaged
(daily or station-averaged) noise spectra.
Attributes
----------
transfunc : :class:`~obstools.atacr.classes.TFNoise.TfDict`
Container Dictionary for all possible transfer functions
Examples
--------
Calculate transfer functions for a DayNoise object
>>> from obstools.atacr import DayNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise = TFNoise(daynoise)
>>> tfnoise.transfer_func()
>>> tfnoise.transfunc.keys()
dict_keys(['ZP', 'Z1', 'Z2-1', 'ZP-21', 'ZH', 'ZP-H'])
Calculate transfer functions for a StaNoise object
>>> from obstools.atacr import StaNoise, TFNoise
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise = TFNoise(stanoise)
>>> tfnoise.transfer_func()
>>> tfnoise.transfunc.keys()
dict_keys(['ZP', 'Z1', 'Z2-1', 'ZP-21'])
"""
transfunc = self.TfDict()
for key, value in self.tf_list.items():
if key == 'ZP':
if value:
tf_ZP = {'TF_ZP': self.cZP/self.cPP}
transfunc.add('ZP', tf_ZP)
elif key == 'Z1':
if value:
tf_Z1 = {'TF_Z1': np.conj(self.c1Z)/self.c11}
transfunc.add('Z1', tf_Z1)
elif key == 'Z2-1':
if value:
lc1c2 = np.conj(self.c12)/self.c11
coh_12 = utils.coherence(self.c12, self.c11, self.c22)
gc2c2_c1 = self.c22*(1. - coh_12)
gc2cZ_c1 = np.conj(self.c2Z) - np.conj(lc1c2*self.c1Z)
lc2cZ_c1 = gc2cZ_c1/gc2c2_c1
tf_Z2_1 = {'TF_21': lc1c2, 'TF_Z2-1': lc2cZ_c1}
transfunc.add('Z2-1', tf_Z2_1)
elif key == 'ZP-21':
if value:
lc1cZ = np.conj(self.c1Z)/self.c11
lc1c2 = np.conj(self.c12)/self.c11
lc1cP = np.conj(self.c1P)/self.c11
coh_12 = utils.coherence(self.c12, self.c11, self.c22)
coh_1P = utils.coherence(self.c1P, self.c11, self.cPP)
gc2c2_c1 = self.c22*(1. - coh_12)
gcPcP_c1 = self.cPP*(1. - coh_1P)
gc2cZ_c1 = np.conj(self.c2Z) - np.conj(lc1c2*self.c1Z)
gcPcZ_c1 = self.cZP - np.conj(lc1cP*self.c1Z)
gc2cP_c1 = np.conj(self.c2P) - np.conj(lc1c2*self.c1P)
lc2cP_c1 = gc2cP_c1/gc2c2_c1
lc2cZ_c1 = gc2cZ_c1/gc2c2_c1
coh_c2cP_c1 = utils.coherence(gc2cP_c1, gc2c2_c1,
gcPcP_c1)
gcPcP_c1c2 = gcPcP_c1*(1. - coh_c2cP_c1)
gcPcZ_c1c2 = gcPcZ_c1 - np.conj(lc2cP_c1)*gc2cZ_c1
lcPcZ_c2c1 = gcPcZ_c1c2/gcPcP_c1c2
tf_ZP_21 = {'TF_Z1': lc1cZ, 'TF_21': lc1c2,
'TF_P1': lc1cP, 'TF_P2-1': lc2cP_c1,
'TF_Z2-1': lc2cZ_c1, 'TF_ZP-21': lcPcZ_c2c1}
transfunc.add('ZP-21', tf_ZP_21)
elif key == 'ZH':
if value:
tf_ZH = {'TF_ZH': np.conj(self.cHZ)/self.cHH}
transfunc.add('ZH', tf_ZH)
elif key == 'ZP-H':
if value:
lcHcP = np.conj(self.cHP)/self.cHH
coh_HP = utils.coherence(self.cHP, self.cHH, self.cPP)
gcPcP_cH = self.cPP*(1. - coh_HP)
gcPcZ_cH = self.cZP - np.conj(lcHcP*self.cHZ)
lcPcZ_cH = gcPcZ_cH/gcPcP_cH
tf_ZP_H = {'TF_PH': lcHcP, 'TF_ZP-H': lcPcZ_cH}
transfunc.add('ZP-H', tf_ZP_H)
else:
raise Exception('Incorrect key')
self.transfunc = transfunc
[docs]
def save(self, filename):
"""
Method to save the object to file using `~Pickle`.
Parameters
----------
filename : str
File name
Examples
--------
Run demo through all methods
>>> from obstools.atacr import DayNoise, StaNoise, TFNoise
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise_day = TFNoise(daynoise)
>>> tfnoise_day.transfer_func()
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise_sta = TFNoise(stanoise)
>>> tfnoise_sta.transfer_func()
Save object
>>> tfnoise_day.save('tf_daynoise_demo.pkl')
>>> tfnoise_sta.save('tf_stanoise_demo.pkl')
Check that everything has been saved
>>> import glob
>>> glob.glob("./tf_daynoise_demo.pkl")
['./tf_daynoise_demo.pkl']
>>> glob.glob("./tf_stanoise_demo.pkl")
['./tf_stanoise_demo.pkl']
"""
if not self.transfunc:
print("Warning: saving before having calculated the transfer " +
"functions")
# Remove traces to save disk space
del self.c11
del self.c22
del self.cZZ
del self.cPP
del self.cHH
del self.cHZ
del self.cHP
del self.c12
del self.c1Z
del self.c1P
del self.c2Z
del self.c2P
del self.cZP
file = open(filename, 'wb')
pickle.dump(self, file)
file.close()
[docs]
class EventStream(object):
"""
An EventStream object contains attributes that store station-event
metadata and methods for applying the transfer functions to the various
components and produce corrected/cleaned vertical components.
Note
----
The object is initialized with :class:`~obspy.core.Trace` objects for
H1, H2, HZ and P components. Traces can be empty if data are not
available. Based on the available components, a list of
possible corrections is determined automatically.
Attributes
----------
tr1, tr2, trZ, trP : :class:`~obspy.core.Trace` object
Corresponding trace objects for components H1, H2, HZ and HP.
Traces can be empty (i.e., ``Trace()``) for missing components.
key : str
Station key for current object
evtime : :class:`~obspy.core.UTCDateTime`
Origin time of seismic event.
tstamp : str
Time stamp for event
prefix : str
Associated prefix of SAC files
npts : int
Number of points in time series.
fs : float
Sampling frequency (in Hz).
dt : float
Sampling distance in seconds.
ncomp : int
Number of available components (either 2, 3 or 4). Obtained from
non-empty ``Trace`` objects
ev_list : Dict
Dictionary of possible transfer functions given the available
components. This is determined during initialization.
Examples
--------
Get demo earthquake data as EventStream object
>>> from obstools.atacr import EventStream
>>> evstream = EventStream('demo')
Uploading demo earthquake data - March 09, 2012, station 7D.M08A
>>> evstream.__dict__.keys()
dict_keys(['tr1', 'tr2', 'trZ', 'trP, 'key', 'evtime',
'tstamp', 'prefix', 'npts', fs', 'dt', 'ncomp', 'ev_list'])
Plot the raw traces
>>> import obstools.atacr.plotting as atplot
>>> figure = atplot.fig_event_raw(evstream, fmin=1./150., fmax=2.)
>>> figure.show()
.. figure:: ../obstools/examples/figures/Figure_11.png
:align: center
"""
def __init__(self, tr1=Trace(), tr2=Trace(), trZ=Trace(), trP=Trace(),
correct=False):
if tr1 == 'demo':
print("Uploading demo earthquake data - March 09, 2012, " +
"station 7D.M08A")
exmpl_path = Path(resource_filename('obstools', 'examples'))
fn = exmpl_path / 'event' / '2012.069*.SAC'
st = read(str(fn))
tr1 = st.select(component='1')[0]
tr2 = st.select(component='2')[0]
trZ = st.select(component='Z')[0]
trP = st.select(component='H')[0]
ncomp = np.sum([np.any(tr.data) for tr in [tr1, tr2, trZ, trP]])
if ncomp <= 1 or len(trZ.data) == 0:
raise Exception(
"Incorrect initialization of EventStream object: " +
"missing a vertical component or too few components")
self.tr1 = tr1
self.tr2 = tr2
self.trZ = trZ
self.trP = trP
zstats = trZ.stats
self.key = zstats.network + '.' + zstats.station
self.evtime = zstats.starttime
# Time stamp
tstamp = str(self.evtime.year).zfill(4)+'.' + \
str(self.evtime.julday).zfill(3)+'.'
tstamp = tstamp + str(self.evtime.hour).zfill(2) + \
'.'+str(self.evtime.minute).zfill(2)
self.tstamp = tstamp
self.prefix = self.key + '.' + self.tstamp
self.npts = zstats.npts
self.fs = zstats.sampling_rate
self.dt = zstats.delta
self.ncomp = ncomp
# Build list of available transfer functions for future use
if self.ncomp == 2:
self.ev_list = {'ZP': True, 'Z1': False, 'Z2-1': False,
'ZP-21': False, 'ZH': False, 'ZP-H': False}
elif self.ncomp == 3:
self.ev_list = {'ZP': False, 'Z1': True, 'Z2-1': True,
'ZP-21': False, 'ZH': True, 'ZP-H': False}
else:
self.ev_list = {'ZP': True, 'Z1': True, 'Z2-1': True,
'ZP-21': True, 'ZH': True, 'ZP-H': True}
[docs]
class CorrectDict(dict):
def __init__(self):
self = dict()
def add(self, key, value):
self[key] = value
[docs]
def correct_data(self, tfnoise):
"""
Method to apply transfer functions between multiple components (and
component combinations) to produce corrected/cleaned vertical
components.
Parameters
----------
tfnoise : :class:`~obstools.atacr.classes.TFNoise`
Object that contains the noise transfer functions used in the
correction
Attributes
----------
correct : :class:`~obstools.atacr.classes.EventStream.CorrectDict`
Container Dictionary for all possible corrections from the
transfer functions
Examples
--------
Let's carry through the correction of the vertical component for a
single day of noise, say corresponding to the noise recorded on March
04, 2012. In practice, the DayNoise object should correspond to the
same day at that of the recorded earthquake to avoid bias in the
correction.
>>> from obstools.atacr import DayNoise, TFNoise, EventStream
>>> daynoise = DayNoise('demo')
Uploading demo data - March 04, 2012, station 7D.M08A
>>> daynoise.QC_daily_spectra()
>>> daynoise.average_daily_spectra()
>>> tfnoise_day = TFNoise(daynoise)
>>> tfnoise_day.transfer_func()
>>> evstream = EventStream('demo')
Uploading demo earthquake data - March 09, 2012, station 7D.M08A
>>> evstream.correct_data(tfnoise_day)
Plot the corrected traces
>>> import obstools.atacr.plotting as atplot
>>> figure = atplot.fig_event_corrected(evstream, tfnoise_day.tf_list)
>>> figure.show()
.. figure:: ../obstools/examples/figures/Figure_corrected_march04.png
:align: center
Carry out the same exercise but this time using a StaNoise object
>>> from obstools.atacr import StaNoise, TFNoise, EventStream
>>> stanoise = StaNoise('demo')
Uploading demo data - March 01 to 04, 2012, station 7D.M08A
>>> stanoise.QC_sta_spectra()
>>> stanoise.average_sta_spectra()
>>> tfnoise_sta = TFNoise(stanoise)
>>> tfnoise_sta.transfer_func()
>>> evstream = EventStream('demo')
Uploading demo earthquake data - March 09, 2012, station 7D.M08A
>>> evstream.correct_data(tfnoise_sta)
Plot the corrected traces
>>> import obstools.atacr.plot as plot
>>> plot.fig_event_corrected(evstream, tfnoise_sta.tf_list)
.. figure:: ../obstools/examples/figures/Figure_corrected_sta.png
:align: center
.. warning::
If the noise window and event window are not identical, they cannot
be compared on the same frequency axis and the code will exit. Make
sure you are using identical time samples in both the noise and
event windows.
"""
if not tfnoise.transfunc:
raise Exception("Error: Object TFNoise has no transfunc " +
"attribute - aborting")
correct = self.CorrectDict()
# Extract list and transfer functions available
tf_list = tfnoise.tf_list
transfunc = tfnoise.transfunc
# Extract traces
trZ = self.trZ
tr1 = self.tr1
tr2 = self.tr2
trP = self.trP
# Get Fourier spectra
ft1 = None
ft2 = None
ftZ = None
ftP = None
ft1 = None
ft2 = None
ftZ = None
ftP = None
ftZ = np.fft.fft(trZ, n=self.npts)
if self.ncomp == 2 or self.ncomp == 4:
ftP = np.fft.fft(trP, n=self.npts)
if self.ncomp == 3 or self.ncomp == 4:
ft1 = np.fft.fft(tr1, n=self.npts)
ft2 = np.fft.fft(tr2, n=self.npts)
f = np.fft.fftfreq(self.npts, d=self.dt)
if not np.allclose(f, tfnoise.f):
raise(Exception(
'Frequency axes are different: ', f, tfnoise.f,
' - the noise and event windows are not the same, aborting'))
for key, value in tf_list.items():
if key == 'ZP' and self.ev_list[key]:
if value and tf_list[key]:
TF_ZP = transfunc[key]['TF_ZP']
corrspec = ftZ - TF_ZP*ftP
corrtime = np.real(np.fft.ifft(corrspec))
correct.add('ZP', corrtime)
if key == 'Z1' and self.ev_list[key]:
if value and tf_list[key]:
TF_Z1 = transfunc[key]['TF_Z1']
corrspec = ftZ - TF_Z1*ft1
corrtime = np.real(np.fft.ifft(corrspec))
correct.add('Z1', corrtime)
if key == 'Z2-1' and self.ev_list[key]:
if value and tf_list[key]:
TF_Z1 = transfunc['Z1']['TF_Z1']
TF_21 = transfunc[key]['TF_21']
TF_Z2_1 = transfunc[key]['TF_Z2-1']
corrspec = ftZ - TF_Z1*ft1 - (ft2 - ft1*TF_21)*TF_Z2_1
corrtime = np.real(np.fft.ifft(corrspec))
correct.add('Z2-1', corrtime)
if key == 'ZP-21' and self.ev_list[key]:
if value and tf_list[key]:
TF_Z1 = transfunc[key]['TF_Z1']
TF_21 = transfunc[key]['TF_21']
TF_Z2_1 = transfunc[key]['TF_Z2-1']
TF_P1 = transfunc[key]['TF_P1']
TF_P2_1 = transfunc[key]['TF_P2-1']
TF_ZP_21 = transfunc[key]['TF_ZP-21']
corrspec = ftZ - TF_Z1*ft1 - \
(ft2 - ft1*TF_21)*TF_Z2_1 - \
(ftP - ft1*TF_P1 -
(ft2 - ft1*TF_21)*TF_P2_1)*TF_ZP_21
corrtime = np.real(np.fft.ifft(corrspec))
correct.add('ZP-21', corrtime)
if key == 'ZH' and self.ev_list[key]:
if value and tf_list[key]:
# Rotate horizontals
ftH = utils.rotate_dir(ft1, ft2, tfnoise.tilt)
TF_ZH = transfunc[key]['TF_ZH']
corrspec = ftZ - TF_ZH*ftH
corrtime = np.real(np.fft.ifft(corrspec))
correct.add('ZH', corrtime)
if key == 'ZP-H' and self.ev_list[key]:
if value and tf_list[key]:
# Rotate horizontals
ftH = utils.rotate_dir(ft1, ft2, tfnoise.tilt)
TF_ZH = transfunc['ZH']['TF_ZH']
TF_PH = transfunc[key]['TF_PH']
TF_ZP_H = transfunc[key]['TF_ZP-H']
corrspec = ftZ - TF_ZH*ftH - (ftP - ftH*TF_PH)*TF_ZP_H
corrtime = np.real(np.fft.ifft(corrspec))
correct.add('ZP-H', corrtime)
self.correct = correct
[docs]
def save(self, filename):
"""
Method to save the object to file using `~Pickle`.
Parameters
----------
filename : str
File name
Examples
--------
Following from the example outlined in method
:func:`~obstools.atacr.classes.EventStream.correct_data`, we simply
save the final object
>>> evstream.save('evstream_demo.pkl')
Check that object has been saved
>>> import glob
>>> glob.glob("./evstream_demo.pkl")
['./evstream_demo.pkl']
"""
if not hasattr(self, 'correct'):
print("Warning: saving EventStream object before having done " +
"the corrections")
file = open(filename, 'wb')
pickle.dump(self, file)
file.close()