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

Merge branch 'develop'

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2014-12-18 13:36:05 +01:00
commit 6be85efa1f
6 changed files with 199 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
QtPyLoT.py
View File

@ -143,18 +143,25 @@ class MainWindow(QMainWindow):
filt = """Supported event formats (*.mat *.qml *.xml *.kor
*.evt)"""
caption = 'Select event to open'
fname = QFileDialog().getOpenFileName(self, caption=caption,
self.fname = QFileDialog().getOpenFileName(self, caption=caption,
filter=filt)
else:
fname = unicode(action.data().toString())
self.fname = unicode(action.data().toString())
if not self.okToContinue():
return
else:
return
if fname:
self.data = Data(evtdata=fname)
self.fname = fname
self.data = Data(evtdata=self.fname)
def saveData(self):
settings = QSettings()
exform = settings.value('data/exportFormat', 'None')
try:
self.data.exportEvent(self.fname, exform)
except FormatError:
return False
return True
def getComponent(self):

70
pylot/core/pick/CharFuns.py
View File

@ -117,8 +117,8 @@ class CharacteristicFunction(object):
def getTimeArray(self):
if self.getTime1():
incr = self.getARdetStep()
self.TimeArray = np.arange(0, len(self.getCF()) * incr[0], incr[0]) + self.getCut()[0] \
incr = self.getARdetStep()[0]
self.TimeArray = np.arange(0, len(self.getCF()) * incr, incr) + self.getCut()[0] \
+ self.getTime1() + self.getTime2()
else:
incr = self.getIncrement()
@ -143,18 +143,31 @@ class CharacteristicFunction(object):
cutting window
'''
if cut is not None:
if self.cut[0] == 0:
if len(self.orig_data) == 1:
if self.cut[0] == 0 and self.cut[1] == 0:
start = 0
stop = len(self.orig_data[0])
elif self.cut[0] == 0 and self.cut[1] is not 0:
start = 0
stop = self.cut[1] / self.dt
else:
start = self.cut[0] / self.dt
stop = self.cut[1] / self.dt
if len(self.orig_data) == 1:
zz = self.orig_data.copy()
z1 = zz[0].copy()
zz[0].data = z1.data[start:stop]
data = zz
return data
elif len(self.orig_data) == 2:
if self.cut[0] == 0 and self.cut[1] == 0:
start = 0
stop = min([len(self.orig_data[0]), len(self.orig_data[1])])
elif self.cut[0] == 0 and self.cut[1] is not 0:
start = 0
stop = self.cut[1] / self.dt
else:
start = self.cut[0] / self.dt
stop = self.cut[1] / self.dt
hh = self.orig_data.copy()
h1 = hh[0].copy()
h2 = hh[1].copy()
@ -163,6 +176,15 @@ class CharacteristicFunction(object):
data = hh
return data
elif len(self.orig_data) == 3:
if self.cut[0] == 0 and self.cut[1] == 0:
start = 0
stop = min([len(self.orig_data[0]), len(self.orig_data[1]), len(self.orig_data[2])])
elif self.cut[0] == 0 and self.cut[1] is not 0:
start = 0
stop = self.cut[1] / self.dt
else:
start = self.cut[0] / self.dt
stop = self.cut[1] / self.dt
hh = self.orig_data.copy()
h1 = hh[0].copy()
h2 = hh[1].copy()
@ -195,6 +217,9 @@ class AICcf(CharacteristicFunction):
print 'Calculating AIC ...'
x = self.getDataArray()
xnp = x[0].data
nn = np.isnan(xnp)
if len(nn) > 1:
xnp[nn] = 0
datlen = len(xnp)
k = np.arange(1, datlen)
cf = np.zeros(datlen)
@ -224,6 +249,9 @@ class HOScf(CharacteristicFunction):
x = self.getDataArray(self.getCut())
xnp =x[0].data
nn = np.isnan(xnp)
if len(nn) > 1:
xnp[nn] = 0
if self.getOrder() == 3: # this is skewness
print 'Calculating skewness ...'
y = np.power(xnp, 3)
@ -256,6 +284,9 @@ class HOScf(CharacteristicFunction):
elif self.getOrder() == 4:
LTA[j] = lta / np.power(lta1, 2)
nn = np.isnan(LTA)
if len(nn) > 1:
LTA[nn] = 0
self.cf = LTA
@ -266,6 +297,9 @@ class ARZcf(CharacteristicFunction):
print 'Calculating AR-prediction error from single trace ...'
x = self.getDataArray(self.getCut())
xnp = x[0].data
nn = np.isnan(xnp)
if len(nn) > 1:
xnp[nn] = 0
#some parameters needed
#add noise to time series
xnoise = xnp + np.random.normal(0.0, 1.0, len(xnp)) * self.getFnoise() * max(abs(xnp))
@ -288,6 +322,9 @@ class ARZcf(CharacteristicFunction):
#convert list to numpy array
cf = np.asarray(cf)
nn = np.isnan(cf)
if len(nn) > 1:
cf[nn] = 0
self.cf = cf
@ -376,6 +413,12 @@ class ARHcf(CharacteristicFunction):
print 'Calculating AR-prediction error from both horizontal traces ...'
xnp = self.getDataArray(self.getCut())
n0 = np.isnan(xnp[0].data)
if len(n0) > 1:
xnp[0].data[n0] = 0
n1 = np.isnan(xnp[1].data)
if len(n1) > 1:
xnp[1].data[n1] = 0
#some parameters needed
#add noise to time series
@ -390,7 +433,7 @@ class ARHcf(CharacteristicFunction):
cf = []
loopstep = self.getARdetStep()
for i in range(ldet + self.getOrder() - 3, tend - lpred + 1, loopstep[1]):
for i in range(ldet + self.getOrder() - 1, tend - lpred + 1, loopstep[1]):
self.arDetH(Xnoise, self.getOrder(), i-ldet, i)
#AR prediction of waveform using calculated AR coefficients
self.arPredH(xnp, self.arpara, i + 1, lpred)
@ -401,6 +444,9 @@ class ARHcf(CharacteristicFunction):
#convert list to numpy array
cf = np.asarray(cf)
nn = np.isnan(cf)
if len(nn) > 1:
cf[nn] = 0
self.cf = cf
def arDetH(self, data, order, rind, ldet):
@ -493,6 +539,15 @@ class AR3Ccf(CharacteristicFunction):
print 'Calculating AR-prediction error from all 3 components ...'
xnp = self.getDataArray(self.getCut())
n0 = np.isnan(xnp[0].data)
if len(n0) > 1:
xnp[0].data[n0] = 0
n1 = np.isnan(xnp[1].data)
if len(n1) > 1:
xnp[1].data[n1] = 0
n2 = np.isnan(xnp[2].data)
if len(n2) > 1:
xnp[2].data[n2] = 0
#some parameters needed
#add noise to time series
@ -508,7 +563,7 @@ class AR3Ccf(CharacteristicFunction):
cf = []
loopstep = self.getARdetStep()
for i in range(ldet + self.getOrder() - 3, tend - lpred + 1, loopstep[1]):
for i in range(ldet + self.getOrder() - 1, tend - lpred + 1, loopstep[1]):
self.arDet3C(Xnoise, self.getOrder(), i-ldet, i)
#AR prediction of waveform using calculated AR coefficients
self.arPred3C(xnp, self.arpara, i + 1, lpred)
@ -520,6 +575,9 @@ class AR3Ccf(CharacteristicFunction):
#convert list to numpy array
cf = np.asarray(cf)
nn = np.isnan(cf)
if len(nn) > 1:
cf[nn] = 0
self.cf = cf
def arDet3C(self, data, order, rind, ldet):

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

@ -3,16 +3,16 @@
Created Dec 2014
Implementation of the picking algorithms published and described in:
Küperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group, 2010:
Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group, 2010:
Automated determination of P-phase arrival times at regional and local distances
using higher order statistics, Geophys. J. Int., 181, 1159-1170
Küperkoch, L., Meier, T., Brüstle, 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
autoregressive prediction: application ot local and regional distances, Geophys. J. Int.,
188, 687-702.
:author: MAGS2 EP3 working group / Ludger Küperkoch
:author: MAGS2 EP3 working group / Ludger Kueperkoch
"""
import numpy as np
import matplotlib.pyplot as plt
@ -23,7 +23,7 @@ class AutoPicking(object):
Superclass of different, automated picking algorithms applied on a CF determined
using AIC, HOS, or AR prediction.
'''
def __init__(self, cf, Tslope, aerr, TSNR, PickWindow, peps=None, Tsmooth=None):
def __init__(self, cf, Tslope, aerr, TSNR, PickWindow, aus=None, Tsmooth=None, Pick1=None):
'''
:param: cf, characteristic function, on which the picking algorithm is applied
:type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
@ -41,11 +41,14 @@ class AutoPicking(object):
:param: PickWindow, length of pick window [s]
:type: float
:param: peps, find local minimum at i if aic(i-1)*(1+peps) >= aic(i)
:param: aus ("artificial uplift of samples"), find local minimum at i if aic(i-1)*(1+aus) >= aic(i)
:type: float
:param: Tsmooth, length of moving smoothing window to calculate smoothed CF [s]
:type: float
:param: Pick1, initial (prelimenary) onset time, starting point for PragPicker
:type: float
'''
#assert isinstance(cf, CharFuns), "%s is not a CharacteristicFunction object" % str(cf)
@ -58,8 +61,9 @@ class AutoPicking(object):
self.setaerr(aerr)
self.setTSNR(TSNR)
self.setPickWindow(PickWindow)
self.setpeps(peps)
self.setaus(aus)
self.setTsmooth(Tsmooth)
self.setpick1(Pick1)
self.calcPick()
def __str__(self):
@ -68,15 +72,17 @@ class AutoPicking(object):
aerr:\t{aerr}\n
TSNR:\t\t\t{TSNR}\n
PickWindow:\t{PickWindow}\n
peps:\t{peps}\n
aus:\t{aus}\n
Tsmooth:\t{Tsmooth}\n
Pick1:\t{Pick1}\n
'''.format(name=type(self).__name__,
Tslope=self.getTslope(),
aerr=self.getaerr(),
TSNR=self.getTSNR(),
PickWindow=self.getPickWindow(),
peps=self.getpeps(),
Tsmooth=self.getTsmooth())
aus=self.getaus(),
Tsmooth=self.getTsmooth(),
Pick1=self.getpick1())
def getTslope(self):
return self.Tslope
@ -102,11 +108,11 @@ class AutoPicking(object):
def setPickWindow(self, PickWindow):
self.PickWindow = PickWindow
def getpeps(self):
return self.peps
def getaus(self):
return self.aus
def setpeps(self, peps):
self.peps = peps
def setaus(self, aus):
self.aus = aus
def setTsmooth(self, Tsmooth):
self.Tsmooth = Tsmooth
@ -117,6 +123,12 @@ class AutoPicking(object):
def getpick(self):
return self.Pick
def getpick1(self):
return self.Pick1
def setpick1(self, Pick1):
self.Pick1 = Pick1
def calcPick(self):
self.Pick = None
@ -128,7 +140,7 @@ class AICPicker(AutoPicking):
def calcPick(self):
print 'Get onset (pick) from AIC-CF ...'
print 'Get onset time (pick) from AIC-CF ...'
self.Pick = -1
#taper AIC-CF to get rid off side maxima
@ -155,4 +167,78 @@ class PragPicker(AutoPicking):
def calcPick(self):
print 'Get onset (pick) from HOS- or AR-CF using pragmatic picking algorithm ...'
if self.getpick1() is not None:
print 'Get onset time (pick) from HOS- or AR-CF using pragmatic picking algorithm ...'
self.Pick = -1
#smooth CF
ismooth = round(self.Tsmooth / self.dt);
cfsmooth = np.zeros(len(self.cf))
if len(self.cf) < ismooth:
print 'PragPicker: Tsmooth larger than CF!'
return
else:
for i in range(1, len(self.cf)):
if i > ismooth:
ii1 = i - ismooth;
cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth
else:
cfsmooth[i] = np.mean(self.cf[1 : i])
#select picking window
#which is centered around tpick1
ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \
& (self.Tcf <= self.getpick1() + self.PickWindow / 2))
cfipick = self.cf[ipick]
Tcfpick = self.Tcf[ipick]
cfsmoothipick = cfsmooth[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
#prominent trend: decrease aus
#flat: use given aus
cfdiff = np.diff(cfipick);
i0diff = np.where(cfdiff > 0)
cfdiff = cfdiff[i0diff]
minaus = min(cfdiff * (1 + self.aus));
aus1 = max([minaus, self.aus]);
#at first we look to the right until the end of the pick window is reached
flagpick_r = 0
flagpick_l = 0
flagpick = 0
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])):
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 cfpick1 >= self.cf[i]:
pick_r = self.Tcf[i]
self.Pick = pick_r
flagpick_l = 1
cfpick_r = self.cf[i]
break
#now we look to the left
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 cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
if cfpick1 >= self.cf[i]:
pick_l = self.Tcf[i]
self.Pick = pick_l
flagpick_r = 1
cfpick_l = self.cf[i]
break
#now decide which pick: left or right?
if flagpick_l > 0 and flagpick_r > 0:
if cfpick_l <= cfpick_r:
self.Pick = pick_l
else:
self.Pick = pick_r
else:
self.Pick = -1
print 'PragPicker: No initial onset time given! Check input!'
return

22
pylot/core/pick/run_makeCF.py
View File

@ -19,7 +19,7 @@ def run_makeCF(project, database, event, iplot, station=None):
#parameters for CF calculation
t2 = 7 #length of moving window for HOS calculation [sec]
p = 4 #order of statistics
cuttimes = [10, 40] #start and end time vor CF calculation
cuttimes = [5, 40] #start and end time for CF calculation
bpz = [2, 30] #corner frequencies of bandpass filter, vertical component
bph = [2, 15] #corner frequencies of bandpass filter, horizontal components
tdetz= 1.2 #length of AR-determination window [sec], vertical component
@ -59,9 +59,6 @@ def run_makeCF(project, database, event, iplot, station=None):
#calculate HOS-CF using subclass HOScf of class CharacteristicFunction
hoscf = HOScf(st_copy, cuttimes, t2, p) #instance of HOScf
##############################################################
#get onset time from HOS-CF using class Picker
#hospick = PragPicker(hoscf, 2, 70, [1, 0.5, 0.2], 2, 0.001, 0.2)
##############################################################
#calculate AIC-HOS-CF using subclass AICcf of class CharacteristicFunction
#class needs stream object => build it
tr_aic = tr_filt.copy()
@ -69,9 +66,12 @@ def run_makeCF(project, database, event, iplot, station=None):
st_copy[0].data = tr_aic.data
aiccf = AICcf(st_copy, cuttimes, t2) #instance of AICcf
##############################################################
#get 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, 2, 70, [1, 0.5, 0.2], 3)
##############################################################
#get refined onset time from HOS-CF using class Picker
hospick = PragPicker(hoscf, 2, 70, [1, 0.5, 0.2], 2, 0.001, 0.2, aicpick.getpick())
##############################################################
#calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction
#get stream object of filtered data
st_copy[0].data = tr_filt.data
@ -86,6 +86,9 @@ def run_makeCF(project, database, event, iplot, station=None):
##############################################################
#get onset time from AIC-ARZ-CF using subclass AICPicker of class AutoPicking
aicarzpick = AICPicker(araiccf, 2, 70, [1, 0.5, 0.2], 2)
##############################################################
#get refined onset time from ARZ-CF using class Picker
arzpick = PragPicker(arzcf, 2, 70, [1, 0.5, 0.2], 2, 0, 0.2, aicarzpick.getpick())
elif not wfzfiles:
print 'No vertical component data found!'
@ -156,9 +159,15 @@ def run_makeCF(project, database, event, iplot, station=None):
plt.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1], 'b--')
plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [1, 1], 'b')
plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [-1, -1], 'b')
plt.plot([hospick.getpick(), hospick.getpick()], [-1.3, 1.3], 'r--')
plt.plot([hospick.getpick()-0.5, hospick.getpick()+0.5], [1.3, 1.3], 'r')
plt.plot([hospick.getpick()-0.5, hospick.getpick()+0.5], [-1.3, -1.3], 'r')
plt.plot([aicarzpick.getpick(), aicarzpick.getpick()], [-1.2, 1.2], 'y--')
plt.plot([aicarzpick.getpick()-0.5, aicarzpick.getpick()+0.5], [1.2, 1.2], 'y')
plt.plot([aicarzpick.getpick()-0.5, aicarzpick.getpick()+0.5], [-1.2, -1.2], 'y')
plt.plot([arzpick.getpick(), arzpick.getpick()], [-1.4, 1.4], 'g--')
plt.plot([arzpick.getpick()-0.5, arzpick.getpick()+0.5], [1.4, 1.4], 'g')
plt.plot([arzpick.getpick()-0.5, arzpick.getpick()+0.5], [-1.4, -1.4], 'g')
plt.yticks([])
plt.xlabel('Time [s]')
plt.ylabel('Normalized Counts')
@ -173,7 +182,7 @@ def run_makeCF(project, database, event, iplot, station=None):
th2data = np.arange(0, trH2_filt.stats.npts / trH2_filt.stats.sampling_rate, trH2_filt.stats.delta)
tarhcf = np.arange(0, len(arhcf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdeth +tpredh
p21 = plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k')
p22 = plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r')
p22 = plt.plot(arhcf.getTimeArray(), arhcf.getCF()/max(arhcf.getCF()), 'r')
p23 = plt.plot(arhaiccf.getTimeArray(), arhaiccf.getCF()/max(arhaiccf.getCF()))
plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'b--')
plt.plot([aicarhpick.getpick()-0.5, aicarhpick.getpick()+0.5], [1, 1], 'b')
@ -185,7 +194,6 @@ def run_makeCF(project, database, event, iplot, station=None):
plt.legend([p21, p22, p23], ['Data', 'ARH-CF', 'ARHAIC-CF'])
plt.subplot(212)
plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k')
plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r')
plt.plot(arhaiccf.getTimeArray(), arhaiccf.getCF()/max(arhaiccf.getCF()))
plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'b--')
plt.plot([aicarhpick.getpick()-0.5, aicarhpick.getpick()+0.5], [1, 1], 'b')

4
pylot/core/read/data.py
View File

@ -7,6 +7,7 @@ from obspy.core import (read, Stream)
from obspy import readEvents
from obspy.core.event import (Event, Catalog)
from pylot.core.util import fnConstructor
from pylot.core.util.errors import FormatError
class Data(object):
@ -71,7 +72,8 @@ class Data(object):
from pylot.core.util.defaults import OUTPUTFORMATS
if evtformat.strip() not in OUTPUTFORMATS.values():
evtformat = OUTPUTFORMATS.values()[0]
errmsg = 'selected format {0} not available'.format(evtformat)
raise FormatError(errmsg)
if fnout is None:
ID = self.evtdata.getEventID()

3
pylot/core/util/errors.py
View File

@ -8,3 +8,6 @@ Created on Thu Mar 20 09:47:04 2014
class OptionsError(Exception):
pass
class FormatError(Exception):
pass