diff --git a/docs/atacr.rst b/docs/atacr.rst index 7339c9d..9fbb6bd 100644 --- a/docs/atacr.rst +++ b/docs/atacr.rst @@ -57,7 +57,7 @@ horizontal (``?H1,2``) component data. It is of course possible to combine both corrections and apply them sequentially. In this case the tilt noise is removed first, followed by compliance. This analysis requires all four components: three-component -seismic (``?HZ,1,2``) and pressure (``?XH``) data. +seismic (``?HZ,1,2``) and pressure (``?DH``) data. API documentation ***************** @@ -173,7 +173,7 @@ Utility functions Plotting functions ++++++++++++++++++ -.. automodule:: obstools.atacr.plot +.. automodule:: obstools.atacr.plotting :members: Scripts @@ -237,8 +237,8 @@ Usage -C CHANNELS, --channels CHANNELS Specify a comma-separated list of channels for which to perform the transfer function analysis. Possible - options are 'H' (for horizontal channels) or 'P' (for - pressure channel). Specifying 'H' allows for tilt + options are '12' (for horizontal channels) or 'P' (for + pressure channel). Specifying '12' allows for tilt correction. Specifying 'P' allows for compliance correction. [Default looks for both horizontal and pressure and allows for both tilt AND compliance @@ -566,8 +566,8 @@ Usage -C CHANNELS, --channels CHANNELS Specify a comma-separated list of channels for which to perform the transfer function analysis. Possible - options are 'H' (for horizontal channels) or 'P' (for - pressure channel). Specifying 'H' allows for tilt + options are '12' (for horizontal channels) or 'P' (for + pressure channel). Specifying '12' allows for tilt correction. Specifying 'P' allows for compliance correction. [Default looks for both horizontal and pressure and allows for both tilt AND compliance @@ -767,7 +767,7 @@ To check the station info for M08A, use the program ``ls_stdb``: We wish to download one month of data for the station M08A for March 2012. The various options above allow us to select the additional channels to use -(e.g., ``-C H,P`` for both horizontal and pressure data - which is the default +(e.g., ``-C 12,P`` for both horizontal and pressure data - which is the default setting). Default frequency settings for data pre-processing match those of the Matlab ``ATaCR`` software and can therefore be ignored when calling the program. Since the file ``M08A.pkl`` contains only one station, it is not necessary to specify a @@ -810,7 +810,7 @@ An example log printed on the terminal will look like: *********************************************************** * Downloading day-long data for key 7D.M08A and day 2012.61 * - * Channels selected: ['H', 'P'] and vertical + * Channels selected: ['12', 'P'] and vertical * 2012.061.*SAC * -> Downloading Seismic data... * ...done @@ -826,7 +826,7 @@ An example log printed on the terminal will look like: *********************************************************** * Downloading day-long data for key 7D.M08A and day 2012.62 * - * Channels selected: ['H', 'P'] and vertical + * Channels selected: ['12', 'P'] and vertical * 2012.062.*SAC * -> Downloading Seismic data... @@ -911,7 +911,7 @@ pass the quality control are colored red. .. figure:: ../obstools/examples/figures/Figure_3b.png :align: center -Figure 3a: Daily average PSD of bad (red) and good (black) windows. +Figure 3b: Daily average PSD of bad (red) and good (black) windows. 3. QC for clean station averages ++++++++++++++++++++++++++++++++ @@ -1025,12 +1025,12 @@ Once the ``StaNoise`` objects have been produced and saved to disk, the transfer functions across all available components can be calculated. By default the software will calculate the ones for which the spectral averages are available. -For compliance only (i.e., only ``?HZ`` and ``?XH?`` components are available), +For compliance only (i.e., only ``?HZ`` and ``?DH?`` components are available), the only transfer function possible is: - ``ZP`` -For tilt only (i.e., all of ``?HZ,1,2`` components are available, but not ``?XH``), +For tilt only (i.e., all of ``?HZ,1,2`` components are available, but not ``?DH``), the transfer functions are: - ``Z1`` @@ -1150,7 +1150,7 @@ Vanuatu earthquake (be conservative with the options), type in a terminal: * Dep: 33.70; Mag: 6.6 * M08A -> Ev: 9651.91 km; 86.80 deg; 239.43; 40.95 - * Channels selected: ['H', 'P'] and vertical + * Channels selected: ['12', 'P'] and vertical * 2012.069.07.09 * -> Downloading Seismic data... * ...done diff --git a/docs/comply.rst b/docs/comply.rst index ca9d362..95270ce 100644 --- a/docs/comply.rst +++ b/docs/comply.rst @@ -165,7 +165,7 @@ Once the ``StaNoise`` objects have been produced and saved to disk, the complian functions across all available components can be calculated. By default the software will calculate the ones for which the spectral averages are available. -If no horizontal components are available (i.e., only ``?HZ`` and ``?XH?`` components are available), +If no horizontal components are available (i.e., only ``?HZ`` and ``?DH?`` components are available), the only compliance function possible is: - ``ZP`` @@ -229,8 +229,10 @@ argument ``--f0``) .. figure:: ../obstools/examples/figures/Figure_comply.png :align: center -Figure 1: Coherence and compliance functions for the component combinations of interest, -as indicated in the title of each subplot. This example is for the month of March 2012 -for station M08A. The daily compliance functions are shown in grey and the average -calculated for the whole month is shown in black. +Figure 1: Coherence and compliance functions for the component combinations +of interest, as indicated in the title of each subplot. This example is for +the month of March 2012 for station M08A. The QC'ed daily compliance functions +are shown in grey. The QC'ed stations average calculated for the whole month is +shown in black. Note that *all* daily compliance functions are calculated and +displayed, even those that did not pass the QC threshold for the station average. diff --git a/obstools/atacr/classes.py b/obstools/atacr/classes.py index 7b1259d..b90a076 100644 --- a/obstools/atacr/classes.py +++ b/obstools/atacr/classes.py @@ -22,7 +22,7 @@ import sys -from scipy.signal import spectrogram, stft, detrend +from scipy.signal import stft, detrend from scipy.linalg import norm import matplotlib.pyplot as plt import numpy as np @@ -32,7 +32,7 @@ from pkg_resources import resource_filename from pathlib import Path np.seterr(all='ignore') -#np.set_printoptions(threshold=sys.maxsize) +# np.set_printoptions(threshold=sys.maxsize) class Power(object): @@ -145,6 +145,9 @@ class DayNoise(object): 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 @@ -169,22 +172,6 @@ class DayNoise(object): tf_list : Dict Dictionary of possible transfer functions given the available components. - goodwins : list - List of booleans representing whether a window is good (True) or - not (False). - This attribute is returned from the method - :func:`~obstools.atacr.classes.DayNoise.QC_daily_spectra` - power : :class:`~obstools.atacr.classes.Power` - Container for daily spectral power for all available components - cross : :class:`~obstools.atacr.classes.Cross` - Container for daily cross spectral power for all available components - rotation : :class:`~obstools.atacr.classes.Rotation` - Container for daily rotated (cross) spectral power for all available - components - f : :class:`~numpy.ndarray` - Frequency axis for corresponding time sampling parameters. Determined - from method - :func:`~obstools.atacr.classes.DayNoise.average_daily_spectra` Examples -------- @@ -198,10 +185,10 @@ class DayNoise(object): Now check its main attributes >>> print(*[daynoise.tr1, daynoise.tr2, daynoise.trZ, daynoise.trP], sep="\n") - 7D.M08A..1 | 2012-03-04T00:00:00.005500Z - 2012-03-04T23:59:59.805500Z | 5.0 Hz, 432000 samples - 7D.M08A..2 | 2012-03-04T00:00:00.005500Z - 2012-03-04T23:59:59.805500Z | 5.0 Hz, 432000 samples - 7D.M08A..P | 2012-03-04T00:00:00.005500Z - 2012-03-04T23:59:59.805500Z | 5.0 Hz, 432000 samples - 7D.M08A..Z | 2012-03-04T00:00:00.005500Z - 2012-03-04T23:59:59.805500Z | 5.0 Hz, 432000 samples + 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 @@ -227,7 +214,7 @@ def __init__(self, tr1=None, tr2=None, trZ=None, trP=None, window=7200., tr1 = st.select(component='1')[0] tr2 = st.select(component='2')[0] trZ = st.select(component='Z')[0] - trP = st.select(component='P')[0] + trP = st.select(component='H')[0] window = 7200. overlap = 0.3 key = '7D.M08A' @@ -248,11 +235,12 @@ def __init__(self, tr1=None, tr2=None, trZ=None, trP=None, window=7200., self.key = key # Get trace attributes - self.dt = self.trZ.stats.delta - self.npts = self.trZ.stats.npts - self.fs = self.trZ.stats.sampling_rate - self.year = self.trZ.stats.starttime.year - self.julday = self.trZ.stats.starttime.julday + 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 @@ -277,7 +265,7 @@ def __init__(self, tr1=None, tr2=None, trZ=None, trP=None, window=7200., self.av = False def QC_daily_spectra(self, pd=[0.004, 0.2], tol=1.5, alpha=0.05, - smooth=True, fig_QC=False, debug=False, save=False, + smooth=True, fig_QC=False, debug=False, save=None, form='png'): """ Method to determine daily time windows for which the spectra are @@ -301,9 +289,18 @@ def QC_daily_spectra(self, pd=[0.004, 0.2], tol=1.5, alpha=0.05, 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) @@ -342,23 +339,64 @@ def QC_daily_spectra(self, pd=[0.004, 0.2], tol=1.5, alpha=0.05, 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 - f, t, psdZ = spectrogram( - self.trZ.data, self.fs, window=wind, nperseg=ws, noverlap=ss) - self.f = f - print(f) + + # 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: - f, t, psdP = spectrogram( - self.trP.data, self.fs, window=wind, nperseg=ws, noverlap=ss) + psdP = np.abs(ftP)**2*2/self.dt + psdP = psdP[0:faxis, :] if self.ncomp == 3 or self.ncomp == 4: - f, t, psd1 = spectrogram( - self.tr1.data, self.fs, window=wind, nperseg=ws, noverlap=ss) - f, t, psd2 = spectrogram( - self.tr2.data, self.fs, window=wind, nperseg=ws, noverlap=ss) + 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: @@ -579,7 +617,7 @@ def QC_daily_spectra(self, pd=[0.004, 0.2], tol=1.5, alpha=0.05, self.QC = True def average_daily_spectra(self, calc_rotation=True, fig_average=False, - fig_coh_ph=False, save=False, form='png'): + fig_coh_ph=False, save=None, form='png'): """ Method to average the daily spectra for good windows. By default, the method will attempt to calculate the azimuth of maximum coherence @@ -597,6 +635,10 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False, fig_coh_ph : boolean Whether or not to produce a figure showing the maximum coherence between H and Z + save : :class:`~pathlib.Path` object + Relative path to figures folder + form : str + File format (e.g., 'png', 'jpg', 'eps') Attributes ---------- @@ -628,8 +670,8 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False, >>> daynoise.__dict__.keys() dict_keys(['tr1', 'tr2', 'trZ', 'trP', 'window', 'overlap', 'key', - 'dt', 'npts', 'fs', 'year', 'julday', 'ncomp', 'tf_list', 'QC', 'av', - 'goodwins', 'f', 'power', 'cross', 'rotation']) + 'dt', 'npts', 'fs', 'year', 'julday', 'tkey', 'ncomp', 'tf_list', + 'QC', 'av', 'f', 'goodwins', 'power', 'cross', 'rotation']) """ @@ -638,61 +680,29 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False, "QC_daily_spectra using default values") self.QC_daily_spectra() - # Points in window - ws = int(self.window/self.dt) - - # Number of points in step - 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] - - ft1 = None - ft2 = None - ftZ = None - ftP = None - - _f, _t, ftZ = stft( - self.trZ.data, self.fs, return_onesided=False, boundary=None, padded=False, - window=wind, nperseg=ws, noverlap=ss) - 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) - 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) - _f, _t, ft2 = stft( - self.tr2.data, self.fs, return_onesided=False, boundary=None, padded=False, - window=wind, nperseg=ws, noverlap=ss) - ft1 = ft1.T - ft2 = ft2.T - # Extract good windows c11 = None c22 = None cZZ = None cPP = None cZZ = np.abs( - np.mean(ftZ[self.goodwins, :]*np.conj(ftZ[self.goodwins, :]), - axis=0))[0:len(self.f)] + 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(ftP[self.goodwins, :]*np.conj(ftP[self.goodwins, :]), - axis=0))[0:len(self.f)] + 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(ft1[self.goodwins, :]*np.conj(ft1[self.goodwins, :]), - axis=0))[0:len(self.f)] + np.mean( + self.ft1[self.goodwins, :]*np.conj( + self.ft1[self.goodwins, :]), axis=0)) c22 = np.abs( - np.mean(ft2[self.goodwins, :]*np.conj(ft2[self.goodwins, :]), - axis=0))[0:len(self.f)] + np.mean( + self.ft2[self.goodwins, :]*np.conj( + self.ft2[self.goodwins, :]), axis=0)) # Extract bad windows bc11 = None @@ -700,20 +710,24 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False, bcZZ = None bcPP = None if np.sum(~self.goodwins) > 0: - bcZZ = np.abs(np.mean( - ftZ[~self.goodwins, :]*np.conj(ftZ[~self.goodwins, :]), - axis=0))[0:len(self.f)] + 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( - ftP[~self.goodwins, :]*np.conj(ftP[~self.goodwins, :]), - axis=0))[0:len(self.f)] + 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( - ft1[~self.goodwins, :]*np.conj(ft1[~self.goodwins, :]), - axis=0))[0:len(self.f)] - bc22 = np.abs(np.mean( - ft2[~self.goodwins, :]*np.conj(ft2[~self.goodwins, :]), - axis=0))[0:len(self.f)] + 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 @@ -723,20 +737,26 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False, c2P = None cZP = None if self.ncomp == 3 or self.ncomp == 4: - c12 = np.mean(ft1[self.goodwins, :] * - np.conj(ft2[self.goodwins, :]), axis=0)[0:len(self.f)] - c1Z = np.mean(ft1[self.goodwins, :] * - np.conj(ftZ[self.goodwins, :]), axis=0)[0:len(self.f)] - c2Z = np.mean(ft2[self.goodwins, :] * - np.conj(ftZ[self.goodwins, :]), axis=0)[0:len(self.f)] + 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(ft1[self.goodwins, :] * - np.conj(ftP[self.goodwins, :]), axis=0)[0:len(self.f)] - c2P = np.mean(ft2[self.goodwins, :] * - np.conj(ftP[self.goodwins, :]), axis=0)[0:len(self.f)] + 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(ftZ[self.goodwins, :] * - np.conj(ftP[self.goodwins, :]), axis=0)[0:len(self.f)] + 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) @@ -758,7 +778,8 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False, if calc_rotation and self.ncomp >= 3: cHH, cHZ, cHP, coh, ph, direc, tilt, coh_value, phase_value = \ utils.calculate_tilt( - ft1, ft2, ftZ, ftP, self.f, self.goodwins) + self.ft1, self.ft2, self.ftZ, self.ftP, self.f, + self.goodwins) self.rotation = Rotation( cHH, cHZ, cHP, coh, ph, tilt, coh_value, phase_value, direc) @@ -838,9 +859,11 @@ class StaNoise(object): 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 + `QC_sta_spectra()` method), each individual spectral quantity is unpacked as an object attribute and the original `DayNoise` objects are removed - from memory. In addition, all spectral quantities associated with the + 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. @@ -888,7 +911,7 @@ class StaNoise(object): , , ] - >>> sta.initialized + >>> stanoise.initialized False """ @@ -903,7 +926,7 @@ def _load_dn(day): tr1 = st.select(component='1')[0] tr2 = st.select(component='2')[0] trZ = st.select(component='Z')[0] - trP = st.select(component='P')[0] + trP = st.select(component='H')[0] window = 7200. overlap = 0.3 key = '7D.M08A' @@ -1024,15 +1047,13 @@ def init(self): Check that `daylist` attribute has been deleted >>> stanoise.daylist - --------------------------------------------------------------------------- - AttributeError Traceback (most recent call last) - in - ----> 1 stanoise.daylist + Traceback (most recent call last): + File "", line 1, in 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', 'direc', 'tilt', 'f', 'nwins', - 'ncomp', 'key', 'tf_list', 'QC', 'av']) + dict_keys(['initialized', 'c11', 'c22', 'cZZ', 'cPP', 'c12', 'c1Z', + 'c1P', 'c2Z', 'c2P', 'cZP', 'cHH', 'cHZ', 'cHP', 'direc', 'tilt', 'f', + 'nwins', 'ncomp', 'key', 'tf_list', 'QC', 'av']) """ @@ -1089,7 +1110,7 @@ def init(self): del self.daylist def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05, - fig_QC=False, debug=False, save=False, form='png'): + 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 @@ -1110,6 +1131,10 @@ def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05, 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 ---------- @@ -1121,7 +1146,7 @@ def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05, -------- Import demo data, call method and generate final figure - >>> obstools.atacr import StaNoise + >>> 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) @@ -1140,6 +1165,9 @@ def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05, # 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 @@ -1170,29 +1198,29 @@ def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05, if self.ncomp == 2: plt.figure(2) plt.subplot(2, 1, 1) - plt.semilogx(self.f, sl_cZZ, 'g', lw=0.5) + plt.semilogx(self.f[faxis], sl_cZZ[faxis], 'g', lw=0.5) plt.subplot(2, 1, 2) - plt.semilogx(self.f, sl_cPP, 'k', lw=0.5) + 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, sl_c11, 'r', lw=0.5) + plt.semilogx(self.f[faxis], sl_c11[faxis], 'r', lw=0.5) plt.subplot(3, 1, 2) - plt.semilogx(self.f, sl_c22, 'b', lw=0.5) + plt.semilogx(self.f[faxis], sl_c22[faxis], 'b', lw=0.5) plt.subplot(3, 1, 3) - plt.semilogx(self.f, sl_cZZ, 'g', lw=0.5) + 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, sl_c11, 'r', lw=0.5) + plt.semilogx(self.f[faxis], sl_c11[faxis], 'r', lw=0.5) plt.subplot(4, 1, 2) - plt.semilogx(self.f, sl_c22, 'b', lw=0.5) + plt.semilogx(self.f[faxis], sl_c22[faxis], 'b', lw=0.5) plt.subplot(4, 1, 3) - plt.semilogx(self.f, sl_cZZ, 'g', lw=0.5) + plt.semilogx(self.f[faxis], sl_cZZ[faxis], 'g', lw=0.5) plt.subplot(4, 1, 4) - plt.semilogx(self.f, sl_cPP, 'k', lw=0.5) + plt.semilogx(self.f[faxis], sl_cPP[faxis], 'k', lw=0.5) plt.tight_layout() if debug: plt.show() @@ -1261,7 +1289,7 @@ def QC_sta_spectra(self, pd=[0.004, 0.2], tol=2.0, alpha=0.05, else: plot.show() - def average_sta_spectra(self, fig_average=False, save=False, form='png'): + def average_sta_spectra(self, fig_average=False, save=None, form='png'): r""" Method to average the daily station spectra for good windows. @@ -1270,6 +1298,10 @@ def average_sta_spectra(self, fig_average=False, save=False, form='png'): 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 ---------- @@ -1739,48 +1771,36 @@ class EventStream(object): Note ---- - An ``EventStream`` object is defined as the data - (:class:`~obspy.core.Stream` object) are read from file or downloaded - from an ``obspy`` Client. Based on the available components, a list of + 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 ---------- - sta : :class:`~stdb.StdbElement` - An instance of an stdb object + 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 - sth : :class:`~obspy.core.Stream` - Stream containing three-component seismic data (traces are empty if - data are not available) - stp : :class:`~obspy.core.Stream` - Stream containing pressure data (trace is empty if data are not - available) + evtime : :class:`~obspy.core.UTCDateTime` + Origin time of seismic event. tstamp : str Time stamp for event - evlat : float - Latitude of seismic event - evlon : float - Longitude of seismic event - evtime : :class:`~obspy.core.UTCDateTime` - Origin time of seismic event - window : float - Length of time window in seconds + prefix : str + Associated prefix of SAC files + npts : int + Number of points in time series. fs : float - Sampling frequency (in Hz) + Sampling frequency (in Hz). dt : float - Sampling distance in seconds - npts : int - Number of points in time series + Sampling distance in seconds. ncomp : int - Number of available components (either 2, 3 or 4) + 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. - correct : :class:`~obstools.atacr.classes.EventStream.CorrectDict` - Container Dictionary for all possible corrections from the transfer - functions. This is calculated from the method - :func:`~obstools.atacr.classes.EventStream.correct_data` Examples -------- @@ -1791,63 +1811,59 @@ class EventStream(object): >>> evstream = EventStream('demo') Uploading demo earthquake data - March 09, 2012, station 7D.M08A >>> evstream.__dict__.keys() - dict_keys(['sta', 'key', 'sth', 'stp', 'tstamp', 'evlat', 'evlon', 'evtime', - 'window', 'fs', 'dt', 'ncomp', 'ev_list']) + dict_keys(['tr1', 'tr2', 'trZ', 'trP, 'key', 'evtime', + 'tstamp', 'prefix', 'npts', fs', 'dt', 'ncomp', 'ev_list']) Plot the raw traces - >>> import obstools.atacr.plot as plot - >>> plot.fig_event_raw(evstream, fmin=1./150., fmax=2.) + >>> 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, sta=None, sth=None, stp=None, tstamp=None, lat=None, - lon=None, time=None, window=None, sampling_rate=None, - ncomp=None, correct=False): + def __init__(self, tr1=Trace(), tr2=Trace(), trZ=Trace(), trP=Trace(), + correct=False): - if sta == 'demo' or sta == 'Demo': + if tr1 == 'demo': print("Uploading demo earthquake data - March 09, 2012, " + "station 7D.M08A") exmpl_path = Path(resource_filename('obstools', 'examples')) - fn = '2012.069.07.09.event.pkl' - fn = exmpl_path / 'event' / fn - evstream = pickle.load(open(fn, 'rb')) - sta = evstream.sta - key = evstream.key - sth = evstream.sth - stp = evstream.stp - tstamp = evstream.tstamp - lat = evstream.evlat - lon = evstream.evlon - time = evstream.evtime - window = evstream.window - sampling_rate = evstream.fs - ncomp = evstream.ncomp - correct = evstream.correct - - if any(value == None for value in [sta, sth, stp, tstamp, lat, lon, - time, window, sampling_rate, - ncomp]): + 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( - "Error: Initializing EventStream object with None values - " + - "aborting")) + "Incorrect initialization of EventStream object: " + + "missing a vertical component or too few components")) - self.sta = sta - self.key = sta.network+'.'+sta.station - self.sth = sth - self.stp = stp + 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.evlat = lat - self.evlon = lon - self.evtime = time - self.window = window - self.fs = sampling_rate - self.dt = 1./sampling_rate + self.prefix = self.key + '.' + self.tstamp + self.npts = zstats.npts + self.fs = zstats.sampling_rate + self.dt = zstats.delta self.ncomp = ncomp - self.correct = False # Build list of available transfer functions for future use if self.ncomp == 2: @@ -1908,8 +1924,9 @@ def correct_data(self, tfnoise): Plot the corrected traces - >>> import obstools.atacr.plot as plot - >>> plot.fig_event_corrected(evstream, tfnoise_day.tf_list) + >>> 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 @@ -1935,6 +1952,13 @@ def correct_data(self, tfnoise): .. 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: @@ -1948,99 +1972,75 @@ def correct_data(self, tfnoise): tf_list = tfnoise.tf_list transfunc = tfnoise.transfunc - # Points in window - ws = int(self.window/self.dt) - # Extract traces - trZ = Trace() - tr1 = Trace() - tr2 = Trace() - trP = Trace() - trZ = self.sth.select(component='Z')[0] - if self.ncomp == 2 or self.ncomp == 4: - trP = self.stp[0] - if self.ncomp == 3 or self.ncomp == 4: - tr1 = self.sth.select(component='1')[0] - tr2 = self.sth.select(component='2')[0] + trZ = self.trZ + tr1 = self.tr1 + tr2 = self.tr2 + trP = self.trP # Get Fourier spectra ft1 = None ft2 = None ftZ = None ftP = None - ftZ = np.fft.fft(trZ, n=ws) + + 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=ws) + ftP = np.fft.fft(trP, n=self.npts) if self.ncomp == 3 or self.ncomp == 4: - ft1 = np.fft.fft(tr1, n=ws) - ft2 = np.fft.fft(tr2, n=ws) + ft1 = np.fft.fft(tr1, n=self.npts) + ft2 = np.fft.fft(tr2, n=self.npts) - # Use one-sided frequency axis to match spectrogram - f = np.fft.rfftfreq(ws, d=self.dt) + 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)) + 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'] - fTF_ZP = np.hstack( - (TF_ZP, np.conj(TF_ZP[::-1][1:len(f)-1]))) - corrspec = ftZ - fTF_ZP*ftP - corrtime = np.real(np.fft.ifft(corrspec))[0:ws] + 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'] - fTF_Z1 = np.hstack( - (TF_Z1, np.conj(TF_Z1[::-1][1:len(f)-1]))) - corrspec = ftZ - fTF_Z1*ft1 - corrtime = np.real(np.fft.ifft(corrspec))[0:ws] + 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'] - fTF_Z1 = np.hstack( - (TF_Z1, np.conj(TF_Z1[::-1][1:len(f)-1]))) TF_21 = transfunc[key]['TF_21'] - fTF_21 = np.hstack( - (TF_21, np.conj(TF_21[::-1][1:len(f)-1]))) TF_Z2_1 = transfunc[key]['TF_Z2-1'] - fTF_Z2_1 = np.hstack( - (TF_Z2_1, np.conj(TF_Z2_1[::-1][1:len(f)-1]))) - corrspec = ftZ - fTF_Z1*ft1 - (ft2 - ft1*fTF_21)*fTF_Z2_1 - corrtime = np.real(np.fft.ifft(corrspec))[0:ws] + 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'] - fTF_Z1 = np.hstack( - (TF_Z1, np.conj(TF_Z1[::-1][1:len(f)-1]))) TF_21 = transfunc[key]['TF_21'] - fTF_21 = np.hstack( - (TF_21, np.conj(TF_21[::-1][1:len(f)-1]))) TF_Z2_1 = transfunc[key]['TF_Z2-1'] - fTF_Z2_1 = np.hstack( - (TF_Z2_1, np.conj(TF_Z2_1[::-1][1:len(f)-1]))) TF_P1 = transfunc[key]['TF_P1'] - fTF_P1 = np.hstack( - (TF_P1, np.conj(TF_P1[::-1][1:len(f)-1]))) TF_P2_1 = transfunc[key]['TF_P2-1'] - fTF_P2_1 = np.hstack( - (TF_P2_1, np.conj(TF_P2_1[::-1][1:len(f)-1]))) TF_ZP_21 = transfunc[key]['TF_ZP-21'] - fTF_ZP_21 = np.hstack( - (TF_ZP_21, np.conj(TF_ZP_21[::-1][1:len(f)-1]))) - corrspec = ftZ - fTF_Z1*ft1 - \ - (ft2 - ft1*fTF_21)*fTF_Z2_1 - \ - (ftP - ft1*fTF_P1 - - (ft2 - ft1*fTF_21)*fTF_P2_1)*fTF_ZP_21 - corrtime = np.real(np.fft.ifft(corrspec))[0:ws] + 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]: @@ -2050,10 +2050,8 @@ def correct_data(self, tfnoise): ftH = utils.rotate_dir(ft1, ft2, tfnoise.tilt) TF_ZH = transfunc[key]['TF_ZH'] - fTF_ZH = np.hstack( - (TF_ZH, np.conj(TF_ZH[::-1][1:len(f)-1]))) - corrspec = ftZ - fTF_ZH*ftH - corrtime = np.real(np.fft.ifft(corrspec))[0:ws] + 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]: @@ -2063,16 +2061,10 @@ def correct_data(self, tfnoise): ftH = utils.rotate_dir(ft1, ft2, tfnoise.tilt) TF_ZH = transfunc['ZH']['TF_ZH'] - fTF_ZH = np.hstack( - (TF_ZH, np.conj(TF_ZH[::-1][1:len(f)-1]))) TF_PH = transfunc[key]['TF_PH'] - fTF_PH = np.hstack( - (TF_PH, np.conj(TF_PH[::-1][1:len(f)-1]))) TF_ZP_H = transfunc[key]['TF_ZP-H'] - fTF_ZP_H = np.hstack( - (TF_ZP_H, np.conj(TF_ZP_H[::-1][1:len(f)-1]))) - corrspec = ftZ - fTF_ZH*ftH - (ftP - ftH*fTF_PH)*fTF_ZP_H - corrtime = np.real(np.fft.ifft(corrspec))[0:ws] + 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 @@ -2103,7 +2095,7 @@ def save(self, filename): """ - if not self.correct: + if hasattr(self, 'correct'): print("Warning: saving EventStream object before having done " + "the corrections") diff --git a/obstools/atacr/plotting.py b/obstools/atacr/plotting.py index 921bda0..bda305e 100644 --- a/obstools/atacr/plotting.py +++ b/obstools/atacr/plotting.py @@ -28,12 +28,13 @@ import numpy as np from matplotlib import pyplot as plt from obstools.atacr import utils +from obspy import Trace def fig_QC(f, power, gooddays, ncomp, key=''): """ - Function to plot the Quality-Control step of the analysis. This function is used - in both the `obs_daily_spectra.py` or `obs_clean_spectra.py` scripts. + Function to plot the Quality-Control step of the analysis. This function + is used in both the `atacr_daily_spectra` or `atacr_clean_spectra` scripts. Parameters ---------- @@ -42,7 +43,8 @@ def fig_QC(f, power, gooddays, ncomp, key=''): power : :class:`~obstools.classes.Power` Container for the Power spectra gooddays : List - List of booleans representing whether a window is good (True) or not (False) + List of booleans representing whether a window is good (True) or not + (False) ncomp : int Number of components used in analysis (can be 2, 3 or 4) key : str @@ -71,12 +73,15 @@ def fig_QC(f, power, gooddays, ncomp, key=''): 'HZ component, Station: '+key, 'HP component, Station: '+key] + # Extract only positive frequencies + faxis = f > 0 + fig = plt.figure(6) for i, sl in enumerate(sls): ax = fig.add_subplot(ncomp, 1, i+1) - ax.semilogx(f, sl[:, gooddays], 'k', lw=0.5) + ax.semilogx(f[faxis], sl[:, gooddays][faxis], 'k', lw=0.5) if np.sum(~gooddays) > 0: - plt.semilogx(f, sl[:, ~gooddays], 'r', lw=0.5) + plt.semilogx(f[faxis], sl[:, ~gooddays][faxis], 'r', lw=0.5) ax.set_title(title[i], fontdict={'fontsize': 8}) if i == len(sls)-1: plt.xlabel('Frequency (Hz)', fontdict={'fontsize': 8}) @@ -89,7 +94,7 @@ def fig_average(f, power, bad, gooddays, ncomp, key=''): """ Function to plot the averaged spectra (those qualified as 'good' in the QC step). This function is used - in both the `obs_daily_spectra.py` or `obs_clean_spectra.py` scripts. + in both the `atacr_daily_spectra` or `atacr_clean_spectra` scripts. Parameters ---------- @@ -100,7 +105,8 @@ def fig_average(f, power, bad, gooddays, ncomp, key=''): bad : :class:`~obstools.classes.Power` Container for the *bad* Power spectra gooddays : List - List of booleans representing whether a window is good (True) or not (False) + List of booleans representing whether a window is good (True) or not + (False) ncomp : int Number of components used in analysis (can be 2, 3 or 4) key : str @@ -136,12 +142,17 @@ def fig_average(f, power, bad, gooddays, ncomp, key=''): 'Average HZ, Station: '+key, 'Average HP, Station: '+key] + # Extract only positive frequencies + faxis = f > 0 + plt.figure() for i, (cc, bc) in enumerate(zip(ccs, bcs)): ax = plt.subplot(ncomp, 1, i+1) - ax.semilogx(f, utils.smooth(np.log(cc), 50), 'k', lw=0.5) + ax.semilogx( + f[faxis], utils.smooth(np.log(cc)[faxis], 50), 'k', lw=0.5) if np.sum(~gooddays) > 0: - ax.semilogx(f, utils.smooth(np.log(bc), 50), 'r', lw=0.5) + ax.semilogx( + f[faxis], utils.smooth(np.log(bc)[faxis], 50), 'r', lw=0.5) ax.set_title(title[i], fontdict={'fontsize': 8}) if i == len(ccs)-1: plt.xlabel('Frequency (Hz)', fontdict={'fontsize': 8}) @@ -153,8 +164,8 @@ def fig_average(f, power, bad, gooddays, ncomp, key=''): def fig_av_cross(f, field, gooddays, ftype, ncomp, key='', save=False, fname='', form='png', **kwargs): """ - Function to plot the averaged cross-spectra (those qualified as 'good' in the - QC step). This function is used in the `obs_daily_spectra.py` script. + Function to plot the averaged cross-spectra (those qualified as 'good' in + the QC step). This function is used in the `atacr_daily_spectra` script. Parameters ---------- @@ -163,9 +174,11 @@ def fig_av_cross(f, field, gooddays, ftype, ncomp, key='', field : :class:`~obstools.classes.Rotation` Container for the Power spectra gooddays : List - List of booleans representing whether a window is good (True) or not (False) + List of booleans representing whether a window is good (True) or not + (False) ftype : str - Type of plot to be displayed. If ftype is Admittance, plot is loglog. Otherwise semilogx + Type of plot to be displayed. If ftype is Admittance, plot is loglog. + Otherwise semilogx key : str String corresponding to the station key under analysis **kwargs : None @@ -173,6 +186,9 @@ def fig_av_cross(f, field, gooddays, ftype, ncomp, key='', """ + # Extract only positive frequencies + faxis = f > 0 + if ncomp == 2: fieldZP = field.cZP.T fields = [fieldZP] @@ -200,13 +216,17 @@ def fig_av_cross(f, field, gooddays, ftype, ncomp, key='', ax = fig.add_subplot(len(fields), 1, i+1) # Extact field if ftype == 'Admittance': - ax.loglog(f, field[:, gooddays], color='gray', **kwargs) + ax.loglog( + f[faxis], field[:, gooddays][faxis], color='gray', **kwargs) if np.sum(~gooddays) > 0: - ax.loglog(f, field[:, ~gooddays], color='r', **kwargs) + ax.loglog( + f[faxis], field[:, ~gooddays][faxis], color='r', **kwargs) else: - ax.semilogx(f, field[:, gooddays], color='gray', **kwargs) + ax.semilogx( + f[faxis], field[:, gooddays][faxis], color='gray', **kwargs) if np.sum(~gooddays) > 0: - ax.semilogx(f, field[:, ~gooddays], color='r', **kwargs) + ax.semilogx( + f[faxis], field[:, ~gooddays][faxis], color='r', **kwargs) plt.ylabel(ftype, fontdict={'fontsize': 8}) plt.title(key+' '+ftype+title[i], fontdict={'fontsize': 8}) if i == len(fields)-1: @@ -219,8 +239,8 @@ def fig_av_cross(f, field, gooddays, ftype, ncomp, key='', def fig_coh_ph(coh, ph, direc): """ - Function to plot the coherence and phase between the rotated H and Z components, - used to characterize the tilt direction. + Function to plot the coherence and phase between the rotated H and Z + components, used to characterize the tilt direction. Parameters ---------- @@ -269,15 +289,22 @@ def fig_TF(f, day_trfs, day_list, sta_trfs, sta_list, skey=''): Frequency axis (in Hz) day_trfs : Dict Dictionary containing the transfer functions for the daily averages + day_list : Dict + Dictionary containing the list of daily transfer functions sta_trfs : Dict Dictionary containing the transfer functions for the station averages - key : str + sta_list : Dict + Dictionary containing the list of average transfer functions + skey : str String corresponding to the station key under analysis """ import matplotlib.ticker as mtick + # Extract only positive frequencies + faxis = f > 0 + # Get max number of TFs to plot ntf = max(sum(day_list.values()), sum(sta_list.values())) @@ -301,9 +328,14 @@ def fig_TF(f, day_trfs, day_list, sta_trfs, sta_list, skey=''): if day_list[key]: for i in range(len(day_trfs)): ax.loglog( - f, np.abs(day_trfs[i][key]['TF_'+key]), 'gray', lw=0.5) + f[faxis], + np.abs(day_trfs[i][key]['TF_'+key][faxis]), + 'gray', lw=0.5) if sta_list[key]: - ax.loglog(f, np.abs(sta_trfs[key]['TF_'+key]), 'k', lw=0.5) + ax.loglog( + f[faxis], + np.abs(sta_trfs[key]['TF_'+key][faxis]), + 'k', lw=0.5) if key == 'ZP': ax.set_ylim(1.e-20, 1.e0) ax.set_xlim(1.e-4, 2.5) @@ -343,7 +375,8 @@ def fig_TF(f, day_trfs, day_list, sta_trfs, sta_list, skey=''): return plt -def fig_comply(f, day_comps, day_list, sta_comps, sta_list, sta, f_0): +def fig_comply(f, day_comps, day_list, sta_comps, sta_list, skey=None, + elev=-1000., f_0=None): """ Function to plot the transfer functions available. @@ -353,22 +386,32 @@ def fig_comply(f, day_comps, day_list, sta_comps, sta_list, sta, f_0): Frequency axis (in Hz) day_comps : Dict Dictionary containing the compliance functions for the daily averages + day_list : Dict + Dictionary containing the list of daily transfer functions sta_comps : Dict Dictionary containing the compliance functions for the station averages - key : str + sta_list : Dict + Dictionary containing the list of average transfer functions + skey : str String corresponding to the station key under analysis + elev : float + Station elevation in meters (OBS stations have negative elevations) + f_0 : float + Lowest frequency to consider in plot (Hz) """ import matplotlib.ticker as mtick import matplotlib.pyplot as plt - # Get station information - sta_key = sta.network + '.' + sta.station - sta_H = -1.*sta.elevation*1.e3 + # Extract only positive frequencies + faxis = f > 0 + + # Positive station elevation for frequency limit calc + elev = -1.*elev # Calculate theoretical frequency limit for infra-gravity waves - f_c = np.sqrt(9.81/np.pi/sta_H)/2. + f_c = np.sqrt(9.81/np.pi/elev)/2. # Define all possible combinations comp_list = {'ZP': True, 'ZP-21': True, 'ZP-H': True} @@ -395,26 +438,33 @@ def fig_comply(f, day_comps, day_list, sta_comps, sta_list, sta, f_0): if day_list[key]: for i in range(len(day_comps)): compliance = np.abs(day_comps[i][key][0]) - ax.plot(f, compliance, 'gray', alpha=0.3, lw=0.5) + ax.plot( + f[faxis], + compliance[faxis], + 'gray', alpha=0.3, lw=0.5) ax.set_xlim(f_0, f_c) ytop = np.max(compliance[(f > f_0) & (f < f_c)]) ybot = np.min(compliance[(f > f_0) & (f < f_c)]) ax.set_ylim(ybot, ytop) if sta_list[key]: - ax.plot(f, np.abs(sta_comps[key][0]), 'k', lw=0.5) + ax.plot( + f[faxis], + np.abs(sta_comps[key][0][faxis]), + 'k', lw=0.5) if key == 'ZP': - ax.set_title(sta_key+' Compliance: ZP', + ax.set_title(skey+' Compliance: ZP', fontdict={'fontsize': 8}) elif key == 'ZP-21': - ax.set_title(sta_key+' Compliance: ZP-21', + ax.set_title(skey+' Compliance: ZP-21', fontdict={'fontsize': 8}) elif key == 'ZP-H': - ax.set_title(sta_key+' Compliance: ZP-H', + ax.set_title(skey+' Compliance: ZP-H', fontdict={'fontsize': 8}) - ax.axvline(f_0, ls='--', c='k', lw=0.75) + if f_0: + ax.axvline(f_0, ls='--', c='k', lw=0.75) ax.axvline(f_c, ls='--', c='k', lw=0.75) ax = fig.add_subplot(ncomps, 2, j*2+2) @@ -423,22 +473,27 @@ def fig_comply(f, day_comps, day_list, sta_comps, sta_list, sta, f_0): if day_list[key]: for i in range(len(day_comps)): ax.semilogx( - f, np.abs(day_comps[i][key][1]), 'gray', - alpha=0.3, lw=0.5) + f[faxis], + np.abs(day_comps[i][key][1][faxis]), + 'gray', alpha=0.3, lw=0.5) if sta_list[key]: - ax.semilogx(f, np.abs(sta_comps[key][1]), 'k', lw=0.5) + ax.semilogx( + f[faxis], + np.abs(sta_comps[key][1][faxis]), + 'k', lw=0.5) if key == 'ZP': - ax.set_title(sta_key+' Coherence: ZP', + ax.set_title(skey+' Coherence: ZP', fontdict={'fontsize': 8}) elif key == 'ZP-21': - ax.set_title(sta_key+' Coherence: ZP-21', + ax.set_title(skey+' Coherence: ZP-21', fontdict={'fontsize': 8}) elif key == 'ZP-H': - ax.set_title(sta_key+' Coherence: ZP-H', + ax.set_title(skey+' Coherence: ZP-H', fontdict={'fontsize': 8}) - ax.axvline(f_0, ls='--', c='k', lw=0.75) + if f_0: + ax.axvline(f_0, ls='--', c='k', lw=0.75) ax.axvline(f_c, ls='--', c='k', lw=0.75) axes = plt.gcf().get_axes() @@ -450,7 +505,7 @@ def fig_comply(f, day_comps, day_list, sta_comps, sta_list, sta, f_0): return plt -def fig_event_raw(evstream, fmin, fmax): +def fig_event_raw(evstream, fmin=1./150., fmax=2.): """ Function to plot the raw (although bandpassed) seismograms. @@ -465,49 +520,55 @@ def fig_event_raw(evstream, fmin, fmax): """ - import matplotlib as mpl - evstream.sth.filter('bandpass', freqmin=fmin, - freqmax=fmax, corners=2, zerophase=True) - evstream.stp.filter('bandpass', freqmin=fmin, - freqmax=fmax, corners=2, zerophase=True) - sr = evstream.sth[0].stats.sampling_rate - taxis = np.arange(0., 7200., 1./sr) + from obspy import Stream + + # Unpack traces + tr1 = evstream.tr1.copy() + tr2 = evstream.tr2.copy() + trZ = evstream.trZ.copy() + trP = evstream.trP.copy() + st = Stream(traces=[tr for tr in [tr1, tr2, trZ, trP] if np.any(tr.data)]) + st.filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + sr = trZ.stats.sampling_rate + + taxis = np.arange(0., trZ.stats.npts/sr, 1./sr) fig = plt.figure(figsize=(6, 6)) ax = fig.add_subplot(4, 1, 1) - ax.plot(taxis, evstream.sth.select(component='Z')[0].data, 'k', lw=0.5) - ax.set_title(evstream.key+' '+evstream.tstamp + + ax.plot(taxis, trZ.data, 'k', lw=0.5) + ax.set_title(evstream.key + ' ' + evstream.tstamp + ': Z', fontdict={'fontsize': 8}) ax.ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - ax.set_xlim((0., 7200.)) + ax.set_xlim((0., trZ.stats.npts/sr)) - if len(evstream.sth) > 1: + if len(tr1.data > 0): ax = fig.add_subplot(4, 1, 2) - ax.plot(taxis, evstream.sth.select(component='1')[0].data, 'k', lw=0.5) + ax.plot(taxis, tr1.data, 'k', lw=0.5) ax.set_xlim((0., 7200.)) - ax.set_title(evstream.tstamp+': 1', fontdict={'fontsize': 8}) + ax.set_title(evstream.tstamp + ': 1', fontdict={'fontsize': 8}) ax.ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) ax = fig.add_subplot(4, 1, 3) - ax.plot(taxis, evstream.sth.select(component='2')[0].data, 'k', lw=0.5) - ax.set_xlim((0., 7200.)) - ax.set_title(evstream.tstamp+': 2', fontdict={'fontsize': 8}) + ax.plot(taxis, tr2.data, 'k', lw=0.5) + ax.set_xlim((0., trZ.stats.npts/sr)) + ax.set_title(evstream.tstamp + ': 2', fontdict={'fontsize': 8}) ax.ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - if evstream.stp: - if len(evstream.sth) > 1: + if len(trP.data > 0): + if len(tr1.data > 0): ax = fig.add_subplot(4, 1, 4) else: ax = fig.add_subplot(4, 1, 2) - ax.plot(taxis, evstream.stp[0].data, 'k', lw=0.5) + ax.plot(taxis, trP.data, 'k', lw=0.5) ax.ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - ax.set_xlim((0., 7200.)) - ax.set_title(evstream.tstamp+': P', fontdict={'fontsize': 8}) + ax.set_xlim((0., trZ.stats.npts/sr)) + ax.set_title(evstream.tstamp + ': P', fontdict={'fontsize': 8}) plt.xlabel('Time since earthquake (sec)') plt.tight_layout() @@ -515,7 +576,7 @@ def fig_event_raw(evstream, fmin, fmax): return plt -def fig_event_corrected(evstream, TF_list): +def fig_event_corrected(evstream, TF_list, fmin=1./150., fmax=2.): """ Function to plot the corrected vertical component seismograms. @@ -529,76 +590,99 @@ def fig_event_corrected(evstream, TF_list): """ - import matplotlib as mpl - evstream.sth.filter('bandpass', freqmin=1./150., - freqmax=1./10., corners=2, zerophase=True) - evstream.stp.filter('bandpass', freqmin=1./150., - freqmax=1./10., corners=2, zerophase=True) - sr = evstream.sth[0].stats.sampling_rate - taxis = np.arange(0., 7200., 1./sr) + # Unpack vertical trace and filter + trZ = evstream.trZ.copy() + trZ.filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + sr = trZ.stats.sampling_rate + taxis = np.arange(0., trZ.stats.npts/sr, 1./sr) plt.figure(figsize=(8, 8)) plt.subplot(611) plt.plot( - taxis, evstream.sth.select(component='Z')[0].data, 'lightgray', lw=0.5) + taxis, trZ.data, 'lightgray', lw=0.5) if TF_list['Z1']: - plt.plot(taxis, evstream.correct['Z1'], 'k', lw=0.5) - plt.title(evstream.key+' '+evstream.tstamp + + tr = Trace( + data=evstream.correct['Z1'], + header=trZ.stats).filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + plt.plot(taxis, tr.data, 'k', lw=0.5) + plt.title(evstream.key + ' ' + evstream.tstamp + ': Z1', fontdict={'fontsize': 8}) plt.gca().ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - plt.xlim((0., 7200.)) + plt.xlim((0., trZ.stats.npts/sr)) plt.subplot(612) plt.plot( - taxis, evstream.sth.select(component='Z')[0].data, 'lightgray', lw=0.5) + taxis, trZ.data, 'lightgray', lw=0.5) if TF_list['Z2-1']: - plt.plot(taxis, evstream.correct['Z2-1'], 'k', lw=0.5) - plt.title(evstream.tstamp+': Z2-1', fontdict={'fontsize': 8}) + tr = Trace( + data=evstream.correct['Z2-1'], + header=trZ.stats).filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + plt.plot(taxis, tr.data, 'k', lw=0.5) + plt.title(evstream.tstamp + ': Z2-1', fontdict={'fontsize': 8}) plt.gca().ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - plt.xlim((0., 7200.)) + plt.xlim((0., trZ.stats.npts/sr)) plt.subplot(613) plt.plot( - taxis, evstream.sth.select(component='Z')[0].data, 'lightgray', lw=0.5) + taxis, trZ.data, 'lightgray', lw=0.5) if TF_list['ZP-21']: - plt.plot(taxis, evstream.correct['ZP-21'], 'k', lw=0.5) - plt.title(evstream.tstamp+': ZP-21', fontdict={'fontsize': 8}) + tr = Trace( + data=evstream.correct['ZP-21'], + header=trZ.stats).filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + plt.plot(taxis, tr.data, 'k', lw=0.5) + plt.title(evstream.tstamp + ': ZP-21', fontdict={'fontsize': 8}) plt.gca().ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - plt.xlim((0., 7200.)) + plt.xlim((0., trZ.stats.npts/sr)) plt.subplot(614) plt.plot( - taxis, evstream.sth.select(component='Z')[0].data, 'lightgray', lw=0.5) + taxis, trZ.data, 'lightgray', lw=0.5) if TF_list['ZH']: - plt.plot(taxis, evstream.correct['ZH'], 'k', lw=0.5) - plt.title(evstream.tstamp+': ZH', fontdict={'fontsize': 8}) + tr = Trace( + data=evstream.correct['ZH'], + header=trZ.stats).filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + plt.plot(taxis, tr.data, 'k', lw=0.5) + plt.title(evstream.tstamp + ': ZH', fontdict={'fontsize': 8}) plt.gca().ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - plt.xlim((0., 7200.)) + plt.xlim((0., trZ.stats.npts/sr)) plt.subplot(615) plt.plot( - taxis, evstream.sth.select(component='Z')[0].data, 'lightgray', lw=0.5) + taxis, trZ.data, 'lightgray', lw=0.5) if TF_list['ZP-H']: - plt.plot(taxis, evstream.correct['ZP-H'], 'k', lw=0.5) - plt.title(evstream.tstamp+': ZP-H', fontdict={'fontsize': 8}) + tr = Trace( + data=evstream.correct['ZP-H'], + header=trZ.stats).filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + plt.plot(taxis, tr.data, 'k', lw=0.5) + plt.title(evstream.tstamp + ': ZP-H', fontdict={'fontsize': 8}) plt.gca().ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - plt.xlim((0., 7200.)) + plt.xlim((0., trZ.stats.npts/sr)) plt.subplot(616) plt.plot( - taxis, evstream.sth.select(component='Z')[0].data, 'lightgray', lw=0.5) + taxis, trZ.data, 'lightgray', lw=0.5) if TF_list['ZP']: - plt.plot(taxis, evstream.correct['ZP'], 'k', lw=0.5) - plt.title(evstream.tstamp+': ZP', fontdict={'fontsize': 8}) + tr = Trace( + data=evstream.correct['ZP'], + header=trZ.stats).filter( + 'bandpass', freqmin=fmin, freqmax=fmax, corners=2, zerophase=True) + plt.plot(taxis, tr.data, 'k', lw=0.5) + plt.title(evstream.tstamp + ': ZP', fontdict={'fontsize': 8}) plt.gca().ticklabel_format(axis='y', style='sci', useOffset=True, scilimits=(-3, 3)) - plt.xlim((0., 7200.)) + plt.xlim((0., trZ.stats.npts/sr)) plt.xlabel('Time since earthquake (sec)') plt.tight_layout() diff --git a/obstools/atacr/utils.py b/obstools/atacr/utils.py index 75dcfd2..59b47cd 100644 --- a/obstools/atacr/utils.py +++ b/obstools/atacr/utils.py @@ -31,19 +31,29 @@ class methods of `~obstools.atacr.classes`. import numpy as np import fnmatch from matplotlib import pyplot as plt -from obspy.core import read, Stream, Trace, AttribDict -from scipy.signal import savgol_filter - - -def floor_decimal(n, decimals=0): - multiplier = 10 ** decimals - return math.floor(n * multiplier) / multiplier +from obspy.core import read, Stream, Trace, AttribDict, UTCDateTime def traceshift(trace, tt): """ Function to shift traces in time given travel time + + Parameters + ---------- + + trace : :class:`~obspy.core.Trace` object + Trace object to update + tt : float + Time shift in seconds + + Returns + ------- + + rtrace : :class:`~obspy.core.Trace` object + Updated trace object + + """ # Define frequencies @@ -69,6 +79,32 @@ def traceshift(trace, tt): def QC_streams(start, end, st): + """ + Function for quality control of traces, which compares the + start and end times that were requested, as well as the total n + length of the traces. + + + Parameters + ---------- + + start : :class:`~obspy.core.UTCDateTime` object + Start time of requested stream + end : :class:`~obspy.core.UTCDateTime` object + End time of requested stream + st : :class:`~obspy.core.Stream` object + Stream object with all trace data + + Returns + ------- + + (pass): bool + Whether the QC test has passed + st : :class:`~obspy.core.Stream` object + Updated stream object + + + """ # Check start times if not np.all([tr.stats.starttime == start for tr in st]): @@ -83,19 +119,11 @@ def QC_streams(start, end, st): traces=[traceshift(tr, dt) for tr, dt in zip(st, delay)]) st = st_shifted.copy() - # # Check sampling rate - # sr = st[0].stats.sampling_rate - # sr_round = float(floor_decimal(sr, 0)) - # if not sr == sr_round: - # print("* Sampling rate is not an integer value: ", sr) - # print("* -> Resampling") - # st.resample(sr_round, no_filter=False) - # Try trimming dt = st[0].stats.delta try: st.trim(start, end-dt, fill_value=0., pad=True) - except: + except Exception: print("* Unable to trim") print("* -> Skipping") print("**************************************************") @@ -125,7 +153,7 @@ def QC_streams(start, end, st): return True, st -def update_stats(tr, stla, stlo, stel, cha): +def update_stats(tr, stla, stlo, stel, cha, evla=None, evlo=None): """ Function to include SAC metadata to :class:`~obspy.core.Trace` objects @@ -138,8 +166,14 @@ def update_stats(tr, stla, stlo, stel, cha): Latitude of station stlo : float Longitude of station + stel : float + Station elevation (m) cha : str Channel for component + evla : float, optional + Latitude of event + evlo : float, optional + Longitute of event Returns ------- @@ -155,6 +189,9 @@ def update_stats(tr, stla, stlo, stel, cha): tr.stats.sac.stel = stel tr.stats.sac.kcmpnm = cha tr.stats.channel = cha + if evla is not None and evlo is not None: + tr.stats.sac.evla = evla + tr.stats.sac.evlo = evlo return tr @@ -175,8 +212,8 @@ def get_data(datapath, tstart, tend): Returns ------- tr1, tr2, trZ, trP : :class:`~obspy.core.Trace` object - Corresponding trace objects for components H1, H2, HZ and HP. Returns empty traces - for missing components. + Corresponding trace objects for components H1, H2, HZ and HP. Returns + empty traces for missing components. """ @@ -244,7 +281,8 @@ def get_data(datapath, tstart, tend): def get_event(eventpath, tstart, tend): """ - Function to grab all available earthquake data given a path and data time range + Function to grab all available earthquake data given a path and data time + range Parameters ---------- @@ -258,45 +296,47 @@ def get_event(eventpath, tstart, tend): Returns ------- tr1, tr2, trZ, trP : :class:`~obspy.core.Trace` object - Corresponding trace objects for components H1, H2, HZ and HP. Returns empty traces - for missing components. + Corresponding trace objects for components H1, H2, HZ and HP. Returns + empty traces for missing components. """ + # Find out how many events from Z.SAC files + eventfiles = list(eventpath.glob('*Z.SAC')) + if not eventfiles: + raise(Exception("No event found in folder "+str(eventpath))) + + # Extract events from time stamps + prefix = [file.name.split('.') for file in eventfiles] + evstamp = [p[0]+'.'+p[1]+'.'+p[2]+'.'+p[3]+'.' for p in prefix] + evDateTime = [UTCDateTime(p[0]+'-'+p[1]+'T'+p[2]+":"+p[3]) for p in prefix] + # Define empty streams tr1 = Stream() tr2 = Stream() trZ = Stream() trP = Stream() - # Time iterator - t1 = tstart - - # Cycle through each day withing time range - while t1 < tend: - - # Time stamp used in file name - tstamp = str(t1.year).zfill(4)+'.'+str(t1.julday).zfill(3)+'.' - - # Cycle through directory and load files - p = eventpath.glob('*.*') - files = [x for x in p if x.is_file()] - for file in files: - if fnmatch.fnmatch(str(file), '*' + tstamp + '*1.SAC'): - tr = read(str(file)) - tr1.append(tr[0]) - elif fnmatch.fnmatch(str(file), '*' + tstamp + '*2.SAC'): - tr = read(str(file)) - tr2.append(tr[0]) - elif fnmatch.fnmatch(str(file), '*' + tstamp + '*Z.SAC'): - tr = read(str(file)) - trZ.append(tr[0]) - elif fnmatch.fnmatch(str(file), '*' + tstamp + '*H.SAC'): - tr = read(str(file)) - trP.append(tr[0]) - - # Increase increment - t1 += 3600.*24. + # Cycle over all available files in time range + for event, tstamp in zip(evDateTime, evstamp): + if event >= tstart and event <= tend: + + # Cycle through directory and load files + p = list(eventpath.glob('*.SAC')) + files = [x for x in p if x.is_file()] + for file in files: + if fnmatch.fnmatch(str(file), '*' + tstamp + '*1.SAC'): + tr = read(str(file)) + tr1.append(tr[0]) + elif fnmatch.fnmatch(str(file), '*' + tstamp + '*2.SAC'): + tr = read(str(file)) + tr2.append(tr[0]) + elif fnmatch.fnmatch(str(file), '*' + tstamp + '*Z.SAC'): + tr = read(str(file)) + trZ.append(tr[0]) + elif fnmatch.fnmatch(str(file), '*' + tstamp + '*H.SAC'): + tr = read(str(file)) + trP.append(tr[0]) # Fill with empty traces if components are not found ntr = len(trZ) @@ -336,20 +376,24 @@ def calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreq=[0.005, 0.035]): f : :class:`~numpy.ndarray` Frequency axis in Hz goodwins : list - List of booleans representing whether a window is good (True) or not (False). - This attribute is returned from the method :func:`~obstools.atacr.classes.DayNoise.QC_daily_spectra` - tiltfreq : list - Two floats representing the frequency band at which the tilt is calculated + List of booleans representing whether a window is good (True) or not + (False). This attribute is returned from the method + :func:`~obstools.atacr.classes.DayNoise.QC_daily_spectra` + tiltfreq : list, optional + Two floats representing the frequency band at which the tilt is + calculated Returns ------- cHH, cHZ, cHP : :class:`~numpy.ndarray` - Arrays of power and cross-spectral density functions of components HH (rotated H1 - in direction of maximum tilt), HZ, and HP + Arrays of power and cross-spectral density functions of components HH + (rotated H1 in direction of maximum tilt), HZ, and HP coh : :class:`~numpy.ndarray` - Coherence value between rotated H and Z components, as a function of directions (azimuths) + Coherence value between rotated H and Z components, as a function of + directions (azimuths) ph : :class:`~numpy.ndarray` - Phase value between rotated H and Z components, as a function of directions (azimuths) + Phase value between rotated H and Z components, as a function of + directions (azimuths) direc : :class:`~numpy.ndarray` Array of directions (azimuths) considered tilt : float @@ -448,45 +492,6 @@ def calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreq=[0.005, 0.035]): return cHH, cHZ, cHP, coh, ph, direc, tilt, coh_value, phase_value -def calculate_windowed_fft(trace, ws, ss=None, hann=True): - """ - Calculates windowed Fourier transform - - Parameters - ---------- - trace : :class:`~obspy.core.Trace` - Input trace data - ws : int - Window size, in number of samples - ss : int - Step size, or number of samples until next window - han : bool - Whether or not to apply a Hanning taper to data - - Returns - ------- - ft : :class:`~numpy.ndarray` - Fourier transform of trace - f : :class:`~numpy.ndarray` - Frequency axis in Hz - - """ - - n2 = _npow2(ws) - f = trace.stats.sampling_rate/2. * np.linspace(0., 1., int(n2/2) + 1) - - # Extract sliding windows - tr, nd = sliding_window(trace.data, ws, ss) - - # Fourier transform - ft = np.fft.fft(tr, n=n2) - - return ft, f - -# def smooth(data, np, poly=0, axis=0): -# return savgol_filter(data, np, poly, axis=axis, mode='wrap') - - def smooth(data, nd, axis=0): """ Function to smooth power spectral density functions from the convolution @@ -498,7 +503,7 @@ def smooth(data, nd, axis=0): Real-valued array to smooth (PSD) nd : int Number of samples over which to smooth - axis : int + axis : int, optional axis over which to perform the smoothing Returns @@ -596,57 +601,6 @@ def phase(Gxy): return None -def sliding_window(a, ws, ss=None, hann=True): - """ - Function to split a data array into overlapping, possibly tapered sub-windows - - Parameters - ---------- - a : :class:`~numpy.ndarray` - 1D array of data to split - ws : int - Window size in samples - ss : int - Step size in samples. If not provided, window and step size - are equal. - - Returns - ------- - out : :class:`~numpy.ndarray` - 1D array of windowed data - nd : int - Number of windows - - """ - - if ss is None: - # no step size was provided. Return non-overlapping windows - ss = ws - - # Calculate the number of windows to return, ignoring leftover samples, and - # allocate memory to contain the samples - valid = len(a) - ss - nd = (valid) // ss - out = np.ndarray((nd, ws), dtype=a.dtype) - - if nd == 0: - if hann: - out = a * np.hanning(ws) - else: - out = a - - for i in range(nd): - # "slide" the window along the samples - start = i * ss - stop = start + ws - if hann: - out[i] = a[start: stop] * np.hanning(ws) - else: - out[i] = a[start: stop] - - return out, nd - - def rotate_dir(tr1, tr2, direc): d = -direc*np.pi/180.+np.pi/2. @@ -679,7 +633,3 @@ def ftest(res1, pars1, res2, pars2): P = 1. - (f_dist.cdf(Fobs, dof1, dof2) - f_dist.cdf(1./Fobs, dof1, dof2)) return P - - -def _npow2(x): - return 1 if x == 0 else 2**(x-1).bit_length() diff --git a/obstools/comply/classes.py b/obstools/comply/classes.py index 1133b2e..df0c510 100644 --- a/obstools/comply/classes.py +++ b/obstools/comply/classes.py @@ -42,8 +42,8 @@ class Comply(object): Attributes ---------- - sta : :class:`~stdb.StdbElement` - An instance of an stdb object + elev : float + Station elevation in meters (OBS stations have negative elevations) f : :class:`~numpy.ndarray` Frequency axis for corresponding time sampling parameters c11 : `numpy.ndarray` @@ -61,9 +61,9 @@ class Comply(object): """ - def __init__(self, sta=None, objnoise=None): + def __init__(self, objnoise=None, elev=None): - if any(value == None for value in [sta, objnoise]): + if any(value is None for value in [elev, objnoise]): raise(Exception( "Error: Initializing EventStream object with None values - " + "aborting")) @@ -79,7 +79,7 @@ def __init__(self, sta=None, objnoise=None): "Error: Noise object has not been processed (QC and " + "averaging) - aborting")) - self.sta = sta + self.elevation = elev self.f = objnoise.f self.c11 = objnoise.power.c11 self.c22 = objnoise.power.c22 @@ -118,7 +118,9 @@ def calculate_compliance(self): Examples -------- - Calculate compliance and coherence functions for a DayNoise object + Calculate compliance and coherence functions for a DayNoise object. + In these examples, station elevation is extracted from the IRIS + metadata aggregator site: http://ds.iris.edu/mda/7D/M08A/ >>> from obstools.atacr import DayNoise >>> from obstools.comply import Comply @@ -126,7 +128,7 @@ def calculate_compliance(self): Uploading demo data - March 04, 2012, station 7D.M08A >>> daynoise.QC_daily_spectra() >>> daynoise.average_daily_spectra() - >>> daycomply = Comply(objnoise=daynoise, sta=sta) + >>> daycomply = Comply(objnoise=daynoise, elev=-126.4) >>> daycomply.calculate_compliance() >>> tfnoise.complyfunc.keys() dict_keys(['ZP', 'ZP-21', 'ZP-H']) @@ -139,7 +141,7 @@ def calculate_compliance(self): Uploading demo data - March 01 to 04, 2012, station 7D.M08A >>> stanoise.QC_sta_spectra() >>> stanoise.average_sta_spectra() - >>> stacomply = Comply(objnoise=stanoise, sta=sta) + >>> stacomply = Comply(objnoise=stanoise, elev=-126.4) >>> stacomply.calculate_compliance() >>> stacomply.complyfunc.keys() dict_keys(['ZP', 'ZP-21']) @@ -189,7 +191,7 @@ def wavenumber(omega, H): k[i] = 0. else: - a0 = -27 * om**2 / g # constant terms + a0 = -27 * om**2 / g # constant terms a1 = 0. # no linear terms a2 = 27 * H - (9 * om**2 * H**2)/g # quadratic terms a3 = 0. # no cubic terms @@ -210,7 +212,7 @@ def wavenumber(omega, H): return k # Calculate wavenumber - careful here, elevation is negative - k = wavenumber(2.*np.pi*self.f, -1.*self.sta.elevation*1.e3) + k = wavenumber(2.*np.pi*self.f, -1.*self.elevation) # Initialize empty dictionary complyfunc = self.ComplyDict() diff --git a/obstools/examples/data/2012.061..BDH.SAC b/obstools/examples/data/2012.061..BDH.SAC new file mode 100644 index 0000000..07c75c2 Binary files /dev/null and b/obstools/examples/data/2012.061..BDH.SAC differ diff --git a/obstools/examples/data/2012.061..BH1.SAC b/obstools/examples/data/2012.061..BH1.SAC index f56cfc8..958c552 100644 Binary files a/obstools/examples/data/2012.061..BH1.SAC and b/obstools/examples/data/2012.061..BH1.SAC differ diff --git a/obstools/examples/data/2012.061..BH2.SAC b/obstools/examples/data/2012.061..BH2.SAC index 5ad3ad9..b10e732 100644 Binary files a/obstools/examples/data/2012.061..BH2.SAC and b/obstools/examples/data/2012.061..BH2.SAC differ diff --git a/obstools/examples/data/2012.061..BHH.SAC b/obstools/examples/data/2012.061..BHH.SAC deleted file mode 100644 index d3ddd8b..0000000 Binary files a/obstools/examples/data/2012.061..BHH.SAC and /dev/null differ diff --git a/obstools/examples/data/2012.061..BHZ.SAC b/obstools/examples/data/2012.061..BHZ.SAC index b3c18f7..664e84a 100644 Binary files a/obstools/examples/data/2012.061..BHZ.SAC and b/obstools/examples/data/2012.061..BHZ.SAC differ diff --git a/obstools/examples/data/2012.062..BDH.SAC b/obstools/examples/data/2012.062..BDH.SAC new file mode 100644 index 0000000..3dbd41b Binary files /dev/null and b/obstools/examples/data/2012.062..BDH.SAC differ diff --git a/obstools/examples/data/2012.062..BH1.SAC b/obstools/examples/data/2012.062..BH1.SAC index bea85b8..08881db 100644 Binary files a/obstools/examples/data/2012.062..BH1.SAC and b/obstools/examples/data/2012.062..BH1.SAC differ diff --git a/obstools/examples/data/2012.062..BH2.SAC b/obstools/examples/data/2012.062..BH2.SAC index fbe6fc9..825e2f5 100644 Binary files a/obstools/examples/data/2012.062..BH2.SAC and b/obstools/examples/data/2012.062..BH2.SAC differ diff --git a/obstools/examples/data/2012.062..BHH.SAC b/obstools/examples/data/2012.062..BHH.SAC deleted file mode 100644 index bb42d7c..0000000 Binary files a/obstools/examples/data/2012.062..BHH.SAC and /dev/null differ diff --git a/obstools/examples/data/2012.062..BHZ.SAC b/obstools/examples/data/2012.062..BHZ.SAC index 6562996..acd3a62 100644 Binary files a/obstools/examples/data/2012.062..BHZ.SAC and b/obstools/examples/data/2012.062..BHZ.SAC differ diff --git a/obstools/examples/data/2012.063..BDH.SAC b/obstools/examples/data/2012.063..BDH.SAC new file mode 100644 index 0000000..2b30c2a Binary files /dev/null and b/obstools/examples/data/2012.063..BDH.SAC differ diff --git a/obstools/examples/data/2012.063..BH1.SAC b/obstools/examples/data/2012.063..BH1.SAC index f6df5b9..1104066 100644 Binary files a/obstools/examples/data/2012.063..BH1.SAC and b/obstools/examples/data/2012.063..BH1.SAC differ diff --git a/obstools/examples/data/2012.063..BH2.SAC b/obstools/examples/data/2012.063..BH2.SAC index 514ee92..68e9a85 100644 Binary files a/obstools/examples/data/2012.063..BH2.SAC and b/obstools/examples/data/2012.063..BH2.SAC differ diff --git a/obstools/examples/data/2012.063..BHH.SAC b/obstools/examples/data/2012.063..BHH.SAC deleted file mode 100644 index af3a97b..0000000 Binary files a/obstools/examples/data/2012.063..BHH.SAC and /dev/null differ diff --git a/obstools/examples/data/2012.063..BHZ.SAC b/obstools/examples/data/2012.063..BHZ.SAC index 3ac5dbc..ec968a8 100644 Binary files a/obstools/examples/data/2012.063..BHZ.SAC and b/obstools/examples/data/2012.063..BHZ.SAC differ diff --git a/obstools/examples/data/2012.064..BDH.SAC b/obstools/examples/data/2012.064..BDH.SAC new file mode 100644 index 0000000..0f38e6f Binary files /dev/null and b/obstools/examples/data/2012.064..BDH.SAC differ diff --git a/obstools/examples/data/2012.064..BH1.SAC b/obstools/examples/data/2012.064..BH1.SAC index 25bc2e9..f6bcce1 100644 Binary files a/obstools/examples/data/2012.064..BH1.SAC and b/obstools/examples/data/2012.064..BH1.SAC differ diff --git a/obstools/examples/data/2012.064..BH2.SAC b/obstools/examples/data/2012.064..BH2.SAC index 6e21865..5e5e95d 100644 Binary files a/obstools/examples/data/2012.064..BH2.SAC and b/obstools/examples/data/2012.064..BH2.SAC differ diff --git a/obstools/examples/data/2012.064..BHH.SAC b/obstools/examples/data/2012.064..BHH.SAC deleted file mode 100644 index 14a8d20..0000000 Binary files a/obstools/examples/data/2012.064..BHH.SAC and /dev/null differ diff --git a/obstools/examples/data/2012.064..BHZ.SAC b/obstools/examples/data/2012.064..BHZ.SAC index 4f198e1..3e72408 100644 Binary files a/obstools/examples/data/2012.064..BHZ.SAC and b/obstools/examples/data/2012.064..BHZ.SAC differ diff --git a/obstools/examples/event/2012.069.07.09.BDH.SAC b/obstools/examples/event/2012.069.07.09.BDH.SAC new file mode 100644 index 0000000..c030bfb Binary files /dev/null and b/obstools/examples/event/2012.069.07.09.BDH.SAC differ diff --git a/obstools/examples/event/2012.069.07.09.BH1.SAC b/obstools/examples/event/2012.069.07.09.BH1.SAC new file mode 100644 index 0000000..96ee5a3 Binary files /dev/null and b/obstools/examples/event/2012.069.07.09.BH1.SAC differ diff --git a/obstools/examples/event/2012.069.07.09.BH2.SAC b/obstools/examples/event/2012.069.07.09.BH2.SAC new file mode 100644 index 0000000..6d339a1 Binary files /dev/null and b/obstools/examples/event/2012.069.07.09.BH2.SAC differ diff --git a/obstools/examples/event/2012.069.07.09.BHZ.SAC b/obstools/examples/event/2012.069.07.09.BHZ.SAC new file mode 100644 index 0000000..660bba2 Binary files /dev/null and b/obstools/examples/event/2012.069.07.09.BHZ.SAC differ diff --git a/obstools/examples/event/2012.069.07.09.event.pkl b/obstools/examples/event/2012.069.07.09.event.pkl deleted file mode 100644 index c05d6f5..0000000 Binary files a/obstools/examples/event/2012.069.07.09.event.pkl and /dev/null differ diff --git a/obstools/examples/event/2012.069.1.SAC b/obstools/examples/event/2012.069.1.SAC deleted file mode 100644 index 0e3847e..0000000 Binary files a/obstools/examples/event/2012.069.1.SAC and /dev/null differ diff --git a/obstools/examples/event/2012.069.2.SAC b/obstools/examples/event/2012.069.2.SAC deleted file mode 100644 index a338af6..0000000 Binary files a/obstools/examples/event/2012.069.2.SAC and /dev/null differ diff --git a/obstools/examples/event/2012.069.H.SAC b/obstools/examples/event/2012.069.H.SAC deleted file mode 100644 index c6c8c64..0000000 Binary files a/obstools/examples/event/2012.069.H.SAC and /dev/null differ diff --git a/obstools/examples/event/2012.069.Z.SAC b/obstools/examples/event/2012.069.Z.SAC deleted file mode 100644 index d0f0e60..0000000 Binary files a/obstools/examples/event/2012.069.Z.SAC and /dev/null differ diff --git a/obstools/examples/figures/Figure_10.png b/obstools/examples/figures/Figure_10.png index 807a903..1121f62 100644 Binary files a/obstools/examples/figures/Figure_10.png and b/obstools/examples/figures/Figure_10.png differ diff --git a/obstools/examples/figures/Figure_11.png b/obstools/examples/figures/Figure_11.png index b172e5c..cb6421a 100644 Binary files a/obstools/examples/figures/Figure_11.png and b/obstools/examples/figures/Figure_11.png differ diff --git a/obstools/examples/figures/Figure_12a.png b/obstools/examples/figures/Figure_12a.png index 3e4a292..c4fc78e 100644 Binary files a/obstools/examples/figures/Figure_12a.png and b/obstools/examples/figures/Figure_12a.png differ diff --git a/obstools/examples/figures/Figure_12b.png b/obstools/examples/figures/Figure_12b.png index ac1be18..209fb4f 100644 Binary files a/obstools/examples/figures/Figure_12b.png and b/obstools/examples/figures/Figure_12b.png differ diff --git a/obstools/examples/figures/Figure_2.png b/obstools/examples/figures/Figure_2.png index 14fa5f4..97cf785 100644 Binary files a/obstools/examples/figures/Figure_2.png and b/obstools/examples/figures/Figure_2.png differ diff --git a/obstools/examples/figures/Figure_3a.png b/obstools/examples/figures/Figure_3a.png index 304e59d..1dc0111 100644 Binary files a/obstools/examples/figures/Figure_3a.png and b/obstools/examples/figures/Figure_3a.png differ diff --git a/obstools/examples/figures/Figure_3b.png b/obstools/examples/figures/Figure_3b.png index d783245..32d800c 100644 Binary files a/obstools/examples/figures/Figure_3b.png and b/obstools/examples/figures/Figure_3b.png differ diff --git a/obstools/examples/figures/Figure_4.png b/obstools/examples/figures/Figure_4.png index 696ae3d..45b4cb1 100644 Binary files a/obstools/examples/figures/Figure_4.png and b/obstools/examples/figures/Figure_4.png differ diff --git a/obstools/examples/figures/Figure_6.png b/obstools/examples/figures/Figure_6.png index d7e42c3..79fc0bd 100644 Binary files a/obstools/examples/figures/Figure_6.png and b/obstools/examples/figures/Figure_6.png differ diff --git a/obstools/examples/figures/Figure_7.png b/obstools/examples/figures/Figure_7.png index 427bcf3..ed3fdee 100644 Binary files a/obstools/examples/figures/Figure_7.png and b/obstools/examples/figures/Figure_7.png differ diff --git a/obstools/examples/figures/Figure_8.png b/obstools/examples/figures/Figure_8.png index 8a68d47..ec71a5d 100644 Binary files a/obstools/examples/figures/Figure_8.png and b/obstools/examples/figures/Figure_8.png differ diff --git a/obstools/examples/figures/Figure_9.png b/obstools/examples/figures/Figure_9.png index 75faf05..9fed72c 100644 Binary files a/obstools/examples/figures/Figure_9.png and b/obstools/examples/figures/Figure_9.png differ diff --git a/obstools/examples/figures/Figure_comply.png b/obstools/examples/figures/Figure_comply.png index 7375f83..0ff8ecb 100644 Binary files a/obstools/examples/figures/Figure_comply.png and b/obstools/examples/figures/Figure_comply.png differ diff --git a/obstools/scripts/atacr_clean_spectra.py b/obstools/scripts/atacr_clean_spectra.py index 4441fc1..f71983d 100644 --- a/obstools/scripts/atacr_clean_spectra.py +++ b/obstools/scripts/atacr_clean_spectra.py @@ -36,11 +36,13 @@ from obspy import UTCDateTime from numpy import nan + def get_cleanspec_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. - Calling options for the script `obs_clean_spectra.py` that accompany this package. + Calling options for the script `obs_clean_spectra.py` that accompany this + package. """ @@ -208,7 +210,7 @@ def get_cleanspec_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from start time: " + args.startT) @@ -219,7 +221,7 @@ def get_cleanspec_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from end time: " + args.endT) @@ -251,7 +253,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop @@ -276,7 +278,9 @@ def main(args=None): # Path where spectra are located specpath = Path('SPECTRA') / stkey if not specpath.is_dir(): - raise(Exception("Path to "+str(specpath)+" doesn`t exist - aborting")) + raise(Exception( + "Path to " + str(specpath) + + " doesn`t exist - aborting")) # Path where average spectra will be saved avstpath = Path('AVG_STA') / stkey @@ -441,32 +445,32 @@ def main(args=None): try: ph_12_all.append( 180./np.pi*utils.phase(daynoise.cross.c12)) - except: + except Exception: ph_12_all.append(None) try: ph_1Z_all.append( 180./np.pi*utils.phase(daynoise.cross.c1Z)) - except: + except Exception: ph_1Z_all.append(None) try: ph_1P_all.append( 180./np.pi*utils.phase(daynoise.cross.c1P)) - except: + except Exception: ph_1P_all.append(None) try: ph_2Z_all.append( 180./np.pi*utils.phase(daynoise.cross.c2Z)) - except: + except Exception: ph_2Z_all.append(None) try: ph_2P_all.append( 180./np.pi*utils.phase(daynoise.cross.c2P)) - except: + except Exception: ph_2P_all.append(None) try: ph_ZP_all.append( 180./np.pi*utils.phase(daynoise.cross.cZP)) - except: + except Exception: ph_ZP_all.append(None) # Admittance @@ -527,32 +531,38 @@ def main(args=None): if args.fig_av_cross: fname = stkey + '.' + 'av_coherence' - plot = plotting.fig_av_cross(stanoise.f, coh, stanoise.gooddays, - 'Coherence', stanoise.ncomp, key=stkey, lw=0.5) + plot = plotting.fig_av_cross( + stanoise.f, coh, stanoise.gooddays, + 'Coherence', stanoise.ncomp, key=stkey, lw=0.5) # if plotpath.is_dir(): if plotpath: - plot.savefig(str(plotpath / (fname + '.' + args.form)), - dpi=300, bbox_inches='tight', format=args.form) + plot.savefig( + str(plotpath / (fname + '.' + args.form)), + dpi=300, bbox_inches='tight', format=args.form) else: plot.show() fname = stkey + '.' + 'av_admittance' - plot = plotting.fig_av_cross(stanoise.f, ad, stanoise.gooddays, - 'Admittance', stanoise.ncomp, key=stkey, lw=0.5) + plot = plotting.fig_av_cross( + stanoise.f, ad, stanoise.gooddays, + 'Admittance', stanoise.ncomp, key=stkey, lw=0.5) if plotpath: - plot.savefig(str(plotpath / (fname + '.' + args.form)), - dpi=300, bbox_inches='tight', format=args.form) + plot.savefig( + str(plotpath / (fname + '.' + args.form)), + dpi=300, bbox_inches='tight', format=args.form) else: plot.show() fname = stkey + '.' + 'av_phase' - plot = plotting.fig_av_cross(stanoise.f, ph, stanoise.gooddays, - 'Phase', stanoise.ncomp, key=stkey, marker=',', lw=0) + plot = plotting.fig_av_cross( + stanoise.f, ph, stanoise.gooddays, + 'Phase', stanoise.ncomp, key=stkey, marker=',', lw=0) if plotpath: - plot.savefig(str(plotpath / (fname + '.' + args.form)), - dpi=300, bbox_inches='tight', format=args.form) + plot.savefig( + str(plotpath / (fname + '.' + args.form)), + dpi=300, bbox_inches='tight', format=args.form) else: plot.show() @@ -560,8 +570,9 @@ def main(args=None): fname = stkey + '.' + 'coh_ph' plot = plotting.fig_coh_ph(coh_all, ph_all, stanoise.direc) if plotpath: - plot.savefig(str(plotpath / (fname + '.' + args.form)), - dpi=300, bbox_inches='tight', format=args.form) + plot.savefig( + str(plotpath / (fname + '.' + args.form)), + dpi=300, bbox_inches='tight', format=args.form) else: plot.show() diff --git a/obstools/scripts/atacr_correct_event.py b/obstools/scripts/atacr_correct_event.py index 2bcfb76..2303691 100644 --- a/obstools/scripts/atacr_correct_event.py +++ b/obstools/scripts/atacr_correct_event.py @@ -25,10 +25,10 @@ # Import modules and functions import numpy as np -from obspy import UTCDateTime +from obspy import UTCDateTime, Stream import pickle import stdb -from obstools.atacr import StaNoise, Power, Cross, Rotation, TFNoise +from obstools.atacr import EventStream from obstools.atacr import utils, plotting from pathlib import Path @@ -41,7 +41,8 @@ def get_correct_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. - Calling options for the script `obs_correct_event.py` that accompany this package. + Calling options for the script `obs_correct_event.py` that accompanies this + package. """ @@ -201,7 +202,7 @@ def get_correct_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from " + "start time: " + args.startT) @@ -212,7 +213,7 @@ def get_correct_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from " + "end time: " + args.endT) @@ -238,7 +239,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop @@ -325,48 +326,33 @@ def main(args=None): sta.enddate.strftime("%Y-%m-%d %H:%M:%S"))) print("|-----------------------------------------------|") - # Find all files in directories - p = eventpath.glob('*.pkl') - event_files = [x for x in p if x.is_file()] - p = transpath.glob('*.*') + # Get all components + trE1, trE2, trEZ, trEP = utils.get_event(eventpath, tstart, tend) + + # Find all TF files in directory + p = list(transpath.glob('*.*')) trans_files = [x for x in p if x.is_file()] # Check if folders contain anything - if not event_files: - raise(Exception("There are no events in folder " + str(eventpath))) - if not trans_files: raise(Exception("There are no transfer functions in folder " + str(transpath))) - # Cycle through available files - for eventfile in event_files: - - # Skip hidden files and folders - if eventfile.name[0] == '.': - continue + # Cycle through available data + for tr1, tr2, trZ, trP in zip(trE1, trE2, trEZ, trEP): - evprefix = eventfile.name.split('.') - evstamp = evprefix[0]+'.'+evprefix[1]+'.' - - evDateTime = UTCDateTime(evprefix[0]+'-'+evprefix[1]) - if evDateTime >= tstart and evDateTime <= tend: - - # Load event file - try: - file = open(eventfile, 'rb') - eventstream = pickle.load(file) - file.close() - except: - print("File "+str(eventfile) + - " exists but cannot be loaded") - continue + eventstream = EventStream(tr1, tr2, trZ, trP) - else: - continue + # Check if Trace is from SAC file with event info + evlo = None + evla = None + if hasattr(trZ.stats, 'sac'): + if hasattr(trZ.stats.sac, 'evlo'): + evlo = trZ.stats.sac.evlo + evla = trZ.stats.sac.evla if args.fig_event_raw: - fname = stkey + '.' + evstamp + 'raw' + fname = stkey + '.' + eventstream.tstamp + 'raw' plot = plotting.fig_event_raw( eventstream, fmin=args.fmin, fmax=args.fmax) @@ -386,17 +372,20 @@ def main(args=None): continue tfprefix = transfile.name.split('transfunc')[0] + print(tfprefix) # This case refers to the "cleaned" spectral averages if len(tfprefix) > 9: if not args.skip_clean: + yr1 = tfprefix.split('-')[0].split('.')[0] jd1 = tfprefix.split('-')[0].split('.')[1] yr2 = tfprefix.split('-')[1].split('.')[0] jd2 = tfprefix.split('-')[1].split('.')[1] date1 = UTCDateTime(yr1+'-'+jd1) date2 = UTCDateTime(yr2+'-'+jd2) - dateev = UTCDateTime(evprefix[0]+'-'+evprefix[1]) + dateev = eventstream.evtime + if dateev >= date1 and dateev <= date2: print(str(transfile) + " file found - applying transfer functions") @@ -405,7 +394,7 @@ def main(args=None): file = open(transfile, 'rb') tfaverage = pickle.load(file) file.close() - except: + except Exception: print("File "+str(transfile) + " exists but cannot be loaded") continue @@ -416,7 +405,7 @@ def main(args=None): correct_sta = eventstream.correct if args.fig_plot_corrected: - fname = stkey + '.' + evstamp + 'sta_corrected' + fname = eventstream.prefix + '.sta_corrected' plot = plotting.fig_event_corrected( eventstream, tfaverage.tf_list) # Save or show figure @@ -432,27 +421,37 @@ def main(args=None): correctpath = eventpath / 'CORRECTED' if not correctpath.is_dir(): correctpath.mkdir(parents=True) - file = correctpath / eventfile.stem - eventstream.save(str(file) + '.day.pkl') + file = correctpath / eventstream.prefix + eventstream.save(str(file) + '.sta.pkl') # Now save as SAC files for key, value in tfaverage.tf_list.items(): if value and eventstream.ev_list[key]: - nameZ = '.sta.' + key + '.' - nameZ += sta.channel+'Z.SAC' - fileZ = correctpath / (eventfile.stem + nameZ) - trZ = eventstream.sth.select(component='Z')[0].copy() - trZ.data = eventstream.correct[key] - trZ = utils.update_stats(trZ, - sta.latitude, sta.longitude, - sta.elevation, 'Z') - trZ.write(str(fileZ), format='SAC') + # Postfix + nameZ = '.sta.' + key + '.' + nameZ += sta.channel + 'Z.SAC' + + # Add Prefix and Postfix + fileZ = str(file) + nameZ + + # Select Z component and update trace + trZ = eventstream.trZ.copy() + trZ.data = eventstream.correct[key] + trZ = utils.update_stats( + trZ, sta.latitude, sta.longitude, + sta.elevation, sta.channel+'Z', + evla=evla, + evlo=evlo) + # Save as SAC file + trZ.write(str(fileZ), format='SAC') # This case refers to the "daily" spectral averages else: if not args.skip_daily: + evprefix = eventstream.tstamp.split('.') + evstamp = evprefix[0]+'.'+evprefix[1]+'.' if tfprefix == evstamp: print(str(transfile) + " file found - applying transfer functions") @@ -461,7 +460,7 @@ def main(args=None): file = open(transfile, 'rb') tfaverage = pickle.load(file) file.close() - except: + except Exception: print("File "+str(transfile) + " exists but cannot be loaded") continue @@ -472,7 +471,7 @@ def main(args=None): correct_day = eventstream.correct if args.fig_plot_corrected: - fname = stkey + '.' + evstamp + 'day_corrected' + fname = eventstream.prefix + '.day_corrected' plot = plotting.fig_event_corrected( eventstream, tfaverage.tf_list) # Save or show figure @@ -488,20 +487,30 @@ def main(args=None): correctpath = eventpath / 'CORRECTED' if not correctpath.is_dir(): correctpath.mkdir(parents=True) - file = correctpath / eventfile.stem - eventstream.save(str(file) + '.sta.pkl') + file = correctpath / eventstream.prefix + eventstream.save(str(file) + '.day.pkl') # Now save as SAC files for key, value in tfaverage.tf_list.items(): if value and eventstream.ev_list[key]: - nameZ = '.day.' + key + '.' + + # Postfix + nameZ = '.day.' + key + '.' nameZ += sta.channel+'Z.SAC' - fileZ = correctpath / (eventfile.stem + nameZ) - trZ = eventstream.sth.select(component='Z')[0].copy() + + # Add Prefix and Postfix + fileZ = str(file) + nameZ + + # Select Z component and update trace + trZ = eventstream.trZ.copy() trZ.data = eventstream.correct[key] - trZ = utils.update_stats(trZ, - sta.latitude, sta.longitude, - sta.elevation, 'Z') + trZ = utils.update_stats( + trZ, sta.latitude, sta.longitude, + sta.elevation, sta.channel+'Z', + evla=evla, + evlo=evlo) + + # Save as SAC file trZ.write(str(fileZ), format='SAC') diff --git a/obstools/scripts/atacr_daily_spectra.py b/obstools/scripts/atacr_daily_spectra.py index 0920ceb..2e7b6f4 100644 --- a/obstools/scripts/atacr_daily_spectra.py +++ b/obstools/scripts/atacr_daily_spectra.py @@ -41,7 +41,8 @@ def get_dailyspec_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. - Calling options for the script `obs_daily_spectra.py` that accompany this package. + Calling options for the script `obs_daily_spectra.py` that accompany this + package. """ @@ -252,7 +253,7 @@ def get_dailyspec_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from start time: " + args.startT) @@ -263,7 +264,7 @@ def get_dailyspec_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from end time: " + args.endT) @@ -296,7 +297,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop diff --git a/obstools/scripts/atacr_download_data.py b/obstools/scripts/atacr_download_data.py index 35cb25d..076ba10 100644 --- a/obstools/scripts/atacr_download_data.py +++ b/obstools/scripts/atacr_download_data.py @@ -42,7 +42,8 @@ def get_daylong_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. - Calling options for the script `obs_download_data.py` that accompany this package. + Calling options for the script `obs_download_data.py` that accompany this + package. """ @@ -82,10 +83,10 @@ def get_daylong_arguments(argv=None): default="", help="Specify a comma-separated list of channels for " + "which to perform the transfer function analysis. " + - "Possible options are 'H' (for horizontal channels) or 'P' " + - "(for pressure channel). Specifying 'H' allows " + + "Possible options are '12' (for horizontal channels 1 and 2) " + + "and/or 'P' (for pressure channel). Specifying '12' allows " + "for tilt correction. Specifying 'P' allows for compliance " + - "correction. [Default looks for both horizontal and " + + "correction. [Default '12,P' looks for both horizontal and " + "pressure and allows for both tilt AND compliance corrections]") parser.add_argument( "-O", "--overwrite", @@ -214,16 +215,16 @@ def get_daylong_arguments(argv=None): if len(args.channels) > 0: args.channels = args.channels.split(',') else: - args.channels = ["H", "P"] + args.channels = ["12", "P"] for cha in args.channels: - if cha not in ["H", "P"]: + if cha not in ["12", "P"]: parser.error("Error: Channel not recognized " + str(cha)) # construct start time if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from start time: " + args.startT) @@ -234,7 +235,7 @@ def get_daylong_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from end time: " + args.endT) @@ -293,7 +294,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop @@ -432,7 +433,7 @@ def main(args=None): starttime=t1, endtime=t2, attach_response=True) print("* ...done") - except: + except Exception: print(" Error: Unable to download ?H? components - " + "continuing") t1 += dt @@ -441,7 +442,7 @@ def main(args=None): st = sth - elif "H" not in args.channels: + elif "12" not in args.channels: # If data files exist, continue if fileZ.exists() and fileP.exists(): @@ -467,7 +468,7 @@ def main(args=None): starttime=t1, endtime=t2, attach_response=True) print("* ...done") - except: + except Exception: print(" Error: Unable to download ?H? components - " + "continuing") t1 += dt @@ -483,14 +484,16 @@ def main(args=None): if len(stp) > 1: print("WARNING: There are more than one ?DH trace") print("* -> Keeping the highest sampling rate") - print("* -> Renaming channel to "+sta.channel[0]+"DH") + print( + "* -> Renaming channel to " + + sta.channel[0]+"DH") if stp[0].stats.sampling_rate > \ stp[1].stats.sampling_rate: stp = Stream(traces=stp[0]) else: stp = Stream(traces=stp[1]) - except: + except Exception: print(" Error: Unable to download ?DH component - " + "continuing") t1 += dt @@ -527,7 +530,7 @@ def main(args=None): starttime=t1, endtime=t2, attach_response=True) print("* ...done") - except: + except Exception: print(" Error: Unable to download ?H? components - " + "continuing") t1 += dt @@ -543,14 +546,16 @@ def main(args=None): if len(stp) > 1: print("WARNING: There are more than one ?DH trace") print("* -> Keeping the highest sampling rate") - print("* -> Renaming channel to "+sta.channel[0]+"DH") + print( + "* -> Renaming channel to " + + sta.channel[0]+"DH") if stp[0].stats.sampling_rate > \ stp[1].stats.sampling_rate: stp = Stream(traces=stp[0]) else: stp = Stream(traces=stp[1]) - except: + except Exception: print(" Error: Unable to download ?DH component - " + "continuing") t1 += dt @@ -562,7 +567,8 @@ def main(args=None): # Detrend, filter st.detrend('demean') st.detrend('linear') - st.filter('lowpass', freq=0.5*args.new_sampling_rate, + st.filter( + 'lowpass', freq=0.5*args.new_sampling_rate, corners=2, zerophase=True) st.resample(args.new_sampling_rate) @@ -581,19 +587,19 @@ def main(args=None): # Extract traces - Z trZ = sth.select(component='Z')[0] trZ = utils.update_stats( - trZ, sta.latitude, sta.longitude, sta.elevation, + trZ, sta.latitude, sta.longitude, sta.elevation, sta.channel+'Z') trZ.write(str(fileZ), format='SAC') # Extract traces - H - if "H" in args.channels: + if "12" in args.channels: tr1 = sth.select(component='1')[0] tr2 = sth.select(component='2')[0] tr1 = utils.update_stats( - tr1, sta.latitude, sta.longitude, sta.elevation, + tr1, sta.latitude, sta.longitude, sta.elevation, sta.channel+'1') tr2 = utils.update_stats( - tr2, sta.latitude, sta.longitude, sta.elevation, + tr2, sta.latitude, sta.longitude, sta.elevation, sta.channel+'2') tr1.write(str(file1), format='SAC') tr2.write(str(file2), format='SAC') @@ -605,8 +611,8 @@ def main(args=None): stp.remove_response(pre_filt=args.pre_filt) trP = stp[0] trP = utils.update_stats( - trP, sta.latitude, sta.longitude, sta.elevation, - sta.channel[0]+'P') + trP, sta.latitude, sta.longitude, sta.elevation, + sta.channel[0]+'DH') trP.write(str(fileP), format='SAC') t1 += dt diff --git a/obstools/scripts/atacr_download_event.py b/obstools/scripts/atacr_download_event.py index 984b5e5..4e8e1b4 100644 --- a/obstools/scripts/atacr_download_event.py +++ b/obstools/scripts/atacr_download_event.py @@ -86,10 +86,10 @@ def get_event_arguments(argv=None): default="", help="Specify a comma-separated list of channels for " + "which to perform the transfer function analysis. " + - "Possible options are 'H' (for horizontal channels) or 'P' " + - "(for pressure channel). Specifying 'H' allows " + + "Possible options are '12' (for horizontal channels) or 'P' " + + "(for pressure channel). Specifying '12' allows " + "for tilt correction. Specifying 'P' allows for compliance " + - "correction. [Default looks for both horizontal and " + + "correction. [Default '12,P' looks for both horizontal and " + "pressure and allows for both tilt AND compliance corrections]") parser.add_argument( "-O", "--overwrite", @@ -123,32 +123,6 @@ def get_event_arguments(argv=None): "(--User-Auth='username:authpassword') to access and download " + "restricted data. [Default no user and password]") - """ - # Database Settings - DataGroup = parser.add_argument_group( - title="Local Data Settings", - description="Settings associated with defining " + - "and using a local data base of pre-downloaded day-long SAC files.") - DataGroup.add_argument( - "--local-data", - action="store", - type=str, - dest="localdata", - default=None, - help="Specify a comma separated list of paths containing " + - "day-long sac files of data already downloaded. " + - "If data exists for a seismogram is already present on disk, " + - "it is selected preferentially over downloading " + - "the data using the Client interface") - DataGroup.add_argument( - "--no-data-zero", - action="store_true", - dest="ndval", - default=False, - help="Specify to force missing data to be set as zero, " + - "rather than default behaviour which sets to nan.") -""" - # Constants Settings FreqGroup = parser.add_argument_group( title='Frequency Settings', @@ -178,6 +152,16 @@ def get_event_arguments(argv=None): help="Specify four comma-separated corner " + "frequencies (float, in Hz) for deconvolution " + "pre-filter. [Default 0.001,0.005,45.,50.]") + FreqGroup.add_argument( + "--window", + action="store", + type=float, + dest="window", + default=7200., + help="Specify window length in seconds. " + + "Default value is highly recommended. " + "Program may not be stable for large deviations " + + "from default value. [Default 7200. (or 2 hours)]") # Event Selection Criteria EventGroup = parser.add_argument_group( @@ -272,17 +256,17 @@ def get_event_arguments(argv=None): if len(args.channels) > 0: args.channels = args.channels.split(',') else: - args.channels = ['H', 'P'] + args.channels = ['12', 'P'] for cha in args.channels: - if cha not in ['H', 'P']: + if cha not in ['12', 'P']: parser.error("Error: Channel not recognized " + str(cha)) # construct start time if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from start time: " + args.startT) @@ -293,7 +277,7 @@ def get_event_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from end time: " + args.endT) @@ -312,18 +296,6 @@ def get_event_arguments(argv=None): else: args.UserAuth = [] - # # Parse Local Data directories - # if args.localdata is not None: - # args.localdata = args.localdata.split(',') - # else: - # args.localdata = [] - - # # Check NoData Value - # if args.ndval: - # args.ndval = 0.0 - # else: - # args.ndval = nan - if args.pre_filt is None: args.pre_filt = [0.001, 0.005, 45., 50.] else: @@ -349,7 +321,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop @@ -461,9 +433,6 @@ def main(args=None): # Extract event ev = cat[iev] - window = 7200. - new_sampling_rate = 5. - time = ev.origins[0].time dep = ev.origins[0].depth lon = ev.origins[0].longitude @@ -492,12 +461,13 @@ def main(args=None): # If distance outside of distance range: if not (gac > args.mindist and gac < args.maxdist): - print("\n* -> Event outside epicentral distance " + + print( + "\n* -> Event outside epicentral distance " + "range - continuing") continue t1 = time - t2 = t1 + window + t2 = t1 + args.window # Time stamp tstamp = str(time.year).zfill(4)+'.' + \ @@ -547,7 +517,7 @@ def main(args=None): location=sta.location[0], channel=channels, starttime=t1, endtime=t2, attach_response=True) print("* ...done") - except: + except Exception: print( " Error: Unable to download ?H? components - " + "continuing") @@ -555,7 +525,7 @@ def main(args=None): st = sth - elif "H" not in args.channels: + elif "12" not in args.channels: # Number of channels ncomp = 2 @@ -573,7 +543,7 @@ def main(args=None): location=sta.location[0], channel=channels, starttime=t1, endtime=t2, attach_response=True) print("* ...done") - except: + except Exception: print( " Error: Unable to download ?H? components - " + "continuing") @@ -588,14 +558,15 @@ def main(args=None): if len(stp) > 1: print("WARNING: There are more than one ?DH trace") print("* -> Keeping the highest sampling rate") - print("* -> Renaming channel to "+ - sta.channel[0]+"DH") + print( + "* -> Renaming channel to " + + sta.channel[0] + "DH") if stp[0].stats.sampling_rate > \ stp[1].stats.sampling_rate: stp = Stream(traces=stp[0]) else: stp = Stream(traces=stp[1]) - except: + except Exception: print(" Error: Unable to download ?DH component - " + "continuing") continue @@ -622,7 +593,7 @@ def main(args=None): location=sta.location[0], channel=channels, starttime=t1, endtime=t2, attach_response=True) print("* ...done") - except: + except Exception: print( " Error: Unable to download ?H? components - " + "continuing") @@ -637,14 +608,15 @@ def main(args=None): if len(stp) > 1: print("WARNING: There are more than one ?DH trace") print("* -> Keeping the highest sampling rate") - print("* -> Renaming channel to "+ - sta.channel[0]+"DH") + print( + "* -> Renaming channel to " + + sta.channel[0] + "DH") if stp[0].stats.sampling_rate > \ stp[1].stats.sampling_rate: stp = Stream(traces=stp[0]) else: stp = Stream(traces=stp[1]) - except: + except Exception: print(" Error: Unable to download ?DH component - " + "continuing") continue @@ -673,21 +645,21 @@ def main(args=None): # Extract traces - Z trZ = sth.select(component='Z')[0] trZ = utils.update_stats( - trZ, sta.latitude, sta.longitude, sta.elevation, - sta.channel+'Z') + trZ, sta.latitude, sta.longitude, sta.elevation, + sta.channel+'Z', evla=lat, evlo=lon) trZ.write(str(fileZ), format='SAC') # Extract traces and write out in SAC format # Seismic channels - if "H" in args.channels: + if "12" in args.channels: tr1 = sth.select(component='1')[0] tr2 = sth.select(component='2')[0] tr1 = utils.update_stats( - tr1, sta.latitude, sta.longitude, sta.elevation, - sta.channel+'1') + tr1, sta.latitude, sta.longitude, sta.elevation, + sta.channel+'1', evla=lat, evlo=lon) tr2 = utils.update_stats( - tr2, sta.latitude, sta.longitude, sta.elevation, - sta.channel+'2') + tr2, sta.latitude, sta.longitude, sta.elevation, + sta.channel+'2', evla=lat, evlo=lon) tr1.write(str(file1), format='SAC') tr2.write(str(file2), format='SAC') @@ -698,19 +670,16 @@ def main(args=None): stp.remove_response(pre_filt=args.pre_filt) trP = stp[0] trP = utils.update_stats( - trP, sta.latitude, sta.longitude, sta.elevation, - sta.channel[0]+'DH') + trP, sta.latitude, sta.longitude, sta.elevation, + sta.channel[0]+'DH', evla=lat, evlo=lon) trP.write(str(fileP), format='SAC') else: stp = Stream() - # Write out EventStream object - eventstream = EventStream( - sta, sth, stp, tstamp, lat, lon, time, window, - args.new_sampling_rate, ncomp) - - eventstream.save(filename) + # # Write out EventStream object + # eventstream = EventStream(sta, sth, stp) + # eventstream.save(filename) if __name__ == "__main__": diff --git a/obstools/scripts/atacr_transfer_functions.py b/obstools/scripts/atacr_transfer_functions.py index 63d653c..aef620b 100644 --- a/obstools/scripts/atacr_transfer_functions.py +++ b/obstools/scripts/atacr_transfer_functions.py @@ -36,11 +36,13 @@ from os.path import exists as exist from numpy import nan + def get_transfer_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. - Calling options for the script `obs_transfer functions.py` that accompany this package. + Calling options for the script `obs_transfer_functions.py` that accompanies + this package. """ @@ -174,7 +176,7 @@ def get_transfer_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from " + "start time: " + args.startT) @@ -185,7 +187,7 @@ def get_transfer_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from " + "end time: " + args.endT) @@ -211,7 +213,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop @@ -395,17 +397,18 @@ def main(args=None): if args.fig_TF: fname = stkey + '.' + 'transfer_functions' - plot = plotting.fig_TF(f, day_transfer_functions, daynoise.tf_list, - sta_transfer_functions, stanoise.tf_list, skey=stkey) + plot = plotting.fig_TF( + f, day_transfer_functions, daynoise.tf_list, + sta_transfer_functions, stanoise.tf_list, skey=stkey) if plotpath: - plot.savefig(plotpath / (fname + '.' + args.form), - dpi=300, bbox_inches='tight', format=args.form) + plot.savefig( + plotpath / (fname + '.' + args.form), + dpi=300, bbox_inches='tight', format=args.form) else: plot.show() - if __name__ == "__main__": # Run main program diff --git a/obstools/scripts/comply_calculate.py b/obstools/scripts/comply_calculate.py index dcb6678..93cb091 100644 --- a/obstools/scripts/comply_calculate.py +++ b/obstools/scripts/comply_calculate.py @@ -41,7 +41,8 @@ def get_comply_arguments(argv=None): """ Get Options from :class:`~optparse.OptionParser` objects. - Calling options for the script `obs_transfer functions.py` that accompany this package. + Calling options for the script `obs_transfer functions.py` that accompanies + this package. """ @@ -190,7 +191,7 @@ def get_comply_arguments(argv=None): if len(args.startT) > 0: try: args.startT = UTCDateTime(args.startT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from " + "start time: " + args.startT) @@ -201,7 +202,7 @@ def get_comply_arguments(argv=None): if len(args.endT) > 0: try: args.endT = UTCDateTime(args.endT) - except: + except Exception: parser.error( "Error: Cannot construct UTCDateTime from " + "end time: " + args.endT) @@ -231,7 +232,7 @@ def main(args=None): db, stkeys = stdb.io.load_db(fname=args.indb, keys=args.stkeys) # stdb=0.1.3 - except: + except Exception: db = stdb.io.load_db(fname=args.indb) # Construct station key loop @@ -378,7 +379,7 @@ def main(args=None): file.close() # Load spectra into TFNoise object - daycomply = Comply(objnoise=daynoise, sta=sta) + daycomply = Comply(objnoise=daynoise, elev=sta.elevation*1.e3) # Calculate the transfer functions daycomply.calculate_compliance() @@ -386,7 +387,7 @@ def main(args=None): # Store the frequency axis f = daycomply.f - # Append to list of transfer functions + # Append to list of transfer functions - for plotting day_comply_functions.append(daycomply.complyfunc) # Save daily transfer functions to file @@ -405,7 +406,7 @@ def main(args=None): if filepkl.exists(): if not args.ovr: print("* -> file " + str(filepkl) + - " exists - continuing") + " exists - loading") # Load Comply object and append to list stacomply = pickle.load(open(filepkl, 'rb')) @@ -424,7 +425,7 @@ def main(args=None): # Load spectra into TFNoise object - no Rotation object # for station averages rotation = Rotation(None, None, None) - stacomply = Comply(objnoise=stanoise, sta=sta) + stacomply = Comply(objnoise=stanoise, elev=sta.elevation*1.e3) # Calculate the transfer functions stacomply.calculate_compliance() @@ -432,7 +433,7 @@ def main(args=None): # Store the frequency axis f = stacomply.f - # Extract the transfer functions + # Extract the transfer functions - for plotting sta_comply_functions = stacomply.complyfunc # Save average transfer functions to file @@ -442,8 +443,8 @@ def main(args=None): fname = stkey + '.' + 'compliance' plot = plotting.fig_comply( f, day_comply_functions, daycomply.tf_list, - sta_comply_functions, stacomply.tf_list, sta=sta, - f_0=args.f0) + sta_comply_functions, stacomply.tf_list, skey=stkey, + elev=sta.elevation*1.e3, f_0=args.f0) if plotpath: plot.savefig(plotpath / (fname + '.' + args.form), diff --git a/obstools/tests/test_0_imports.py b/obstools/tests/test_0_imports.py index 5ea6953..de96b02 100644 --- a/obstools/tests/test_0_imports.py +++ b/obstools/tests/test_0_imports.py @@ -1,9 +1,11 @@ def test_stdb_import(): import stdb + def test_obspy_import(): import obspy + def test_obstools_modules(): import obstools import obstools.atacr @@ -13,5 +15,6 @@ def test_obstools_modules(): import matplotlib matplotlib.use('Agg') from obstools.atacr import plotting - from obstools.atacr.classes import DayNoise, StaNoise, TFNoise, EventStream, Power, Cross, Rotation + from obstools.atacr.classes import DayNoise, StaNoise, TFNoise + from obstools.atacr.classes import EventStream, Power, Cross, Rotation from obstools.comply.classes import Comply diff --git a/obstools/tests/test_1_args.py b/obstools/tests/test_1_args.py index 263a9a2..f36248f 100644 --- a/obstools/tests/test_1_args.py +++ b/obstools/tests/test_1_args.py @@ -23,7 +23,7 @@ def test_get_daylong_arguments(): dbfile, '--keys', 'MM08']) # channels args = atacr.get_daylong_arguments([ - dbfile, '--channels', 'H,P']) + dbfile, '--channels', '12,P']) with pytest.raises(SystemExit): atacr.get_daylong_arguments([ dbfile, '--channels', 'J,P']) @@ -76,7 +76,7 @@ def test_get_event_arguments(): dbfile, '--keys', '7D.MM08']) # channels args = atacr.get_event_arguments([ - dbfile, '--channels', 'H,P']) + dbfile, '--channels', '12,P']) with pytest.raises(SystemExit): atacr.get_event_arguments([ dbfile, '--channels', 'J,P']) @@ -92,7 +92,7 @@ def test_get_event_arguments(): with pytest.raises(SystemExit): atacr.get_event_arguments([ dbfile, '--end', 'abcd']) - # user auth + # user auth args = atacr.get_event_arguments([ dbfile, '-U', 'user:name']) with pytest.raises(SystemExit): diff --git a/obstools/tests/test_2_a_scripts_noH.py b/obstools/tests/test_2_a_scripts_noH.py index 272ed6f..d0de416 100644 --- a/obstools/tests/test_2_a_scripts_noH.py +++ b/obstools/tests/test_2_a_scripts_noH.py @@ -1,6 +1,6 @@ import os import matplotlib as mpl -if os.environ.get('DISPLAY','') == '': +if os.environ.get('DISPLAY', '') == '': mpl.use('Agg') import stdb @@ -33,43 +33,56 @@ def test_01_data_noH(): '--start', '2012-03-08', '--end', '2012-03-10', '--sampling-rate', '1.0', '--channels', 'P']) atacr.main(args=args0) + + def test_02_daily_noH(): from obstools.scripts import atacr_daily_spectra as atacr args0 = atacr.get_dailyspec_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--figQC', '--figAverage', '--save-fig']) atacr.main(args=args0) + + def test_03_clean_noH(): from obstools.scripts import atacr_clean_spectra as atacr args0 = atacr.get_cleanspec_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--figCross', '--figCoh']) atacr.main(args=args0) + + def test_04_trans_noH(): from obstools.scripts import atacr_transfer_functions as atacr args0 = atacr.get_transfer_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--figTF']) atacr.main(args=args0) + + def test_05_event_noH(): from obstools.scripts import atacr_download_event as atacr args0 = atacr.get_event_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--start', '2012-03-08', '--end', '2012-03-10', - '--min-mag', '6.3', '--max-mag', '6.7', + '--min-mag', '6.3', '--max-mag', '6.7', '--window', '7200.', '--sampling-rate', '1.0', '--channels', 'P']) atacr.main(args=args0) + + def test_06_correct_noH(): from obstools.scripts import atacr_correct_event as atacr args0 = atacr.get_correct_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--figRaw', '--figClean', '--save-fig', '--save']) atacr.main(args=args0) + + def test_07_comply_noH(): from obstools.scripts import comply_calculate as comply args0 = comply.get_comply_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--fig']) comply.main(args=args0) + def test_08_rmtree(): shutil.rmtree(datadir) shutil.rmtree(avgdir) diff --git a/obstools/tests/test_2_b_scripts_noP.py b/obstools/tests/test_2_b_scripts_noP.py index 86196c8..93b2e12 100644 --- a/obstools/tests/test_2_b_scripts_noP.py +++ b/obstools/tests/test_2_b_scripts_noP.py @@ -26,33 +26,43 @@ def test_11_data_noP(): args0 = atacr.get_daylong_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--start', '2012-03-08', '--end', '2012-03-10', - '--sampling-rate', '1.0', '--channels', 'H']) + '--sampling-rate', '1.0', '--channels', '12']) atacr.main(args=args0) + + def test_12_daily_noP(): from obstools.scripts import atacr_daily_spectra as atacr args0 = atacr.get_dailyspec_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--figQC', '--figAverage', '--save-fig']) atacr.main(args=args0) + + def test_13_clean_noP(): from obstools.scripts import atacr_clean_spectra as atacr args0 = atacr.get_cleanspec_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--figCross', '--figCoh']) atacr.main(args=args0) + + def test_14_trans_noP(): from obstools.scripts import atacr_transfer_functions as atacr args0 = atacr.get_transfer_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--figTF']) atacr.main(args=args0) + + def test_15_event_noP(): from obstools.scripts import atacr_download_event as atacr args0 = atacr.get_event_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--start', '2012-03-08', '--end', '2012-03-10', - '--min-mag', '6.3', '--max-mag', '6.7', - '--sampling-rate', '1.0', '--channels', 'H']) + '--min-mag', '6.3', '--max-mag', '6.7', '--window', '7200.', + '--sampling-rate', '1.0', '--channels', '12']) atacr.main(args=args0) + + def test_16_correct_noP(): from obstools.scripts import atacr_correct_event as atacr args0 = atacr.get_correct_arguments([ @@ -60,6 +70,7 @@ def test_16_correct_noP(): '--figClean', '--save-fig', '--save']) atacr.main(args=args0) + def test_17_rmtree(): shutil.rmtree(datadir) shutil.rmtree(avgdir) diff --git a/obstools/tests/test_2_c_scripts_all.py b/obstools/tests/test_2_c_scripts_all.py index b3cf3de..9015a68 100644 --- a/obstools/tests/test_2_c_scripts_all.py +++ b/obstools/tests/test_2_c_scripts_all.py @@ -28,42 +28,55 @@ def test_21_data_all(): '--start', '2012-03-08', '--end', '2012-03-10', '--sampling-rate', '1.0']) atacr.main(args=args0) -def test_22_daily_all(): + + +def test_22_event_all(): from obstools.scripts import atacr_download_event as atacr args0 = atacr.get_event_arguments([ dbfile, '--keys', '7D.M08A', - '--start', '2012-03-09', '--end', '2012-03-10', - '--min-mag', '6.3', '--max-mag', '6.7', + '--start', '2012-03-08', '--end', '2012-03-10', + '--min-mag', '6.3', '--max-mag', '6.7', '--window', '7200.', '--sampling-rate', '1.0']) atacr.main(args=args0) -def test_23_clean_all(): + + +def test_23_daily_all(): from obstools.scripts import atacr_daily_spectra as atacr args0 = atacr.get_dailyspec_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig']) atacr.main(args=args0) -def test_24_trans_all(): + + +def test_24_clean_all(): from obstools.scripts import atacr_clean_spectra as atacr args0 = atacr.get_cleanspec_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--figCross', '--figCoh']) atacr.main(args=args0) -def test_25_event_all(): + + +def test_25_trans_all(): from obstools.scripts import atacr_transfer_functions as atacr args0 = atacr.get_transfer_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--figTF']) atacr.main(args=args0) + + def test_26_correct_all(): from obstools.scripts import atacr_correct_event as atacr args0 = atacr.get_correct_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--figRaw', '--figClean', '--save-fig', '--save']) atacr.main(args=args0) + + def test_27_comply_all(): from obstools.scripts import comply_calculate as comply args0 = comply.get_comply_arguments([ dbfile, '--keys', '7D.M08A', '-O', '--save-fig', '--fig']) comply.main(args=args0) + def test_17_rmtree(): shutil.rmtree(avgdir) shutil.rmtree(specdir) diff --git a/obstools/tests/test_3_utils.py b/obstools/tests/test_3_utils.py index 9c726a8..3638bc4 100644 --- a/obstools/tests/test_3_utils.py +++ b/obstools/tests/test_3_utils.py @@ -7,7 +7,8 @@ import os # Path where data are located -exmpl_path = Path(resource_filename('obstools','examples')) +exmpl_path = Path(resource_filename('obstools', 'examples')) + def test_get_data(): datapath = Path('DATA') / '7D.M08A' @@ -20,6 +21,7 @@ def test_get_data(): assert trNP is not None return trN1, trN2, trNZ, trNP + def test_get_H_data(tmp_path): datapath = Path('DATA') / '7D.M08A' @@ -38,6 +40,7 @@ def test_get_H_data(tmp_path): assert len(trNZ) > 0 assert [len(tr.data) == 0 for tr in trNP] + def test_get_P_data(tmp_path): datapath = Path('DATA') / '7D.M08A' @@ -54,23 +57,25 @@ def test_get_P_data(tmp_path): assert len(trNZ) > 0 assert len(trNP) > 0 + def test_get_P_data_sr1(tmp_path): datapath = Path('DATA') / '7D.M08A' - for filename in glob.glob(os.path.join(datapath, '*H.SAC')): + for filename in glob.glob(os.path.join(datapath, '*.SAC')): shutil.copy(filename, tmp_path) for filename in glob.glob(os.path.join(tmp_path, '*H.SAC')): stP = read(filename) - stP[0].resample(1.) + stP[0].resample(0.5) stP[0].write(filename, format='SAC') tstart = UTCDateTime('2012-03-01') tend = UTCDateTime('2012-03-10') trN1, trN2, trNZ, trNP = utils.get_data(tmp_path, tstart, tend) + def test_get_P_data_sr2(tmp_path): datapath = Path('DATA') / '7D.M08A' - for filename in glob.glob(os.path.join(datapath, '*H.SAC')): + for filename in glob.glob(os.path.join(datapath, '*.SAC')): shutil.copy(filename, tmp_path) for filename in glob.glob(os.path.join(tmp_path, '*H.SAC')): stP = read(filename) @@ -80,6 +85,7 @@ def test_get_P_data_sr2(tmp_path): tend = UTCDateTime('2012-03-10') trN1, trN2, trNZ, trNP = utils.get_data(tmp_path, tstart, tend) + def test_get_event(): datapath = Path('EVENTS') / '7D.M08A' tstart = UTCDateTime('2012-03-08') @@ -91,6 +97,7 @@ def test_get_event(): assert len(trP) > 0 return tr1, tr2, trZ, trP + def test_get_H_event(tmp_path): datapath = Path('EVENTS') / '7D.M08A' @@ -109,6 +116,7 @@ def test_get_H_event(tmp_path): assert len(trZ) > 0 assert [len(tr.data) == 0 for tr in trP] + def test_get_P_event(tmp_path): datapath = Path('EVENTS') / '7D.M08A' @@ -129,26 +137,26 @@ def test_get_P_event(tmp_path): def test_get_P_event_sr1(tmp_path): datapath = Path('EVENTS') / '7D.M08A' - for filename in glob.glob(os.path.join(datapath, '*H.SAC')): + for filename in glob.glob(os.path.join(datapath, '*.SAC')): shutil.copy(filename, tmp_path) for filename in glob.glob(os.path.join(tmp_path, '*H.SAC')): stP = read(filename) - stP[0].resample(5.) + stP[0].resample(0.5) stP[0].write(filename, format='SAC') tstart = UTCDateTime('2012-03-08') tend = UTCDateTime('2012-03-10') tr1, tr2, trZ, trP = utils.get_event(tmp_path, tstart, tend) + def test_get_P_event_sr2(tmp_path): datapath = Path('EVENTS') / '7D.M08A' - for filename in glob.glob(os.path.join(datapath, '*H.SAC')): + for filename in glob.glob(os.path.join(datapath, '*.SAC')): shutil.copy(filename, tmp_path) for filename in glob.glob(os.path.join(tmp_path, '*H.SAC')): stP = read(filename) - stP[0].resample(5.) + stP[0].resample(10.) stP[0].write(filename, format='SAC') tstart = UTCDateTime('2012-03-08') tend = UTCDateTime('2012-03-10') tr1, tr2, trZ, trP = utils.get_event(tmp_path, tstart, tend) - diff --git a/obstools/tests/test_4_classes.py b/obstools/tests/test_4_classes.py index 4cbec07..7b95685 100644 --- a/obstools/tests/test_4_classes.py +++ b/obstools/tests/test_4_classes.py @@ -3,14 +3,17 @@ from . import get_meta import pytest + def test_daynoise_demo(): return DayNoise('demo') + def test_day_QC(tmp_path): daynoise = test_daynoise_demo() daynoise.QC_daily_spectra(fig_QC=True, save=tmp_path) return daynoise + def test_day_ncomp_opts(tmp_path): daynoise = test_daynoise_demo() daynoise.ncomp = 3 @@ -18,30 +21,37 @@ def test_day_ncomp_opts(tmp_path): daynoise.ncomp = 2 daynoise.QC_daily_spectra(fig_QC=True, save=tmp_path, smooth=False) + def test_day_average(tmp_path): daynoise = test_daynoise_demo() - daynoise.average_daily_spectra(fig_average=True, fig_coh_ph=True, + daynoise.average_daily_spectra( + fig_average=True, fig_coh_ph=True, save=tmp_path) return daynoise + def test_day_save(tmp_path): daynoise = test_daynoise_demo() d = tmp_path / "tmp" daynoise.save(d) + def test_stanoise_demo(): return StaNoise('demo') + def test_stanoise_day_demo(): daynoise = test_daynoise_demo() stanoise = StaNoise(daylist=[daynoise]) return stanoise + def test_sta_QC(tmp_path): stanoise = test_stanoise_demo() stanoise.QC_sta_spectra(fig_QC=True, save=tmp_path) return stanoise + def test_sta_ncomp_opts(tmp_path): stanoise = test_stanoise_demo() stanoise.ncomp = 3 @@ -51,6 +61,7 @@ def test_sta_ncomp_opts(tmp_path): stanoise.initialized = None stanoise.QC_sta_spectra(fig_QC=True, save=tmp_path) + def test_sta_average(tmp_path): stanoise = test_stanoise_demo() with pytest.raises(Exception): @@ -59,6 +70,7 @@ def test_sta_average(tmp_path): stanoise.average_sta_spectra(fig_average=True, save=tmp_path) return stanoise + def test_sta_operations(): dn1 = test_daynoise_demo() dn2 = test_daynoise_demo() @@ -77,6 +89,7 @@ def test_sta_operations(): with pytest.raises(Exception): assert sn.init() + def test_sta_save(tmp_path): stanoise = test_stanoise_demo() d = tmp_path / "tmp" @@ -86,12 +99,14 @@ def test_sta_save(tmp_path): stanoise.average_sta_spectra() stanoise.save(d) + def test_evstream_demo(tmp_path): eventstream = EventStream('demo') d = tmp_path / "tmp" eventstream.save(d) return eventstream + def test_tfnoise_day_demo(tmp_path): daynoise = test_daynoise_demo() with pytest.raises(TypeError): @@ -105,26 +120,30 @@ def test_tfnoise_day_demo(tmp_path): tfnoise_day.save(d) return tfnoise_day + def test_tfnoise_sta_demo(tmp_path): stanoise = test_sta_average(tmp_path) tfnoise_sta = TFNoise(stanoise) tfnoise_sta.transfer_func() return tfnoise_sta + def test_comply_day_demo(tmp_path): daynoise = test_day_average(tmp_path) sta = get_meta.get_stdb() - comply_day = Comply(objnoise=daynoise, sta=sta) + comply_day = Comply(objnoise=daynoise, elev=sta.elevation*1.e3) comply_day.calculate_compliance() return comply_day + def test_comply_sta_demo(tmp_path): stanoise = test_sta_average(tmp_path) sta = get_meta.get_stdb() - comply_sta = Comply(objnoise=stanoise, sta=sta) + comply_sta = Comply(objnoise=stanoise, elev=sta.elevation*1.e3) comply_sta.calculate_compliance() return comply_sta + def test_comply_fail(tmp_path): import pytest import pickle @@ -135,15 +154,15 @@ def test_comply_fail(tmp_path): sta = get_meta.get_stdb() with pytest.raises(Exception): - assert Comply(sta=sta, objnoise=[]) + assert Comply(elev=sta.elevation*1.e3, objnoise=[]) objnoise = test_day_average(tmp_path) objnoise.av = None with pytest.raises(Exception): - assert Comply(sta=sta, objnoise=objnoise) + assert Comply(elev=sta.elevation*1.e3, objnoise=objnoise) daynoise = test_day_average(tmp_path) - comply_day = Comply(objnoise=daynoise, sta=sta) + comply_day = Comply(objnoise=daynoise, elev=sta.elevation*1.e3) comply_day.calculate_compliance() d = tmp_path / "tmp" comply_day.save(d, form='pkl') diff --git a/setup.py b/setup.py index 0beb93a..4d79f9c 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ def find_version(*paths): fname = os.path.join(os.path.dirname(__file__), *paths) - with open(fname) as fp: + with open(fname, encoding='utf-8') as fp: code = fp.read() match = re.search(r"^__version__ = ['\"]([^'\"]*)['\"]", code, re.M) if match: @@ -32,13 +32,14 @@ def find_version(*paths): 'License :: OSI Approved :: MIT License', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', - 'Programming Language :: Python :: 3.8'], + 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9'], install_requires=['numpy', 'obspy', 'stdb', 'pandas'], python_requires='>=3.6', packages=setuptools.find_packages(), include_package_data=True, package_data={'': ['examples/meta/*.pkl', 'examples/event/*.pkl', - 'examples/data/*.SAC']}, + 'examples/data/*.SAC', 'examples/event/*.SAC']}, # scripts=scripts) entry_points={ 'console_scripts':