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

Switched off warnings.

Browse Source
This commit is contained in:
Ludger Küperkoch 2015-06-22 15:35:16 +02:00
parent 4a911a4ac9
commit f2510ff400
1 changed files with 215 additions and 234 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -3,65 +3,62 @@
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, Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group, 2010:
2010: Automated determination of P-phase arrival times at regional and local Automated determination of P-phase arrival times at regional and local distances
distances using higher order statistics, Geophys. J. Int., 181, 1159-1170 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, autoregressive prediction: application ot local and regional distances, Geophys. J. Int.,
Geophys. J. Int., 188, 687-702. 188, 687-702.
The picks with the above described algorithms are assumed to be the most likely The picks with the above described algorithms are assumed to be the most likely picks.
picks. For each most likely pick the corresponding earliest and latest possible For each most likely pick the corresponding earliest and latest possible picks are
picks are calculated after Diehl & Kissling (2009). 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
import warnings
class AutoPicking(object): class AutoPicking(object):
''' '''
Superclass of different, automated picking algorithms applied on a CF Superclass of different, automated picking algorithms applied on a CF determined
determined using AIC, HOS, or AR prediction. using AIC, HOS, or AR prediction.
''' '''
def __init__(self, cf, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None, warnings.simplefilter('ignore')
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 :param: cf, characteristic function, on which the picking algorithm is applied
applied :type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
:type cf: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
:param TSNR: length of time windows for SNR determination - [s] :param: TSNR, length of time windows around pick used to determine SNR [s]
:type TSNR: tuple (T_noise, T_gap, T_signal) :type: tuple (T_noise, T_gap, T_signal)
:param PickWindow: length of pick window - [s] :param: PickWindow, length of pick window [s]
:type PickWindow: float :type: float
:param iplot: no. of figure window for plotting interims results :param: iplot, no. of figure window for plotting interims results
:type iplot: integer :type: integer
:param aus: aus ("artificial uplift of samples"), find local minimum at :param: aus ("artificial uplift of samples"), find local minimum at i if aic(i-1)*(1+aus) >= aic(i)
i if aic(i-1)*(1+aus) >= aic(i) :type: float
:type aus: float
:param Tsmooth: length of moving window to calculate smoothed CF - [s] :param: Tsmooth, length of moving smoothing window to calculate smoothed CF [s]
:type Tsmooth: float :type: float
:param Pick1: initial (prelimenary) onset time, starting point for :param: Pick1, initial (prelimenary) onset time, starting point for PragPicker and
PragPicker EarlLatePicker
:type Pick1: float :type: float
''' '''
assert isinstance(cf, assert isinstance(cf, CharacteristicFunction), "%s is not a CharacteristicFunction object" % str(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()
@ -87,8 +84,9 @@ 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
@ -118,7 +116,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
@ -142,156 +140,144 @@ 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[ aicsmooth[i] = aicsmooth[i - 1] + (aic[i] - aic[ii1]) / ismooth
ii1]) / ismooth else:
else: aicsmooth[i] = np.mean(aic[1 : i])
aicsmooth[i] = np.mean(aic[1: i]) #remove offset
# 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 #get maximum of 1st derivative of AIC-CF (more stable!) as starting point
# 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 #find minimum in AIC-CF front of maximum
lpickwindow = int(round(self.PickWindow / self.dt)) lpickwindow = int(round(self.PickWindow / self.dt))
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1): for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1):
if aicsmooth[i - 1] >= aicsmooth[i]: if aicsmooth[i - 1] >= aicsmooth[i]:
self.Pick = self.Tcf[i] self.Pick = self.Tcf[i]
break break
# if no minimum could be found: #if no minimum could be found:
# search in 1st derivative of AIC-CF #search in 1st derivative of AIC-CF
if self.Pick is None: if self.Pick is None:
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1): for i in range(icfmax -1, max([icfmax -lpickwindow, 2]), -1):
if diffcf[i - 1] >= diffcf[i]: if diffcf[i -1] >= diffcf[i]:
self.Pick = self.Tcf[i] self.Pick = self.Tcf[i]
break break
# quality assessment using SNR and slope from CF # quality assessment using SNR and slope from CF
if self.Pick is not None: if self.Pick is not None:
# get noise window # get noise window
inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1])
self.TSNR[1]) # check, if these are counts or m/s, important for slope estimation!
# check, if these are counts or m/s, important for slope estimation! # this is quick and dirty, better solution?
# this is quick and dirty, better solution? if max(self.Data[0].data < 1e-3):
if max(self.Data[0].data < 1e-3): self.Data[0].data = self.Data[0].data * 1000000
self.Data[0].data *= 1000000 # get signal window
# get signal window isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2])
isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2]) # calculate SNR from CF
# calculate SNR from CF self.SNR = max(abs(aic[isignal] - np.mean(aic[isignal]))) / max(abs(aic[inoise] \
self.SNR = max(abs(aic[isignal] - np.mean(aic[isignal]))) / \ - np.mean(aic[inoise])))
max(abs(aic[inoise] - np.mean(aic[inoise]))) # calculate slope from CF after initial pick
# calculate slope from CF after initial pick # get slope window
# get slope window tslope = self.TSNR[3] #slope determination window
tslope = self.TSNR[3] # slope determination window islope = np.where((self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)])) \
islope = np.where( & (self.Tcf >= self.Pick))
(self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)])) # find maximum within slope determination window
and (self.Tcf >= self.Pick)) # 'cause slope should be calculated up to first local minimum only!
# find maximum within slope determination window imax = np.argmax(self.Data[0].data[islope])
# 'cause slope should be calculated up to first local minimum only! if imax == 0:
imax = np.argmax(self.Data[0].data[islope]) print 'AICPicker: Maximum for slope determination right at the beginning of the window!'
if imax == 0: print 'Choose longer slope determination window!'
print 'AICPicker: Maximum for slope determination right at ' \ return
'the beginning of the window!' islope = islope[0][0 :imax]
print 'Choose longer slope determination window!' dataslope = self.Data[0].data[islope]
return # calculate slope as polynomal fit of order 1
islope = islope[0][0:imax] xslope = np.arange(0, len(dataslope), 1)
dataslope = self.Data[0].data[islope] P = np.polyfit(xslope, dataslope, 1)
# calculate slope as polynomal fit of order 1 datafit = np.polyval(P, xslope)
xslope = np.arange(0, len(dataslope), 1) if datafit[0] >= datafit[len(datafit) - 1]:
P = np.polyfit(xslope, dataslope, 1) print 'AICPicker: Negative slope, bad onset skipped!'
datafit = np.polyval(P, xslope) return
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] 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', p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2)
linewidth=2) plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick'])
plt.legend([p1, p2, p3], else:
['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick']) plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
else: plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) plt.yticks([])
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) plt.title(self.Data[0].stats.station)
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], p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r')
'r') p14, = plt.plot(self.Tcf[islope], dataslope, 'g--')
p14, = plt.plot(self.Tcf[islope], dataslope, 'g--') p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2)
p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2) plt.legend([p11, p12, p13, p14, p15], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'], \
plt.legend([p11, p12, p13, p14, p15], loc='best')
['Data', 'Noise Window', 'Signal Window', plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \
'Slope Window', 'Slope'], self.SNR, self.slope))
loc='best') plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % ( plt.ylabel('Counts')
self.Data[0].stats.station, ax = plt.gca()
self.SNR, self.slope)) plt.yticks([])
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) ax.set_xlim([self.Tcf[inoise[0][0]] - 5, self.Tcf[isignal[0][len(isignal) - 1]] + 5])
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):
''' '''
@ -301,95 +287,90 @@ 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 ' \ print 'PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...'
'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[ cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth
ii1]) / ismooth else:
else: cfsmooth[i] = np.mean(self.cf[1 : i])
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 >= ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \
(self.getpick1() - self.PickWindow / 2)) and & (self.Tcf <= self.getpick1() + self.PickWindow / 2))
(self.Tcf <= cfipick = self.cf[ipick] - np.mean(self.cf[ipick])
(self.getpick1() + self.PickWindow / 2))) Tcfpick = self.Tcf[ipick]
cfipick = self.cf[ipick] - np.mean(self.cf[ipick]) cfsmoothipick = cfsmooth[ipick]- np.mean(self.cf[ipick])
Tcfpick = self.Tcf[ipick] ipick1 = np.argmin(abs(self.Tcf - self.getpick1()))
cfsmoothipick = cfsmooth[ipick] - np.mean(self.cf[ipick]) cfpick1 = 2 * self.cf[ipick1]
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)), for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])):
min([ipick1 + lpickwindow + 1, len(self.cf) - 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_r = self.Tcf[i]
pick_r = self.Tcf[i] self.Pick = pick_r
self.Pick = pick_r flagpick_l = 1
flagpick_l = 1 cfpick_r = self.cf[i]
cfpick_r = self.cf[i] break
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], p3, = plt.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', linewidth=2)
[min(cfipick), max(cfipick)], 'b', linewidth=2) plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick'])
plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick']) plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) plt.yticks([])
plt.yticks([]) plt.title(self.Data[0].stats.station)
plt.title(self.Data[0].stats.station) plt.show()
plt.show() raw_input()
raw_input() plt.close(p)
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