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

just cleaning up the code to meet coding conventions

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2015-06-12 09:02:00 +02:00
parent c5da8fd994
commit 883fdf6bf5
2 changed files with 643 additions and 518 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

455
pylot/core/pick/Picker.py
View File

@ -3,58 +3,65 @@
Created Dec 2014 to Feb 2015 Created Dec 2014 to Feb 2015
Implementation of the automated picking algorithms published and described in: Implementation of the automated picking algorithms published and described in:
Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group, 2010: Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group,
Automated determination of P-phase arrival times at regional and local distances 2010: Automated determination of P-phase arrival times at regional and local
using higher order statistics, Geophys. J. Int., 181, 1159-1170 distances using higher order statistics, Geophys. J. Int., 181, 1159-1170
Kueperkoch, L., Meier, T., Bruestle, A., Lee, J., Friederich, W., & Egelados Kueperkoch, L., Meier, T., Bruestle, A., Lee, J., Friederich, W., & Egelados
Working Group, 2012: Automated determination of S-phase arrival times using Working Group, 2012: Automated determination of S-phase arrival times using
autoregressive prediction: application ot local and regional distances, Geophys. J. Int., autoregressive prediction: application ot local and regional distances,
188, 687-702. Geophys. J. Int., 188, 687-702.
The picks with the above described algorithms are assumed to be the most likely picks. The picks with the above described algorithms are assumed to be the most likely
For each most likely pick the corresponding earliest and latest possible picks are picks. For each most likely pick the corresponding earliest and latest possible
calculated after Diehl & Kissling (2009). picks are calculated after Diehl & Kissling (2009).
:author: MAGS2 EP3 working group / Ludger Kueperkoch :author: MAGS2 EP3 working group / Ludger Kueperkoch
""" """
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from pylot.core.pick.utils import * from pylot.core.pick.utils import *
from pylot.core.pick.CharFuns import CharacteristicFunction from pylot.core.pick.CharFuns import CharacteristicFunction
class AutoPicking(object): class AutoPicking(object):
''' '''
Superclass of different, automated picking algorithms applied on a CF determined Superclass of different, automated picking algorithms applied on a CF
using AIC, HOS, or AR prediction. determined using AIC, HOS, or AR prediction.
''' '''
def __init__(self, cf, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None, Pick1=None):
def __init__(self, cf, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None,
Pick1=None):
''' '''
:param: cf, characteristic function, on which the picking algorithm is applied :param cf: characteristic function, on which the picking algorithm is
:type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object applied
:type cf: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
:param: TSNR, length of time windows around pick used to determine SNR [s] :param TSNR: length of time windows for SNR determination - [s]
:type: tuple (T_noise, T_gap, T_signal) :type TSNR: tuple (T_noise, T_gap, T_signal)
:param: PickWindow, length of pick window [s] :param PickWindow: length of pick window - [s]
:type: float :type PickWindow: float
:param: iplot, no. of figure window for plotting interims results :param iplot: no. of figure window for plotting interims results
:type: integer :type iplot: integer
:param: aus ("artificial uplift of samples"), find local minimum at i if aic(i-1)*(1+aus) >= aic(i) :param aus: aus ("artificial uplift of samples"), find local minimum at
:type: float i if aic(i-1)*(1+aus) >= aic(i)
:type aus: float
:param: Tsmooth, length of moving smoothing window to calculate smoothed CF [s] :param Tsmooth: length of moving window to calculate smoothed CF - [s]
:type: float :type Tsmooth: float
:param: Pick1, initial (prelimenary) onset time, starting point for PragPicker and :param Pick1: initial (prelimenary) onset time, starting point for
EarlLatePicker PragPicker
:type: float :type Pick1: float
''' '''
assert isinstance(cf, CharacteristicFunction), "%s is not a CharacteristicFunction object" % str(cf) assert isinstance(cf,
CharacteristicFunction), "%s is of wrong type" % str(
cf)
self.cf = cf.getCF() self.cf = cf.getCF()
self.Tcf = cf.getTimeArray() self.Tcf = cf.getTimeArray()
@ -80,9 +87,8 @@ class AutoPicking(object):
PickWindow=self.getPickWindow(), PickWindow=self.getPickWindow(),
aus=self.getaus(), aus=self.getaus(),
Tsmooth=self.getTsmooth(), Tsmooth=self.getTsmooth(),
Pick1=self.getpick1()) Pick1=self.getpick1())
def getTSNR(self): def getTSNR(self):
return self.TSNR return self.TSNR
@ -112,7 +118,7 @@ class AutoPicking(object):
def getSNR(self): def getSNR(self):
return self.SNR return self.SNR
def getSlope(self): def getSlope(self):
return self.slope return self.slope
@ -136,144 +142,156 @@ class AICPicker(AutoPicking):
''' '''
Method to derive the onset time of an arriving phase based on CF Method to derive the onset time of an arriving phase based on CF
derived from AIC. In order to get an impression of the quality of this inital pick, derived from AIC. In order to get an impression of the quality of this inital pick,
a quality assessment is applied based on SNR and slope determination derived from the CF, a quality assessment is applied based on SNR and slope determination derived from the CF,
from which the AIC has been calculated. from which the AIC has been calculated.
''' '''
def calcPick(self): def calcPick(self):
print 'AICPicker: Get initial onset time (pick) from AIC-CF ...' print 'AICPicker: Get initial onset time (pick) from AIC-CF ...'
self.Pick = None self.Pick = None
self.slope = None self.slope = None
self.SNR = None self.SNR = None
#find NaN's # find NaN's
nn = np.isnan(self.cf) nn = np.isnan(self.cf)
if len(nn) > 1: if len(nn) > 1:
self.cf[nn] = 0 self.cf[nn] = 0
#taper AIC-CF to get rid off side maxima # taper AIC-CF to get rid off side maxima
tap = np.hanning(len(self.cf)) tap = np.hanning(len(self.cf))
aic = tap * self.cf + max(abs(self.cf)) aic = tap * self.cf + max(abs(self.cf))
#smooth AIC-CF # smooth AIC-CF
ismooth = int(round(self.Tsmooth / self.dt)) ismooth = int(round(self.Tsmooth / self.dt))
aicsmooth = np.zeros(len(aic)) aicsmooth = np.zeros(len(aic))
if len(aic) < ismooth: if len(aic) < ismooth:
print 'AICPicker: Tsmooth larger than CF!' print 'AICPicker: Tsmooth larger than CF!'
return return
else: else:
for i in range(1, len(aic)): for i in range(1, len(aic)):
if i > ismooth: if i > ismooth:
ii1 = i - ismooth ii1 = i - ismooth
aicsmooth[i] = aicsmooth[i - 1] + (aic[i] - aic[ii1]) / ismooth aicsmooth[i] = aicsmooth[i - 1] + (aic[i] - aic[
else: ii1]) / ismooth
aicsmooth[i] = np.mean(aic[1 : i]) else:
#remove offset aicsmooth[i] = np.mean(aic[1: i])
# remove offset
offset = abs(min(aic) - min(aicsmooth)) offset = abs(min(aic) - min(aicsmooth))
aicsmooth = aicsmooth - offset aicsmooth = aicsmooth - offset
#get maximum of 1st derivative of AIC-CF (more stable!) as starting point # get maximum of 1st derivative of AIC-CF (more stable!) as starting
# point
diffcf = np.diff(aicsmooth) diffcf = np.diff(aicsmooth)
#find NaN's # find NaN's
nn = np.isnan(diffcf) nn = np.isnan(diffcf)
if len(nn) > 1: if len(nn) > 1:
diffcf[nn] = 0 diffcf[nn] = 0
#taper CF to get rid off side maxima # taper CF to get rid off side maxima
tap = np.hanning(len(diffcf)) tap = np.hanning(len(diffcf))
diffcf = tap * diffcf * max(abs(aicsmooth)) diffcf = tap * diffcf * max(abs(aicsmooth))
icfmax = np.argmax(diffcf) icfmax = np.argmax(diffcf)
#find minimum in AIC-CF front of maximum
lpickwindow = int(round(self.PickWindow / self.dt))
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1):
if aicsmooth[i - 1] >= aicsmooth[i]:
self.Pick = self.Tcf[i]
break
#if no minimum could be found:
#search in 1st derivative of AIC-CF
if self.Pick is None:
for i in range(icfmax -1, max([icfmax -lpickwindow, 2]), -1):
if diffcf[i -1] >= diffcf[i]:
self.Pick = self.Tcf[i]
break
#quality assessment using SNR and slope from CF
if self.Pick is not None:
#get noise window
inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1])
#check, if these are counts or m/s, important for slope estimation!
#this is quick and dirty, better solution?
if max(self.Data[0].data < 1e-3):
self.Data[0].data = self.Data[0].data * 1000000
#get signal window
isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2])
#calculate SNR from CF
self.SNR = max(abs(aic[isignal] - np.mean(aic[isignal]))) / max(abs(aic[inoise] \
- np.mean(aic[inoise])))
#calculate slope from CF after initial pick
#get slope window
tslope = self.TSNR[3] #slope determination window
islope = np.where((self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)])) \
& (self.Tcf >= self.Pick))
#find maximum within slope determination window
#'cause slope should be calculated up to first local minimum only!
imax = np.argmax(self.Data[0].data[islope])
if imax == 0:
print 'AICPicker: Maximum for slope determination right at the beginning of the window!'
print 'Choose longer slope determination window!'
return
islope = islope[0][0 :imax]
dataslope = self.Data[0].data[islope]
#calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope), 1)
P = np.polyfit(xslope, dataslope, 1)
datafit = np.polyval(P, xslope)
if datafit[0] >= datafit[len(datafit) - 1]:
print 'AICPicker: Negative slope, bad onset skipped!'
return
self.slope = 1 / tslope * datafit[len(dataslope) - 1] - datafit[0] # find minimum in AIC-CF front of maximum
lpickwindow = int(round(self.PickWindow / self.dt))
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1):
if aicsmooth[i - 1] >= aicsmooth[i]:
self.Pick = self.Tcf[i]
break
# if no minimum could be found:
# search in 1st derivative of AIC-CF
if self.Pick is None:
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1):
if diffcf[i - 1] >= diffcf[i]:
self.Pick = self.Tcf[i]
break
# quality assessment using SNR and slope from CF
if self.Pick is not None:
# get noise window
inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0],
self.TSNR[1])
# check, if these are counts or m/s, important for slope estimation!
# this is quick and dirty, better solution?
if max(self.Data[0].data < 1e-3):
self.Data[0].data *= 1000000
# get signal window
isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2])
# calculate SNR from CF
self.SNR = max(abs(aic[isignal] - np.mean(aic[isignal]))) / \
max(abs(aic[inoise] - np.mean(aic[inoise])))
# calculate slope from CF after initial pick
# get slope window
tslope = self.TSNR[3] # slope determination window
islope = np.where(
(self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)]))
and (self.Tcf >= self.Pick))
# find maximum within slope determination window
# 'cause slope should be calculated up to first local minimum only!
imax = np.argmax(self.Data[0].data[islope])
if imax == 0:
print 'AICPicker: Maximum for slope determination right at ' \
'the beginning of the window!'
print 'Choose longer slope determination window!'
return
islope = islope[0][0:imax]
dataslope = self.Data[0].data[islope]
# calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope), 1)
P = np.polyfit(xslope, dataslope, 1)
datafit = np.polyval(P, xslope)
if datafit[0] >= datafit[len(datafit) - 1]:
print 'AICPicker: Negative slope, bad onset skipped!'
return
self.slope = 1 / tslope * datafit[len(dataslope) - 1] - datafit[0]
else: else:
self.SNR = None self.SNR = None
self.slope = None self.slope = None
if self.iplot > 1: if self.iplot > 1:
p = plt.figure(self.iplot) p = plt.figure(self.iplot)
x = self.Data[0].data x = self.Data[0].data
p1, = plt.plot(self.Tcf, x / max(x), 'k') p1, = plt.plot(self.Tcf, x / max(x), 'k')
p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r') p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r')
if self.Pick is not None: if self.Pick is not None:
p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2) p3, = plt.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b',
plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick']) linewidth=2)
else: plt.legend([p1, p2, p3],
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick'])
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) else:
plt.yticks([]) plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
plt.title(self.Data[0].stats.station) plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.yticks([])
plt.title(self.Data[0].stats.station)
if self.Pick is not None: if self.Pick is not None:
plt.figure(self.iplot + 1) plt.figure(self.iplot + 1)
p11, = plt.plot(self.Tcf, x, 'k') p11, = plt.plot(self.Tcf, x, 'k')
p12, = plt.plot(self.Tcf[inoise], self.Data[0].data[inoise]) p12, = plt.plot(self.Tcf[inoise], self.Data[0].data[inoise])
p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r') p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal],
p14, = plt.plot(self.Tcf[islope], dataslope, 'g--') 'r')
p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2) p14, = plt.plot(self.Tcf[islope], dataslope, 'g--')
plt.legend([p11, p12, p13, p14, p15], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'], \ p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2)
loc='best') plt.legend([p11, p12, p13, p14, p15],
plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \ ['Data', 'Noise Window', 'Signal Window',
self.SNR, self.slope)) 'Slope Window', 'Slope'],
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) loc='best')
plt.ylabel('Counts') plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (
ax = plt.gca() self.Data[0].stats.station,
plt.yticks([]) self.SNR, self.slope))
ax.set_xlim([self.Tcf[inoise[0][0]] - 5, self.Tcf[isignal[0][len(isignal) - 1]] + 5]) plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.ylabel('Counts')
ax = plt.gca()
plt.yticks([])
ax.set_xlim([self.Tcf[inoise[0][0]] - 5,
self.Tcf[isignal[0][len(isignal) - 1]] + 5])
plt.show() plt.show()
raw_input() raw_input()
plt.close(p) plt.close(p)
if self.Pick is None:
print 'AICPicker: Could not find minimum, picking window too short?'
if self.Pick == None:
print 'AICPicker: Could not find minimum, picking window too short?'
class PragPicker(AutoPicking): class PragPicker(AutoPicking):
''' '''
@ -283,90 +301,95 @@ class PragPicker(AutoPicking):
def calcPick(self): def calcPick(self):
if self.getpick1() is not None: if self.getpick1() is not None:
print 'PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...' print 'PragPicker: Get most likely pick from HOS- or AR-CF using ' \
'pragmatic picking algorithm ...'
self.Pick = None self.Pick = None
self.SNR = None self.SNR = None
self.slope = None self.slope = None
#smooth CF # smooth CF
ismooth = int(round(self.Tsmooth / self.dt)) ismooth = int(round(self.Tsmooth / self.dt))
cfsmooth = np.zeros(len(self.cf)) cfsmooth = np.zeros(len(self.cf))
if len(self.cf) < ismooth: if len(self.cf) < ismooth:
print 'PragPicker: Tsmooth larger than CF!' print 'PragPicker: Tsmooth larger than CF!'
return return
else: else:
for i in range(1, len(self.cf)): for i in range(1, len(self.cf)):
if i > ismooth: if i > ismooth:
ii1 = i - ismooth; ii1 = i - ismooth
cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[
else: ii1]) / ismooth
cfsmooth[i] = np.mean(self.cf[1 : i]) else:
cfsmooth[i] = np.mean(self.cf[1: i])
#select picking window # select picking window
#which is centered around tpick1 # which is centered around tpick1
ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \ ipick = np.where((self.Tcf >=
& (self.Tcf <= self.getpick1() + self.PickWindow / 2)) (self.getpick1() - self.PickWindow / 2)) and
cfipick = self.cf[ipick] - np.mean(self.cf[ipick]) (self.Tcf <=
Tcfpick = self.Tcf[ipick] (self.getpick1() + self.PickWindow / 2)))
cfsmoothipick = cfsmooth[ipick]- np.mean(self.cf[ipick]) cfipick = self.cf[ipick] - np.mean(self.cf[ipick])
ipick1 = np.argmin(abs(self.Tcf - self.getpick1())) Tcfpick = self.Tcf[ipick]
cfpick1 = 2 * self.cf[ipick1] cfsmoothipick = cfsmooth[ipick] - np.mean(self.cf[ipick])
ipick1 = np.argmin(abs(self.Tcf - self.getpick1()))
cfpick1 = 2 * self.cf[ipick1]
#check trend of CF, i.e. differences of CF and adjust aus regarding this trend # check trend of CF, i.e. differences of CF and adjust aus regarding this trend
#prominent trend: decrease aus # prominent trend: decrease aus
#flat: use given aus # flat: use given aus
cfdiff = np.diff(cfipick); cfdiff = np.diff(cfipick)
i0diff = np.where(cfdiff > 0) i0diff = np.where(cfdiff > 0)
cfdiff = cfdiff[i0diff] cfdiff = cfdiff[i0diff]
minaus = min(cfdiff * (1 + self.aus)); minaus = min(cfdiff * (1 + self.aus))
aus1 = max([minaus, self.aus]); aus1 = max([minaus, self.aus])
#at first we look to the right until the end of the pick window is reached # at first we look to the right until the end of the pick window is reached
flagpick_r = 0 flagpick_r = 0
flagpick_l = 0 flagpick_l = 0
flagpick = 0 flagpick = 0
lpickwindow = int(round(self.PickWindow / self.dt)) lpickwindow = int(round(self.PickWindow / self.dt))
for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])): for i in range(max(np.insert(ipick, 0, 2)),
if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]: min([ipick1 + lpickwindow + 1, len(self.cf) - 1])):
if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]: if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]:
if cfpick1 >= self.cf[i]: if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
pick_r = self.Tcf[i] if cfpick1 >= self.cf[i]:
self.Pick = pick_r pick_r = self.Tcf[i]
flagpick_l = 1 self.Pick = pick_r
cfpick_r = self.cf[i] flagpick_l = 1
break cfpick_r = self.cf[i]
break
#now we look to the left # now we look to the left
for i in range(ipick1, max([ipick1 - lpickwindow + 1, 2]), -1): for i in range(ipick1, max([ipick1 - lpickwindow + 1, 2]), -1):
if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]: if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]:
if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]: if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
if cfpick1 >= self.cf[i]: if cfpick1 >= self.cf[i]:
pick_l = self.Tcf[i] pick_l = self.Tcf[i]
self.Pick = pick_l self.Pick = pick_l
flagpick_r = 1 flagpick_r = 1
cfpick_l = self.cf[i] cfpick_l = self.cf[i]
break break
#now decide which pick: left or right? # now decide which pick: left or right?
if flagpick_l > 0 and flagpick_r > 0 and cfpick_l <= cfpick_r: if flagpick_l > 0 and flagpick_r > 0 and cfpick_l <= cfpick_r:
self.Pick = pick_l self.Pick = pick_l
elif flagpick_l > 0 and flagpick_r > 0 and cfpick_l >= cfpick_r: elif flagpick_l > 0 and flagpick_r > 0 and cfpick_l >= cfpick_r:
self.Pick = pick_r self.Pick = pick_r
if self.getiplot() > 1: if self.getiplot() > 1:
p = plt.figure(self.getiplot()) p = plt.figure(self.getiplot())
p1, = plt.plot(Tcfpick,cfipick, 'k') p1, = plt.plot(Tcfpick, cfipick, 'k')
p2, = plt.plot(Tcfpick,cfsmoothipick, 'r') p2, = plt.plot(Tcfpick, cfsmoothipick, 'r')
p3, = plt.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', linewidth=2) p3, = plt.plot([self.Pick, self.Pick],
plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick']) [min(cfipick), max(cfipick)], 'b', linewidth=2)
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick'])
plt.yticks([]) plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.title(self.Data[0].stats.station) plt.yticks([])
plt.show() plt.title(self.Data[0].stats.station)
raw_input() plt.show()
plt.close(p) raw_input()
plt.close(p)
else: else:
self.Pick = None self.Pick = None
print 'PragPicker: No initial onset time given! Check input!' print 'PragPicker: No initial onset time given! Check input!'
return

706
pylot/core/pick/run_autopicking.py
View File

@ -3,43 +3,38 @@
""" """
Function to run automated picking algorithms using AIC, Function to run automated picking algorithms using AIC,
HOS and AR prediction. Uses object CharFuns and Picker and HOS and AR prediction. Uses object CharFuns and Picker and
function conglomerate utils. function conglomerate utils.
:author: MAGS2 EP3 working group / Ludger Kueperkoch :author: MAGS2 EP3 working group / Ludger Kueperkoch
""" """
from obspy.core import read
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from pylot.core.pick.CharFuns import *
from pylot.core.pick.Picker import * from pylot.core.pick.Picker import *
from pylot.core.pick.CharFuns import * from pylot.core.pick.CharFuns import *
from pylot.core.pick import utils
def run_autopicking(wfstream, pickparam): def run_autopicking(wfstream, pickparam):
"""
'''
param: wfstream param: wfstream
:type: `~obspy.core.stream.Stream` :type: `~obspy.core.stream.Stream`
param: pickparam param: pickparam
:type: container of picking parameters from input file, :type: container of picking parameters from input file,
usually autoPyLoT.in usually autoPyLoT.in
''' """
# declaring pickparam variables (only for convenience) # declaring pickparam variables (only for convenience)
# read your autoPyLoT.in for details! # read your autoPyLoT.in for details!
#special parameters for P picking # special parameters for P picking
algoP = pickparam.getParam('algoP') algoP = pickparam.getParam('algoP')
iplot = pickparam.getParam('iplot') iplot = pickparam.getParam('iplot')
pstart = pickparam.getParam('pstart') pstart = pickparam.getParam('pstart')
pstop = pickparam.getParam('pstop') pstop = pickparam.getParam('pstop')
thosmw = pickparam.getParam('tlta') thosmw = pickparam.getParam('tlta')
hosorder = pickparam.getParam('hosorder') tsnrz = pickparam.getParam('tsnrz')
tsnrz = pickparam.getParam('tsnrz')
hosorder = pickparam.getParam('hosorder') hosorder = pickparam.getParam('hosorder')
bpz1 = pickparam.getParam('bpz1') bpz1 = pickparam.getParam('bpz1')
bpz2 = pickparam.getParam('bpz2') bpz2 = pickparam.getParam('bpz2')
@ -55,13 +50,13 @@ def run_autopicking(wfstream, pickparam):
minAICPslope = pickparam.getParam('minAICPslope') minAICPslope = pickparam.getParam('minAICPslope')
minAICPSNR = pickparam.getParam('minAICPSNR') minAICPSNR = pickparam.getParam('minAICPSNR')
timeerrorsP = pickparam.getParam('timeerrorsP') timeerrorsP = pickparam.getParam('timeerrorsP')
#special parameters for S picking # special parameters for S picking
algoS = pickparam.getParam('algoS') algoS = pickparam.getParam('algoS')
sstart = pickparam.getParam('sstart') sstart = pickparam.getParam('sstart')
sstop = pickparam.getParam('sstop') sstop = pickparam.getParam('sstop')
bph1 = pickparam.getParam('bph1') bph1 = pickparam.getParam('bph1')
bph2 = pickparam.getParam('bph2') bph2 = pickparam.getParam('bph2')
tsnrh = pickparam.getParam('tsnrh') tsnrh = pickparam.getParam('tsnrh')
pickwinS = pickparam.getParam('pickwinS') pickwinS = pickparam.getParam('pickwinS')
tpred1h = pickparam.getParam('tpred1h') tpred1h = pickparam.getParam('tpred1h')
tdet1h = pickparam.getParam('tdet1h') tdet1h = pickparam.getParam('tdet1h')
@ -76,7 +71,7 @@ def run_autopicking(wfstream, pickparam):
Srecalcwin = pickparam.getParam('Srecalcwin') Srecalcwin = pickparam.getParam('Srecalcwin')
nfacS = pickparam.getParam('nfacS') nfacS = pickparam.getParam('nfacS')
timeerrorsS = pickparam.getParam('timeerrorsS') timeerrorsS = pickparam.getParam('timeerrorsS')
#parameters for first-motion determination # parameters for first-motion determination
minFMSNR = pickparam.getParam('minFMSNR') minFMSNR = pickparam.getParam('minFMSNR')
fmpickwin = pickparam.getParam('fmpickwin') fmpickwin = pickparam.getParam('fmpickwin')
minfmweight = pickparam.getParam('minfmweight') minfmweight = pickparam.getParam('minfmweight')
@ -84,164 +79,192 @@ def run_autopicking(wfstream, pickparam):
# split components # split components
zdat = wfstream.select(component="Z") zdat = wfstream.select(component="Z")
edat = wfstream.select(component="E") edat = wfstream.select(component="E")
if len(edat) == 0: #check for other components if len(edat) == 0: # check for other components
edat = wfstream.select(component="2") edat = wfstream.select(component="2")
ndat = wfstream.select(component="N") ndat = wfstream.select(component="N")
if len(ndat) == 0: #check for other components if len(ndat) == 0: # check for other components
ndat = wfstream.select(component="1") ndat = wfstream.select(component="1")
if algoP == 'HOS' or algoP == 'ARZ' and zdat is not None: if algoP == 'HOS' or algoP == 'ARZ' and zdat is not None:
print '##########################################' print '##########################################'
print 'run_autopicking: Working on P onset of station %s' % zdat[0].stats.station print 'run_autopicking: Working on P onset of station %s' % zdat[
0].stats.station
print 'Filtering vertical trace ...' print 'Filtering vertical trace ...'
print zdat print zdat
z_copy = zdat.copy() z_copy = zdat.copy()
#filter and taper data # filter and taper data
tr_filt = zdat[0].copy() tr_filt = zdat[0].copy()
tr_filt.filter('bandpass', freqmin=bpz1[0], freqmax=bpz1[1], zerophase=False) tr_filt.filter('bandpass', freqmin=bpz1[0], freqmax=bpz1[1],
zerophase=False)
tr_filt.taper(max_percentage=0.05, type='hann') tr_filt.taper(max_percentage=0.05, type='hann')
z_copy[0].data = tr_filt.data z_copy[0].data = tr_filt.data
############################################################## ##############################################################
#check length of waveform and compare with cut times # check length of waveform and compare with cut times
Lc = pstop - pstart Lc = pstop - pstart
Lwf = zdat[0].stats.endtime - zdat[0].stats.starttime Lwf = zdat[0].stats.endtime - zdat[0].stats.starttime
Ldiff = Lwf - Lc Ldiff = Lwf - Lc
if Ldiff < 0: if Ldiff < 0:
print 'run_autopicking: Cutting times are too large for actual waveform!' print 'run_autopicking: Cutting times are too large for actual ' \
print 'Use entire waveform instead!' 'waveform!'
pstart = 0 print 'Use entire waveform instead!'
pstop = len(zdat[0].data) * zdat[0].stats.delta pstart = 0
pstop = len(zdat[0].data) * zdat[0].stats.delta
cuttimes = [pstart, pstop] cuttimes = [pstart, pstop]
if algoP == 'HOS': if algoP == 'HOS':
#calculate HOS-CF using subclass HOScf of class CharacteristicFunction # calculate HOS-CF using subclass HOScf of class
cf1 = HOScf(z_copy, cuttimes, thosmw, hosorder) #instance of HOScf # CharacteristicFunction
cf1 = HOScf(z_copy, cuttimes, thosmw, hosorder) # instance of HOScf
elif algoP == 'ARZ': elif algoP == 'ARZ':
#calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction # calculate ARZ-CF using subclass ARZcf of class
cf1 = ARZcf(z_copy, cuttimes, tpred1z, Parorder, tdet1z, addnoise) #instance of ARZcf # CharcteristicFunction
cf1 = ARZcf(z_copy, cuttimes, tpred1z, Parorder, tdet1z,
addnoise) # instance of ARZcf
############################################################## ##############################################################
#calculate AIC-HOS-CF using subclass AICcf of class CharacteristicFunction # calculate AIC-HOS-CF using subclass AICcf of class
#class needs stream object => build it # CharacteristicFunction
# class needs stream object => build it
tr_aic = tr_filt.copy() tr_aic = tr_filt.copy()
tr_aic.data =cf1.getCF() tr_aic.data = cf1.getCF()
z_copy[0].data = tr_aic.data z_copy[0].data = tr_aic.data
aiccf = AICcf(z_copy, cuttimes) #instance of AICcf aiccf = AICcf(z_copy, cuttimes) # instance of AICcf
############################################################## ##############################################################
#get prelimenary onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking # get prelimenary onset time from AIC-HOS-CF using subclass AICPicker
# of class AutoPicking
aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP) aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP)
############################################################## ##############################################################
#go on with processing if AIC onset passes quality control # go on with processing if AIC onset passes quality control
if aicpick.getSlope() >= minAICPslope and aicpick.getSNR() >= minAICPSNR: if (aicpick.getSlope() >= minAICPslope and
aicPflag = 1 aicpick.getSNR() >= minAICPSNR):
print 'AIC P-pick passes quality control: Slope: %f, SNR: %f' % \ aicPflag = 1
print 'AIC P-pick passes quality control: Slope: %f, SNR: %f' % \
(aicpick.getSlope(), aicpick.getSNR()) (aicpick.getSlope(), aicpick.getSNR())
print 'Go on with refined picking ...' print 'Go on with refined picking ...'
#re-filter waveform with larger bandpass # re-filter waveform with larger bandpass
print 'run_autopicking: re-filtering vertical trace ...' print 'run_autopicking: re-filtering vertical trace ...'
z_copy = zdat.copy() z_copy = zdat.copy()
tr_filt = zdat[0].copy() tr_filt = zdat[0].copy()
tr_filt.filter('bandpass', freqmin=bpz2[0], freqmax=bpz2[1], zerophase=False) tr_filt.filter('bandpass', freqmin=bpz2[0], freqmax=bpz2[1],
tr_filt.taper(max_percentage=0.05, type='hann') zerophase=False)
z_copy[0].data = tr_filt.data tr_filt.taper(max_percentage=0.05, type='hann')
############################################################# z_copy[0].data = tr_filt.data
#re-calculate CF from re-filtered trace in vicinity of initial onset #############################################################
cuttimes2 = [round(max([aicpick.getpick() - Precalcwin, 0])), \ # re-calculate CF from re-filtered trace in vicinity of initial
round(min([len(zdat[0].data) * zdat[0].stats.delta, \ # onset
aicpick.getpick() + Precalcwin]))] cuttimes2 = [round(max([aicpick.getpick() - Precalcwin, 0])),
if algoP == 'HOS': round(min([len(zdat[0].data) * zdat[0].stats.delta,
#calculate HOS-CF using subclass HOScf of class CharacteristicFunction aicpick.getpick() + Precalcwin]))]
cf2 = HOScf(z_copy, cuttimes2, thosmw, hosorder) #instance of HOScf if algoP == 'HOS':
elif algoP == 'ARZ': # calculate HOS-CF using subclass HOScf of class
#calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction # CharacteristicFunction
cf2 = ARZcf(z_copy, cuttimes2, tpred1z, Parorder, tdet1z, addnoise) #instance of ARZcf cf2 = HOScf(z_copy, cuttimes2, thosmw,
############################################################## hosorder) # instance of HOScf
#get refined onset time from CF2 using class Picker elif algoP == 'ARZ':
refPpick = PragPicker(cf2, tsnrz, pickwinP, iplot, ausP, tsmoothP, aicpick.getpick()) # calculate ARZ-CF using subclass ARZcf of class
############################################################# # CharcteristicFunction
#quality assessment cf2 = ARZcf(z_copy, cuttimes2, tpred1z, Parorder, tdet1z,
#get earliest and latest possible pick and symmetrized uncertainty addnoise) # instance of ARZcf
[lpickP, epickP, Perror] = earllatepicker(z_copy, nfacP, tsnrz, refPpick.getpick(), iplot) ##############################################################
# get refined onset time from CF2 using class Picker
refPpick = PragPicker(cf2, tsnrz, pickwinP, iplot, ausP, tsmoothP,
aicpick.getpick())
#############################################################
# quality assessment
# get earliest and latest possible pick and symmetrized uncertainty
[lpickP, epickP, Perror] = earllatepicker(z_copy, nfacP, tsnrz,
refPpick.getpick(), iplot)
#get SNR # get SNR
[SNRP, SNRPdB, Pnoiselevel] = getSNR(z_copy, tsnrz, refPpick.getpick()) [SNRP, SNRPdB, Pnoiselevel] = getSNR(z_copy, tsnrz,
refPpick.getpick())
#weight P-onset using symmetric error # weight P-onset using symmetric error
if Perror <= timeerrorsP[0]: if Perror <= timeerrorsP[0]:
Pweight = 0 Pweight = 0
elif Perror > timeerrorsP[0] and Perror <= timeerrorsP[1]: elif timeerrorsP[0] < Perror <= timeerrorsP[1]:
Pweight = 1 Pweight = 1
elif Perror > timeerrorsP[1] and Perror <= timeerrorsP[2]: elif timeerrorsP[1] < Perror <= timeerrorsP[2]:
Pweight = 2 Pweight = 2
elif Perror > timeerrorsP[2] and Perror <= timeerrorsP[3]: elif timeerrorsP[2] < Perror <= timeerrorsP[3]:
Pweight = 3 Pweight = 3
elif Perror > timeerrorsP[3]: elif Perror > timeerrorsP[3]:
Pweight = 4 Pweight = 4
##############################################################
# get first motion of P onset
# certain quality required
if Pweight <= minfmweight and SNRP >= minFMSNR:
FM = fmpicker(zdat, z_copy, fmpickwin, refPpick.getpick(),
iplot)
else:
FM = 'N'
print 'run_autopicking: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \
'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM)
##############################################################
#get first motion of P onset
#certain quality required
if Pweight <= minfmweight and SNRP >= minFMSNR:
FM = fmpicker(zdat, z_copy, fmpickwin, refPpick.getpick(), iplot)
else:
FM = 'N'
print 'run_autopicking: P-weight: %d, SNR: %f, SNR[dB]: %f, Polarity: %s' % (Pweight, SNRP, SNRPdB, FM)
else: else:
print 'Bad initial (AIC) P-pick, skip this onset!' print 'Bad initial (AIC) P-pick, skip this onset!'
print 'AIC-SNR=', aicpick.getSNR(), 'AIC-Slope=', aicpick.getSlope() print 'AIC-SNR=', aicpick.getSNR(), 'AIC-Slope=', aicpick.getSlope()
Pweight = 4 Pweight = 4
Sweight = 4 Sweight = 4
FM = 'N' FM = 'N'
SNRP = None SNRP = None
SNRPdB = None SNRPdB = None
SNRS = None SNRS = None
SNRSdB = None SNRSdB = None
aicSflag = 0 aicSflag = 0
aicPflag = 0 aicPflag = 0
else: else:
print 'run_autopicking: No vertical component data available, skipping station!' print 'run_autopicking: No vertical component data available, ' \
return 'skipping station!'
return
if edat is not None and ndat is not None and len(edat) > 0 and len(ndat) > 0 and Pweight < 4: if edat is not None and ndat is not None and len(edat) > 0 and len(
print 'Go on picking S onset ...' ndat) > 0 and Pweight < 4:
print 'Go on picking S onset ...'
print '##################################################' print '##################################################'
print 'Working on S onset of station %s' % edat[0].stats.station print 'Working on S onset of station %s' % edat[0].stats.station
print 'Filtering horizontal traces ...' print 'Filtering horizontal traces ...'
#determine time window for calculating CF after P onset # determine time window for calculating CF after P onset
#cuttimesh = [round(refPpick.getpick() + sstart), round(refPpick.getpick() + sstop)] # cuttimesh = [round(refPpick.getpick() + sstart),
cuttimesh = [round(max([refPpick.getpick() + sstart, 0])), \ # round(refPpick.getpick() + sstop)]
cuttimesh = [round(max([refPpick.getpick() + sstart, 0])),
round(min([refPpick.getpick() + sstop, Lwf]))] round(min([refPpick.getpick() + sstop, Lwf]))]
if algoS == 'ARH': if algoS == 'ARH':
print edat, ndat print edat, ndat
#re-create stream object including both horizontal components # re-create stream object including both horizontal components
hdat = edat.copy() hdat = edat.copy()
hdat += ndat hdat += ndat
h_copy = hdat.copy() h_copy = hdat.copy()
#filter and taper data # filter and taper data
trH1_filt = hdat[0].copy() trH1_filt = hdat[0].copy()
trH2_filt = hdat[1].copy() trH2_filt = hdat[1].copy()
trH1_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1], zerophase=False) trH1_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1],
trH2_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1], zerophase=False) zerophase=False)
trH2_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1],
zerophase=False)
trH1_filt.taper(max_percentage=0.05, type='hann') trH1_filt.taper(max_percentage=0.05, type='hann')
trH2_filt.taper(max_percentage=0.05, type='hann') trH2_filt.taper(max_percentage=0.05, type='hann')
h_copy[0].data = trH1_filt.data h_copy[0].data = trH1_filt.data
h_copy[1].data = trH2_filt.data h_copy[1].data = trH2_filt.data
elif algoS == 'AR3': elif algoS == 'AR3':
print zdat, edat, ndat print zdat, edat, ndat
#re-create stream object including both horizontal components # re-create stream object including both horizontal components
hdat = zdat.copy() hdat = zdat.copy()
hdat += edat hdat += edat
hdat += ndat hdat += ndat
h_copy = hdat.copy() h_copy = hdat.copy()
#filter and taper data # filter and taper data
trH1_filt = hdat[0].copy() trH1_filt = hdat[0].copy()
trH2_filt = hdat[1].copy() trH2_filt = hdat[1].copy()
trH3_filt = hdat[2].copy() trH3_filt = hdat[2].copy()
trH1_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1], zerophase=False) trH1_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1],
trH2_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1], zerophase=False) zerophase=False)
trH3_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1], zerophase=False) trH2_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1],
zerophase=False)
trH3_filt.filter('bandpass', freqmin=bph1[0], freqmax=bph1[1],
zerophase=False)
trH1_filt.taper(max_percentage=0.05, type='hann') trH1_filt.taper(max_percentage=0.05, type='hann')
trH2_filt.taper(max_percentage=0.05, type='hann') trH2_filt.taper(max_percentage=0.05, type='hann')
trH3_filt.taper(max_percentage=0.05, type='hann') trH3_filt.taper(max_percentage=0.05, type='hann')
@ -250,210 +273,289 @@ def run_autopicking(wfstream, pickparam):
h_copy[2].data = trH3_filt.data h_copy[2].data = trH3_filt.data
############################################################## ##############################################################
if algoS == 'ARH': if algoS == 'ARH':
#calculate ARH-CF using subclass ARHcf of class CharcteristicFunction # calculate ARH-CF using subclass ARHcf of class
arhcf1 = ARHcf(h_copy, cuttimesh, tpred1h, Sarorder, tdet1h, addnoise) #instance of ARHcf # CharcteristicFunction
arhcf1 = ARHcf(h_copy, cuttimesh, tpred1h, Sarorder, tdet1h,
addnoise) # instance of ARHcf
elif algoS == 'AR3': elif algoS == 'AR3':
#calculate ARH-CF using subclass AR3cf of class CharcteristicFunction # calculate ARH-CF using subclass AR3cf of class
arhcf1 = AR3Ccf(h_copy, cuttimesh, tpred1h, Sarorder, tdet1h, addnoise) #instance of ARHcf # CharcteristicFunction
arhcf1 = AR3Ccf(h_copy, cuttimesh, tpred1h, Sarorder, tdet1h,
addnoise) # instance of ARHcf
############################################################## ##############################################################
#calculate AIC-ARH-CF using subclass AICcf of class CharacteristicFunction # calculate AIC-ARH-CF using subclass AICcf of class
#class needs stream object => build it # CharacteristicFunction
# class needs stream object => build it
tr_arhaic = trH1_filt.copy() tr_arhaic = trH1_filt.copy()
tr_arhaic.data = arhcf1.getCF() tr_arhaic.data = arhcf1.getCF()
h_copy[0].data = tr_arhaic.data h_copy[0].data = tr_arhaic.data
#calculate ARH-AIC-CF # calculate ARH-AIC-CF
haiccf = AICcf(h_copy, cuttimesh) #instance of AICcf haiccf = AICcf(h_copy, cuttimesh) # instance of AICcf
############################################################## ##############################################################
#get prelimenary onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking # get prelimenary onset time from AIC-HOS-CF using subclass AICPicker
aicarhpick = AICPicker(haiccf, tsnrh, pickwinS, iplot, None, aictsmoothS) # of class AutoPicking
aicarhpick = AICPicker(haiccf, tsnrh, pickwinS, iplot, None,
aictsmoothS)
############################################################### ###############################################################
#go on with processing if AIC onset passes quality control # go on with processing if AIC onset passes quality control
if aicarhpick.getSlope() >= minAICSslope and aicarhpick.getSNR() >= minAICSSNR: if (aicarhpick.getSlope() >= minAICSslope and
aicSflag = 1 aicarhpick.getSNR() >= minAICSSNR):
print 'AIC S-pick passes quality control: Slope: %f, SNR: %f' \ aicSflag = 1
% (aicarhpick.getSlope(), aicarhpick.getSNR()) print 'AIC S-pick passes quality control: Slope: %f, SNR: %f' \
print 'Go on with refined picking ...' % (aicarhpick.getSlope(), aicarhpick.getSNR())
#re-calculate CF from re-filtered trace in vicinity of initial onset print 'Go on with refined picking ...'
cuttimesh2 = [round(aicarhpick.getpick() - Srecalcwin), \ # re-calculate CF from re-filtered trace in vicinity of initial
round(aicarhpick.getpick() + Srecalcwin)] # onset
#re-filter waveform with larger bandpass cuttimesh2 = [round(aicarhpick.getpick() - Srecalcwin),
print 'run_autopicking: re-filtering horizontal traces...' round(aicarhpick.getpick() + Srecalcwin)]
h_copy = hdat.copy() # re-filter waveform with larger bandpass
#filter and taper data print 'run_autopicking: re-filtering horizontal traces...'
if algoS == 'ARH': h_copy = hdat.copy()
trH1_filt = hdat[0].copy() # filter and taper data
trH2_filt = hdat[1].copy() if algoS == 'ARH':
trH1_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1], zerophase=False) trH1_filt = hdat[0].copy()
trH2_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1], zerophase=False) trH2_filt = hdat[1].copy()
trH1_filt.taper(max_percentage=0.05, type='hann') trH1_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1],
trH2_filt.taper(max_percentage=0.05, type='hann') zerophase=False)
h_copy[0].data = trH1_filt.data trH2_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1],
h_copy[1].data = trH2_filt.data zerophase=False)
############################################################# trH1_filt.taper(max_percentage=0.05, type='hann')
arhcf2 = ARHcf(h_copy, cuttimesh2, tpred2h, Sarorder, tdet2h, addnoise) #instance of ARHcf trH2_filt.taper(max_percentage=0.05, type='hann')
elif algoS == 'AR3': h_copy[0].data = trH1_filt.data
trH1_filt = hdat[0].copy() h_copy[1].data = trH2_filt.data
trH2_filt = hdat[1].copy() #############################################################
trH3_filt = hdat[2].copy() arhcf2 = ARHcf(h_copy, cuttimesh2, tpred2h, Sarorder, tdet2h,
trH1_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1], zerophase=False) addnoise) # instance of ARHcf
trH2_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1], zerophase=False) elif algoS == 'AR3':
trH3_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1], zerophase=False) trH1_filt = hdat[0].copy()
trH1_filt.taper(max_percentage=0.05, type='hann') trH2_filt = hdat[1].copy()
trH2_filt.taper(max_percentage=0.05, type='hann') trH3_filt = hdat[2].copy()
trH3_filt.taper(max_percentage=0.05, type='hann') trH1_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1],
h_copy[0].data = trH1_filt.data zerophase=False)
h_copy[1].data = trH2_filt.data trH2_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1],
h_copy[2].data = trH3_filt.data zerophase=False)
############################################################# trH3_filt.filter('bandpass', freqmin=bph2[0], freqmax=bph2[1],
arhcf2 = AR3Ccf(h_copy, cuttimesh2, tpred2h, Sarorder, tdet2h, addnoise) #instance of ARHcf zerophase=False)
trH1_filt.taper(max_percentage=0.05, type='hann')
trH2_filt.taper(max_percentage=0.05, type='hann')
trH3_filt.taper(max_percentage=0.05, type='hann')
h_copy[0].data = trH1_filt.data
h_copy[1].data = trH2_filt.data
h_copy[2].data = trH3_filt.data
#############################################################
arhcf2 = AR3Ccf(h_copy, cuttimesh2, tpred2h, Sarorder, tdet2h,
addnoise) # instance of ARHcf
#get refined onset time from CF2 using class Picker # get refined onset time from CF2 using class Picker
refSpick = PragPicker(arhcf2, tsnrh, pickwinS, iplot, ausS, tsmoothS, aicarhpick.getpick()) refSpick = PragPicker(arhcf2, tsnrh, pickwinS, iplot, ausS,
############################################################# tsmoothS, aicarhpick.getpick())
#quality assessment #############################################################
#get earliest and latest possible pick and symmetrized uncertainty # quality assessment
h_copy[0].data = trH1_filt.data # get earliest and latest possible pick and symmetrized uncertainty
[lpickS1, epickS1, Serror1] = earllatepicker(h_copy, nfacS, tsnrh, refSpick.getpick(), iplot) h_copy[0].data = trH1_filt.data
h_copy[0].data = trH2_filt.data [lpickS1, epickS1, Serror1] = earllatepicker(h_copy, nfacS, tsnrh,
[lpickS2, epickS2, Serror2] = earllatepicker(h_copy, nfacS, tsnrh, refSpick.getpick(), iplot) refSpick.getpick(),
if algoS == 'ARH': iplot)
#get earliest pick of both earliest possible picks h_copy[0].data = trH2_filt.data
epick = [epickS1, epickS2] [lpickS2, epickS2, Serror2] = earllatepicker(h_copy, nfacS, tsnrh,
lpick = [lpickS1, lpickS2] refSpick.getpick(),
pickerr = [Serror1, Serror2] iplot)
ipick =np.argmin([epickS1, epickS2]) if algoS == 'ARH':
elif algoS == 'AR3': # get earliest pick of both earliest possible picks
[lpickS3, epickS3, Serror3] = earllatepicker(h_copy, nfacS, tsnrh, refSpick.getpick(), iplot) epick = [epickS1, epickS2]
#get earliest pick of all three picks lpick = [lpickS1, lpickS2]
epick = [epickS1, epickS2, epickS3] pickerr = [Serror1, Serror2]
lpick = [lpickS1, lpickS2, lpickS3] ipick = np.argmin([epickS1, epickS2])
pickerr = [Serror1, Serror2, Serror3] elif algoS == 'AR3':
ipick =np.argmin([epickS1, epickS2, epickS3]) [lpickS3, epickS3, Serror3] = earllatepicker(h_copy, nfacS,
epickS = epick[ipick] tsnrh,
lpickS = lpick[ipick] refSpick.getpick(),
Serror = pickerr[ipick] iplot)
# get earliest pick of all three picks
epick = [epickS1, epickS2, epickS3]
lpick = [lpickS1, lpickS2, lpickS3]
pickerr = [Serror1, Serror2, Serror3]
ipick = np.argmin([epickS1, epickS2, epickS3])
epickS = epick[ipick]
lpickS = lpick[ipick]
Serror = pickerr[ipick]
#get SNR # get SNR
[SNRS, SNRSdB, Snoiselevel] = getSNR(h_copy, tsnrh, refSpick.getpick()) [SNRS, SNRSdB, Snoiselevel] = getSNR(h_copy, tsnrh,
refSpick.getpick())
#weight S-onset using symmetric error # weight S-onset using symmetric error
if Serror <= timeerrorsS[0]: if Serror <= timeerrorsS[0]:
Sweight = 0 Sweight = 0
elif Serror > timeerrorsS[0] and Serror <= timeerrorsS[1]: elif timeerrorsS[0] < Serror <= timeerrorsS[1]:
Sweight = 1 Sweight = 1
elif Perror > timeerrorsS[1] and Serror <= timeerrorsS[2]: elif Perror > timeerrorsS[1] and Serror <= timeerrorsS[2]:
Sweight = 2 Sweight = 2
elif Serror > timeerrorsS[2] and Serror <= timeerrorsS[3]: elif timeerrorsS[2] < Serror <= timeerrorsS[3]:
Sweight = 3 Sweight = 3
elif Serror > timeerrorsS[3]: elif Serror > timeerrorsS[3]:
Sweight = 4 Sweight = 4
print 'run_autopicking: S-weight: %d, SNR: %f, SNR[dB]: %f' % (Sweight, SNRS, SNRSdB) print 'run_autopicking: S-weight: %d, SNR: %f, SNR[dB]: %f' % (
Sweight, SNRS, SNRSdB)
else: else:
print 'Bad initial (AIC) S-pick, skip this onset!' print 'Bad initial (AIC) S-pick, skip this onset!'
print 'AIC-SNR=', aicarhpick.getSNR(), 'AIC-Slope=', aicarhpick.getSlope() print 'AIC-SNR=', aicarhpick.getSNR(), \
Sweight = 4 'AIC-Slope=', aicarhpick.getSlope()
SNRS = None Sweight = 4
SNRSdB = None SNRS = None
aicSflag = 0 SNRSdB = None
aicSflag = 0
else: else:
print 'run_autopicking: No horizontal component data available or bad P onset, skipping S picking!' print 'run_autopicking: No horizontal component data available or ' \
return 'bad P onset, skipping S picking!'
return
############################################################## ##############################################################
if iplot > 0: if iplot > 0:
#plot vertical trace # plot vertical trace
plt.figure() plt.figure()
plt.subplot(3,1,1) plt.subplot(3, 1, 1)
tdata = np.arange(0, zdat[0].stats.npts / tr_filt.stats.sampling_rate, tr_filt.stats.delta) tdata = np.arange(0, zdat[0].stats.npts / tr_filt.stats.sampling_rate,
#check equal length of arrays, sometimes they are different!? tr_filt.stats.delta)
wfldiff = len(tr_filt.data) - len(tdata) # check equal length of arrays, sometimes they are different!?
if wfldiff < 0: wfldiff = len(tr_filt.data) - len(tdata)
tdata = tdata[0:len(tdata) - abs(wfldiff)] if wfldiff < 0:
p1, = plt.plot(tdata, tr_filt.data/max(tr_filt.data), 'k') tdata = tdata[0:len(tdata) - abs(wfldiff)]
if Pweight < 4: p1, = plt.plot(tdata, tr_filt.data / max(tr_filt.data), 'k')
p2, = plt.plot(cf1.getTimeArray(), cf1.getCF() / max(cf1.getCF()), 'b') if Pweight < 4:
if aicPflag == 1: p2, = plt.plot(cf1.getTimeArray(), cf1.getCF() / max(cf1.getCF()),
p3, = plt.plot(cf2.getTimeArray(), cf2.getCF() / max(cf2.getCF()), 'm') 'b')
p4, = plt.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1], 'r') if aicPflag == 1:
plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [1, 1], 'r') p3, = plt.plot(cf2.getTimeArray(),
plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [-1, -1], 'r') cf2.getCF() / max(cf2.getCF()), 'm')
p5, = plt.plot([refPpick.getpick(), refPpick.getpick()], [-1.3, 1.3], 'r', linewidth=2) p4, = plt.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1],
plt.plot([refPpick.getpick()-0.5, refPpick.getpick()+0.5], [1.3, 1.3], 'r', linewidth=2) 'r')
plt.plot([refPpick.getpick()-0.5, refPpick.getpick()+0.5], [-1.3, -1.3], 'r', linewidth=2) plt.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5],
plt.plot([lpickP, lpickP], [-1.1, 1.1], 'r--') [1, 1], 'r')
plt.plot([epickP, epickP], [-1.1, 1.1], 'r--') plt.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5],
plt.legend([p1, p2, p3, p4, p5], ['Data', 'CF1', 'CF2', 'Initial P Onset', 'Final P Pick']) [-1, -1], 'r')
plt.title('%s, %s, P Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f Polarity: %s' % (tr_filt.stats.station, \ p5, = plt.plot([refPpick.getpick(), refPpick.getpick()],
tr_filt.stats.channel, Pweight, SNRP, SNRPdB, FM)) [-1.3, 1.3], 'r', linewidth=2)
else: plt.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5],
plt.legend([p1, p2], ['Data', 'CF1']) [1.3, 1.3], 'r', linewidth=2)
plt.title('%s, P Weight=%d, SNR=None, SNRdB=None' % (tr_filt.stats.channel, Pweight)) plt.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5],
plt.yticks([]) [-1.3, -1.3], 'r', linewidth=2)
plt.ylim([-1.5, 1.5]) plt.plot([lpickP, lpickP], [-1.1, 1.1], 'r--')
plt.ylabel('Normalized Counts') plt.plot([epickP, epickP], [-1.1, 1.1], 'r--')
plt.suptitle(tr_filt.stats.starttime) plt.legend([p1, p2, p3, p4, p5],
['Data', 'CF1', 'CF2', 'Initial P Onset',
'Final P Pick'])
plt.title('%s, %s, P Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f '
'Polarity: %s' % (tr_filt.stats.station,
tr_filt.stats.channel,
Pweight,
SNRP,
SNRPdB,
FM))
else:
plt.legend([p1, p2], ['Data', 'CF1'])
plt.title('%s, P Weight=%d, SNR=None, '
'SNRdB=None' % (tr_filt.stats.channel, Pweight))
plt.yticks([])
plt.ylim([-1.5, 1.5])
plt.ylabel('Normalized Counts')
plt.suptitle(tr_filt.stats.starttime)
#plot horizontal traces # plot horizontal traces
plt.subplot(3,1,2) plt.subplot(3, 1, 2)
th1data = np.arange(0, trH1_filt.stats.npts / trH1_filt.stats.sampling_rate, trH1_filt.stats.delta) th1data = np.arange(0,
#check equal length of arrays, sometimes they are different!? trH1_filt.stats.npts /
wfldiff = len(trH1_filt.data) - len(th1data) trH1_filt.stats.sampling_rate,
if wfldiff < 0: trH1_filt.stats.delta)
th1data = th1data[0:len(th1data) - abs(wfldiff)] # check equal length of arrays, sometimes they are different!?
p21, = plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k') wfldiff = len(trH1_filt.data) - len(th1data)
if Pweight < 4: if wfldiff < 0:
p22, = plt.plot(arhcf1.getTimeArray(), arhcf1.getCF()/max(arhcf1.getCF()), 'b') th1data = th1data[0:len(th1data) - abs(wfldiff)]
if aicSflag == 1: p21, = plt.plot(th1data, trH1_filt.data / max(trH1_filt.data), 'k')
p23, = plt.plot(arhcf2.getTimeArray(), arhcf2.getCF()/max(arhcf2.getCF()), 'm') if Pweight < 4:
p24, = plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'g') p22, = plt.plot(arhcf1.getTimeArray(),
plt.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [1, 1], 'g') arhcf1.getCF() / max(arhcf1.getCF()), 'b')
plt.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [-1, -1], 'g') if aicSflag == 1:
p25, = plt.plot([refSpick.getpick(), refSpick.getpick()], [-1.3, 1.3], 'g', linewidth=2) p23, = plt.plot(arhcf2.getTimeArray(),
plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [1.3, 1.3], 'g', linewidth=2) arhcf2.getCF() / max(arhcf2.getCF()), 'm')
plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [-1.3, -1.3], 'g', linewidth=2) p24, = plt.plot([aicarhpick.getpick(), aicarhpick.getpick()],
plt.plot([lpickS, lpickS], [-1.1, 1.1], 'g--') [-1, 1], 'g')
plt.plot([epickS, epickS], [-1.1, 1.1], 'g--') plt.plot(
plt.legend([p21, p22, p23, p24, p25], ['Data', 'CF1', 'CF2', 'Initial S Onset', 'Final S Pick']) [aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5],
plt.title('%s, S Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f' % (trH1_filt.stats.channel, \ [1, 1], 'g')
Sweight, SNRS, SNRSdB)) plt.plot(
else: [aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5],
plt.legend([p21, p22], ['Data', 'CF1']) [-1, -1], 'g')
plt.title('%s, S Weight=%d, SNR=None, SNRdB=None' % (trH1_filt.stats.channel, Sweight)) p25, = plt.plot([refSpick.getpick(), refSpick.getpick()],
plt.yticks([]) [-1.3, 1.3], 'g', linewidth=2)
plt.ylim([-1.5, 1.5]) plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
plt.ylabel('Normalized Counts') [1.3, 1.3], 'g', linewidth=2)
plt.suptitle(trH1_filt.stats.starttime) plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[-1.3, -1.3], 'g', linewidth=2)
plt.plot([lpickS, lpickS], [-1.1, 1.1], 'g--')
plt.plot([epickS, epickS], [-1.1, 1.1], 'g--')
plt.legend([p21, p22, p23, p24, p25],
['Data', 'CF1', 'CF2', 'Initial S Onset',
'Final S Pick'])
plt.title('%s, S Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f' % (
trH1_filt.stats.channel,
Sweight, SNRS, SNRSdB))
else:
plt.legend([p21, p22], ['Data', 'CF1'])
plt.title('%s, S Weight=%d, SNR=None, SNRdB=None' % (
trH1_filt.stats.channel, Sweight))
plt.yticks([])
plt.ylim([-1.5, 1.5])
plt.ylabel('Normalized Counts')
plt.suptitle(trH1_filt.stats.starttime)
plt.subplot(3,1,3) plt.subplot(3, 1, 3)
th2data = np.arange(0, trH2_filt.stats.npts / trH2_filt.stats.sampling_rate, trH2_filt.stats.delta) th2data = np.arange(0,
#check equal length of arrays, sometimes they are different!? trH2_filt.stats.npts /
wfldiff = len(trH2_filt.data) - len(th2data) trH2_filt.stats.sampling_rate,
if wfldiff < 0: trH2_filt.stats.delta)
th2data = th2data[0:len(th2data) - abs(wfldiff)] # check equal length of arrays, sometimes they are different!?
plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k') wfldiff = len(trH2_filt.data) - len(th2data)
if Pweight < 4: if wfldiff < 0:
p22, = plt.plot(arhcf1.getTimeArray(), arhcf1.getCF()/max(arhcf1.getCF()), 'b') th2data = th2data[0:len(th2data) - abs(wfldiff)]
if aicSflag == 1: plt.plot(th2data, trH2_filt.data / max(trH2_filt.data), 'k')
p23, = plt.plot(arhcf2.getTimeArray(), arhcf2.getCF()/max(arhcf2.getCF()), 'm') if Pweight < 4:
p24, = plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'g') p22, = plt.plot(arhcf1.getTimeArray(),
plt.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [1, 1], 'g') arhcf1.getCF() / max(arhcf1.getCF()), 'b')
plt.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [-1, -1], 'g') if aicSflag == 1:
p25, = plt.plot([refSpick.getpick(), refSpick.getpick()], [-1.3, 1.3], 'g', linewidth=2) p23, = plt.plot(arhcf2.getTimeArray(),
plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [1.3, 1.3], 'g', linewidth=2) arhcf2.getCF() / max(arhcf2.getCF()), 'm')
plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [-1.3, -1.3], 'g', linewidth=2) p24, = plt.plot([aicarhpick.getpick(), aicarhpick.getpick()],
plt.plot([lpickS, lpickS], [-1.1, 1.1], 'g--') [-1, 1], 'g')
plt.plot([epickS, epickS], [-1.1, 1.1], 'g--') plt.plot(
plt.legend([p21, p22, p23, p24, p25], ['Data', 'CF1', 'CF2', 'Initial S Onset', 'Final S Pick']) [aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5],
else: [1, 1], 'g')
plt.legend([p21, p22], ['Data', 'CF1']) plt.plot(
plt.yticks([]) [aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5],
plt.ylim([-1.5, 1.5]) [-1, -1], 'g')
plt.xlabel('Time [s] after %s' % tr_filt.stats.starttime) p25, = plt.plot([refSpick.getpick(), refSpick.getpick()],
plt.ylabel('Normalized Counts') [-1.3, 1.3], 'g', linewidth=2)
plt.title(trH2_filt.stats.channel) plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
plt.show() [1.3, 1.3], 'g', linewidth=2)
raw_input() plt.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
plt.close() [-1.3, -1.3], 'g', linewidth=2)
plt.plot([lpickS, lpickS], [-1.1, 1.1], 'g--')
plt.plot([epickS, epickS], [-1.1, 1.1], 'g--')
plt.legend([p21, p22, p23, p24, p25],
['Data', 'CF1', 'CF2', 'Initial S Onset',
'Final S Pick'])
else:
plt.legend([p21, p22], ['Data', 'CF1'])
plt.yticks([])
plt.ylim([-1.5, 1.5])
plt.xlabel('Time [s] after %s' % tr_filt.stats.starttime)
plt.ylabel('Normalized Counts')
plt.title(trH2_filt.stats.channel)
plt.show()
raw_input()
plt.close()