# Copyright 2019 Pascal Audet
#
# This file is part of OrientPy.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# -*- coding: utf-8 -*-
import os
import numpy as np
import matplotlib.pyplot as plt
from obspy import Trace, Stream, UTCDateTime
from scipy.stats import circmean as cmean
from scipy.stats import circstd as cstd
from scipy.stats import hmean as hm
from orientpy import io, utils, plotting
from obspy.signal.rotate import rotate_rt_ne, rotate_ne_rt
from pkg_resources import resource_filename
[docs]
class Orient(object):
"""
A Orient object contains Class attributes that associate
station information with a single event (i.e., earthquake)
metadata.
Note
----
The object is initialized with the ``sta`` field only, and
other attributes are added to the object as the analysis proceeds.
Parameters
----------
sta : :class:`~stdb.StDbElement`
Object containing station information - from :mod:`~stdb` database.
meta : :class:`~orientpy.classes.Meta`
Object of metadata information for single event (initially set to None)
data : :class:`~obspy.core.Stream`
Stream object containing the three-component seismograms
zcomp : str
Vertical Component Identifier. Should be a single character.
This is different then 'Z' only for fully unknown component
orientation (i.e., components are 1, 2, 3)
"""
def __init__(self, sta, zcomp='Z'):
# Attributes from parameters
self.sta = sta
# Initialize meta and data objects as None
self.meta = None
self.data = None
self.zcomp = zcomp
[docs]
def add_event(self, event, gacmin=None, gacmax=None,
depmax=1000., returned=False):
"""
Adds event metadata to Orient object.
Parameters
----------
event : :class:`~obspy.core.event`
Event XML object
gacmin : float
Minimum great-arc circle distance to consider (deg)
gacmax : float
Maximum great-arc circle distance to consider (deg)
depmax : float
Maximum event depth to consider (km)
returned : bool
Whether or not to return the ``accept`` attribute
Returns
-------
accept : bool
Whether or not the object is accepted for further analysis
Attributes
----------
meta : :class:`~orientpy.classes.Meta`
Meta data object
"""
from obspy.core.event.event import Event
if not isinstance(event, Event):
raise Exception("Event has incorrect type")
# Store as object attributes
self.meta = Meta(sta=self.sta, event=event,
gacmin=gacmin, gacmax=gacmax,
depmax=depmax)
if returned:
return self.meta.accept
[docs]
def download_data(self, client, new_sr=None, t1=None, t2=None,
returned=False, verbose=False):
"""
Downloads seismograms based on a time interval.
Parameters
----------
client : :class:`~obspy.client.fdsn.Client`
Client object
new_sr : float
New sampling rate (Hz)
t1 : float
Start time to consider (sec). Can be float or
:class:`~obspy.core.UTCDateTime` object.
t2 : float
End time to consider (sec). Can be float or
:class:`~obspy.core.UTCDateTime` object.
returned : bool
Whether or not to return the ``accept`` attribute
verbose : bool
Whether or not to print messages to screen during run-time
Returns
-------
accept : bool
Whether or not the object is accepted for further analysis
Attributes
----------
data : :class:`~obspy.core.Stream`
Stream containing :class:`~obspy.core.Trace` objects
"""
if self.meta is None:
raise Exception("Requires event data as attribute - aborting")
if not self.meta.accept:
return
try:
# Define start time for request
tstart = self.meta.time + t1
tend = self.meta.time + t2
except Exception:
raise Exception(" Start and end request times have not been set")
self.meta.accept = False
if returned:
return self.meta.accept
# Get waveforms
if verbose:
print("* Requesting Waveforms: ")
print("* Startime: " + tstart.strftime("%Y-%m-%d %H:%M:%S"))
print("* Endtime: " + tend.strftime("%Y-%m-%d %H:%M:%S"))
# Download data
err, stream = io.download_data(
client=client, sta=self.sta, start=tstart, end=tend,
zcomp=self.zcomp, verbose=verbose)
if err or stream is None:
print("* Waveform Request Failed...Skipping")
self.meta.accept = False
if returned:
return self.meta.accept
# Store as attributes with traces in dictionary
try:
trE = stream.select(component='E')[0]
trN = stream.select(component='N')[0]
trZ = stream.select(component=self.zcomp)[0]
self.data = Stream(traces=[trZ, trN, trE])
# Filter Traces and resample
if new_sr is not None:
self.data.filter('lowpass', freq=0.5*new_sr,
corners=2, zerophase=True)
self.data.resample(new_sr)
# If there is no ZNE, perhaps there is Z12?
except Exception as e:
tr1 = stream.select(component='1')[0]
tr2 = stream.select(component='2')[0]
trZ = stream.select(component=self.zcomp)[0]
self.data = Stream(traces=[trZ, tr1, tr2])
# Filter Traces and resample
if new_sr is not None:
self.data.filter('lowpass', freq=0.5*new_sr,
corners=2, zerophase=True)
self.data.resample(new_sr)
if not np.sum([np.std(tr.data) for tr in self.data]) > 0.:
print('Error: Data are all zeros')
self.meta.accept = False
if returned:
return self.meta.accept
def save(self, file):
import pickle
output = open(file, 'wb')
pickle.dump(self, output)
output.close()
return
[docs]
class BNG(Orient):
"""
A BNG object inherits from :class:`~orientpy.classes.Orient`.
This object contains a method to estimate station orientation
based on P-wave particle motion.
Notes
-----
This class is designed after the method developed by Braunmiller,
Nabelek and Ghods (2020) [1]_. It is slightly more flexible, however,
in that it can handle either regional or teleseismic P-wave data
and provides additional quality-control mesures (e.g., SNR, 1-R/Z).
References
----------
.. [1] Braunmiller, J., Nabelek, J., and Ghods, A., 2020, Sensor orientation
of Iranian broadband seismic stations from P-wave particle motion,
*Seismological Research Letters*, doi:10.1785/0220200019.
Parameters
----------
sta : :class:`~stdb.StDbElement`
Object containing station information - from :mod:`~stdb` database.
meta : :class:`~orientpy.classes.Meta`
Object of metadata information for single event (initially set to None)
data : :class:`~obspy.core.Stream`
Stream object containing the three-component seismograms
"""
def __init__(self, sta, zcomp='Z'):
Orient.__init__(self, sta, zcomp=zcomp)
[docs]
def calc(self, dphi, dts, tt, bp=None, showplot=False):
"""
Method to estimate azimuth of component `?H1` (or `?HN`). This
method minimizes the energy (RMS) of the transverse component of
P-wave data within some bandwidth.
Parameters
----------
dphi : float
Azimuth increment for search (deg)
dts : float
Length of time window on either side of
predicted P-wave arrival time (sec)
tt : list
List of two floats containing the time picks relative
to P-wave time, within which to perform the estimation of
station orientation (sec)
bp : list
List of two floats containing the low- and high-frequency
corners of a bandpass filter (Hz)
showplot : bool
Whether or not to plot waveforms.
Attributes
----------
meta.phi : float
Azimuth of H1 (or HN) component (deg)
meta.cc : float
Cross-correlation coefficient between
vertical and radial component
meta.snr : float
Signal-to-noise ratio of P-wave measured on the
vertical seismogram
meta.TR : float
Measure of the transverse to radial ratio. In reality
this is 1 - T/R
meta.RZ : float
Measure of the radial to vertical ratio. In reality
this is 1 - R/Z
"""
# Work on a copy of the waveform data
stream = self.data.copy()
# Identify components in Stream
comps_id = [tr.stats.channel[2] for tr in stream]
test_sets = {
'ZNE': {'Z', 'N', 'E'},
'Z12': {'Z', '1', '2'},
'Z23': {'Z', '2', '3'},
'312': {'3', '1', '2'},
'123': {'1', '2', '3'}
}
# Probably should raise an exception if set is '123',
# as no correction is estimated for the vertical component
for test_key in test_sets:
test_set = test_sets[test_key]
if test_set.issubset(set(comps_id)):
# Use sets to avoid sorting complications
comps_codes = list(test_key)
break
# Temporarily modify channel codes,
# assuming that N/E are not oriented properly
channel_code_prefix = stream[0].stats.channel[:2]
stream.select(component=comps_codes[1])[0].stats.channel = channel_code_prefix + '1'
stream.select(component=comps_codes[2])[0].stats.channel = channel_code_prefix + '2'
stream.select(component=comps_codes[0])[0].stats.channel = channel_code_prefix + self.zcomp.upper()
# Filter if specified
if bp is not None:
stream.filter('bandpass', freqmin=bp[0],
freqmax=bp[1], zerophase=True)
# Get data and noise based on symmetric waveform wrt arrival
start = stream[0].stats.starttime
stnoise = stream.copy().trim(start, start+dts+tt[0])
stdata = stream.copy().trim(start+dts+tt[0], start+dts+tt[1])
# Define signal and noise
tr1 = stdata.select(component='1')[0].copy()
tr2 = stdata.select(component='2')[0].copy()
trZ = stdata.select(component=self.zcomp.upper())[0].copy()
ntrZ = stnoise.select(component=self.zcomp.upper())[0].copy()
# Calculate and store SNR as attribute
self.meta.snr = 10.*np.log10(
utils.rms(trZ)*utils.rms(trZ)/utils.rms(ntrZ)/utils.rms(ntrZ))
# Search through azimuths from 0 to 180 deg and find best-fit azimuth
ang = np.arange(0., 180., dphi)
cc1 = np.zeros(len(ang))
cc2 = np.zeros(len(ang))
cc3 = np.zeros(len(ang))
cc4 = np.zeros(len(ang))
for k, a in enumerate(ang):
R, T = rotate_ne_rt(tr1.data, tr2.data, a)
covmat = np.corrcoef(R, trZ.data)
cc1[k] = covmat[0, 1]
cc2[k] = 1. - utils.rms(T)/utils.rms(R)
cc3[k] = utils.rms(T)
cc4[k] = 1. - utils.rms(R)/utils.rms(trZ.data)
# Get argument of minimum of cc3 and store useful measures
ia = cc3.argmin()
self.meta.cc = cc1[ia]
self.meta.TR = cc2[ia]
self.meta.RZ = cc4[ia]
# correct for angles above 360
phi = (self.meta.baz - float(ia)*dphi)
# Use azimuth where CC is negative
if self.meta.cc < 0.:
phi += 180.
if phi < 0.:
phi += 360.
if phi >= 360.:
phi -= 360.
# Store the best-fit azimuth
self.meta.phi = phi
# If a plot is requested, rotate Z12 to ZNE and then ZRT
if showplot:
# Now rotate components to proper N, E and then R, T
sttmp = stream.copy()
# Apply filter if defined previously
if bp:
sttmp.filter('bandpass', freqmin=bp[0],
freqmax=bp[1], zerophase=True)
# Copy traces
trN = sttmp.select(component='1')[0].copy()
trE = sttmp.select(component='2')[0].copy()
# Rotating from 1,2 to N,E is the negative of
# rotation from RT to NE, with baz corresponding
# to azim of component 1, or phi previously determined
azim = self.meta.phi
N, E = rotate_rt_ne(trN.data, trE.data, azim)
trN.data = -1.*N
trE.data = -1.*E
# Update stats of streams
trN.stats.channel = trN.stats.channel[:-1] + 'N'
trE.stats.channel = trE.stats.channel[:-1] + 'E'
# Store corrected traces in new stream and rotate to
# R, T using back-azimuth
stcorr = Stream(traces=[trN, trE])
stcorr.rotate('NE->RT', back_azimuth=self.meta.baz)
# Merge original and corrected streams
st = stream + stcorr
# Plot
plot = plotting.plot_bng_waveforms(self, st, dts, tt)
plot.show()
return
[docs]
class DL(Orient):
"""
A DL object inherits from :class:`~orientpy.classes.Orient`. This object
contains a method to estimate station orientation based on Rayleigh-wave
particle motion.
Notes
-----
This algorithm is heavily based on the code `DLOPy` by
Doran and Laske (2017) [2]_. The original code can be found here:
https://igppweb.ucsd.edu/~adoran/DLOPy.html
There are also a couple of versions available on GitHub:
- https://github.com/kschramm-usgs/DLOPy
- https://github.com/jbrussell/DLOPy_v1.0
References
----------
.. [2] Doran, A. K., and G. Laske (2017). Ocean-bottom seismometer
instrument orientations via automated Rayleigh-wave arrival-angle
measurements, *Bulletin of the Seismological Society of America*,
107(2), doi:10.1785/0120160165.
Parameters
----------
sta : :class:`~stdb.StDbElement`
Object containing station information - from :mod:`~stdb` database.
meta : :class:`~orientpy.classes.Meta`
Object of metadata information for single event (initially set to None)
data : :class:`~obspy.core.Stream`
Stream object containing the three-component seismograms
"""
def __init__(self, sta, zcomp='Z'):
Orient.__init__(self, sta, zcomp=zcomp)
[docs]
def calc(self, showplot=False):
"""
Method to estimate azimuth of component `?H1` (or `?HN`). This
method maximizes the normalized covariance between the
Hilbert-transformed vertical component and the radial component of
Rayleigh-wave data measured at multiple periods. This is done for the
shortest (R1) and longest (R2) orbits of fundamental-mode Rayleigh
waves. Window selection is done based on average group velocity
extracted from a global model of Rayleigh-wave dispersion.
Parameters
----------
showplot : bool
Whether or not to plot waveforms.
Attributes
----------
meta.R1phi : list
List of azimuth of H1 (or HN) component (deg) based on R1, measured
at different periods
meta.R1cc : float
Cross-correlation coefficient between hilbert-transformed vertical
and radial component for R1, measured at different periods
meta.R2phi : list
List of azimuth of H1 (or HN) component (deg) based on R2, measured
at different periods
meta.R2cc : float
Cross-correlation coefficient between hilbert-transformed vertical
and radial component for R2, measured at different periods
"""
# Work on a copy of the waveform data
stream = self.data.copy()
# Initialize surface wave arrays
nper = 7
R1phi = np.zeros(nper)
R1cc = np.zeros(nper)
R2phi = np.zeros(nper)
R2cc = np.zeros(nper)
# Load group velocity maps
map10 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.10.txt'))
map15 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.15.txt'))
map20 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.20.txt'))
map25 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.25.txt'))
map30 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.30.txt'))
map35 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.35.txt'))
map40 = np.loadtxt(resource_filename('orientpy',
'dispmaps/R.gv.40.txt'))
# Get parameters for R2
Rearth = 6371.25
circE = 2.*np.pi*Rearth
dist2 = circE - self.meta.epi_dist
baz2 = self.meta.baz + 180.
if baz2 >= 360.:
baz2 -= 360.
# Check data length, data quality
if utils.checklen(self.data, 4.*60.*60.):
raise Exception(" Error: Length Incorrect")
# Get path-averaged group velocities
Ray1, Ray2 = utils.pathvels(
self.sta.latitude, self.sta.longitude,
self.meta.lat, self.meta.lon,
map10, map15, map20, map25, map30, map35, map40)
# Calculate arrival angle for each frequency and orbit
Rf = [40., 35., 30., 25., 20., 15., 10.]
LPF = [0.035, 0.03, 0.025, 0.02, 0.015, 0.01, 0.005]
HPF = [0.045, 0.04, 0.035, 0.03, 0.025, 0.02, 0.015]
winlen1 = [20., 17., 14., 12., 10., 10., 7.]
winlen2 = [24., 20., 16., 13., 10., 10., 7.]
flist = list(zip(Rf, HPF, LPF, winlen1, winlen2))
for k, item in enumerate(flist):
# R1 path
R1phi[k], R1cc[k] = utils.DLcalc(
stream, item[0], item[1],
item[2], self.meta.epi_dist, self.meta.baz, Ray1,
winlen=item[3], ptype=0, zcomp=self.zcomp)
# R2 path
R2phi[k], R2cc[k] = utils.DLcalc(
stream, item[0], item[1],
item[2], dist2, baz2, Ray2,
winlen=item[4], ptype=0, zcomp=self.zcomp)
# Store azimuths and CC values as attributes
self.meta.R1phi = R1phi
self.meta.R2phi = R2phi
self.meta.R1cc = R1cc
self.meta.R2cc = R2cc
return