marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

Re-Added local changes that had been lost due to technical problems ( no access to old machine )

Browse Source
This commit is contained in:
jeldrik
2022-10-04 11:44:31 +02:00
parent 2d33a60421
commit 0c3fca9299
7 changed files with 110044 additions and 109102 deletions
Download Patch File Download Diff File Expand all files Collapse all files
pylot/core/pick/charfuns.py
+25 -4
View File
@@ -16,10 +16,9 @@ autoregressive prediction: application ot local and regional distances, Geophys.
:author: MAGS2 EP3 working group
"""
import numpy as np
from obspy.core import Stream
from scipy import signal
from obspy.core import Stream
class CharacteristicFunction(object):
@@ -259,7 +258,7 @@ class HOScf(CharacteristicFunction):
"""
Function to calculate skewness (statistics of order 3) or kurtosis
(statistics of order 4), using one long moving window, as published
in Kueperkoch et al. (2010).
in Kueperkoch et al. (2010), or order 2, i.e. STA/LTA.
:param data: data, time series (whether seismogram or CF)
:type data: tuple
:return: HOS cf
@@ -276,28 +275,47 @@ class HOScf(CharacteristicFunction):
elif self.getOrder() == 4: # this is kurtosis
y = np.power(xnp, 4)
y1 = np.power(xnp, 2)
elif self.getOrder() == 2: # this is variance, used for STA/LTA processing
y = np.power(xnp, 2)
y1 = np.power(xnp, 2)
# Initialisation
# t2: long term moving window
ilta = int(round(self.getTime2() / self.getIncrement()))
ista = int(round((self.getTime2() / 10) / self.getIncrement())) # TODO: still hard coded!!
lta = y[0]
lta1 = y1[0]
sta = y[0]
# moving windows
LTA = np.zeros(len(xnp))
STA = np.zeros(len(xnp))
for j in range(0, len(xnp)):
if j < 4:
LTA[j] = 0
STA[j] = 0
elif j <= ista:
lta = (y[j] + lta * (j - 1)) / j
if self.getOrder() == 2:
sta = (y[j] + sta * (j - 1)) / j
# elif j < 4:
elif j <= ilta:
lta = (y[j] + lta * (j - 1)) / j
lta1 = (y1[j] + lta1 * (j - 1)) / j
if self.getOrder() == 2:
sta = (y[j] - y[j - ista]) / ista + sta
else:
lta = (y[j] - y[j - ilta]) / ilta + lta
lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1
if self.getOrder() == 2:
sta = (y[j] - y[j - ista]) / ista + sta
# define LTA
if self.getOrder() == 3:
LTA[j] = lta / np.power(lta1, 1.5)
elif self.getOrder() == 4:
LTA[j] = lta / np.power(lta1, 2)
else:
LTA[j] = lta
STA[j] = sta
# remove NaN's with first not-NaN-value,
# so autopicker doesnt pick discontinuity at start of the trace
@@ -306,7 +324,10 @@ class HOScf(CharacteristicFunction):
first = ind[0]
LTA[:first] = LTA[first]
self.cf = LTA
if self.getOrder() > 2:
self.cf = LTA
else: # order 2 means STA/LTA!
self.cf = STA / LTA
self.xcf = x