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

Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

Browse Source
This commit is contained in:
Marcel Paffrath 2015-09-28 12:22:09 +02:00
commit 48c889129a
9 changed files with 323 additions and 219 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
QtPyLoT.py
View File

@ -386,8 +386,8 @@ class MainWindow(QMainWindow):
else: else:
raise DatastructureError('not specified') raise DatastructureError('not specified')
return self.fnames return self.fnames
except DatastructureError, e: except DatastructureError as e:
print e print(e)
props = PropertiesDlg(self) props = PropertiesDlg(self)
if props.exec_() == QDialog.Accepted: if props.exec_() == QDialog.Accepted:
return self.getWFFnames() return self.getWFFnames()
@ -410,7 +410,7 @@ class MainWindow(QMainWindow):
def saveData(self): def saveData(self):
def getSavePath(e): def getSavePath(e):
print 'warning: {0}'.format(e) print('warning: {0}'.format(e))
directory = os.path.join(self.getRoot(), self.getEventFileName()) directory = os.path.join(self.getRoot(), self.getEventFileName())
file_filter = "QuakeML file (*.xml);;VELEST observation file format (*.cnv);;NonLinLoc observation file (*.obs)" file_filter = "QuakeML file (*.xml);;VELEST observation file format (*.cnv);;NonLinLoc observation file (*.obs)"
fname = QFileDialog.getSaveFileName(self, 'Save event data ...', fname = QFileDialog.getSaveFileName(self, 'Save event data ...',
@ -513,7 +513,7 @@ class MainWindow(QMainWindow):
def plotWaveformData(self): def plotWaveformData(self):
zne_text = {'Z': 'vertical', 'N': 'north-south', 'E': 'east-west'} zne_text = {'Z': 'vertical', 'N': 'north-south', 'E': 'east-west'}
comp = self.getComponent() comp = self.getComponent()
title = 'overview: {0} components'.format(zne_text[comp]) title = 'section: {0} components'.format(zne_text[comp])
wfst = self.getData().getWFData().select(component=comp) wfst = self.getData().getWFData().select(component=comp)
self.getPlotWidget().plotWFData(wfdata=wfst, title=title) self.getPlotWidget().plotWFData(wfdata=wfst, title=title)
self.draw() self.draw()
@ -574,8 +574,8 @@ class MainWindow(QMainWindow):
def getFilterOptions(self): def getFilterOptions(self):
try: try:
return self.filteroptions[self.getSeismicPhase()] return self.filteroptions[self.getSeismicPhase()]
except AttributeError, e: except AttributeError as e:
print e print(e)
return FilterOptions(None, None, None) return FilterOptions(None, None, None)
def getFilters(self): def getFilters(self):
@ -592,12 +592,12 @@ class MainWindow(QMainWindow):
settings = QSettings() settings = QSettings()
if settings.value("filterdefaults", if settings.value("filterdefaults",
None) is None and not self.getFilters(): None) is None and not self.getFilters():
for key, value in FILTERDEFAULTS.iteritems(): for key, value in FILTERDEFAULTS.items():
self.setFilterOptions(FilterOptions(**value), key) self.setFilterOptions(FilterOptions(**value), key)
elif settings.value("filterdefaults", None) is not None: elif settings.value("filterdefaults", None) is not None:
for key, value in settings.value("filterdefaults"): for key, value in settings.value("filterdefaults"):
self.setFilterOptions(FilterOptions(**value), key) self.setFilterOptions(FilterOptions(**value), key)
except Exception, e: except Exception as e:
self.updateStatus('Error ...') self.updateStatus('Error ...')
emsg = QErrorMessage(self) emsg = QErrorMessage(self)
emsg.showMessage('Error: {0}'.format(e)) emsg.showMessage('Error: {0}'.format(e))
@ -636,8 +636,12 @@ class MainWindow(QMainWindow):
if pickDlg.exec_(): if pickDlg.exec_():
self.setDirty(True) self.setDirty(True)
self.updateStatus('picks accepted ({0})'.format(station)) self.updateStatus('picks accepted ({0})'.format(station))
self.addPicks(station, pickDlg.getPicks()) replot = self.addPicks(station, pickDlg.getPicks())
self.drawPicks(station) if replot:
self.plotWaveformData()
self.drawPicks()
else:
self.drawPicks(station)
else: else:
self.updateStatus('picks discarded ({0})'.format(station)) self.updateStatus('picks discarded ({0})'.format(station))
@ -667,6 +671,7 @@ class MainWindow(QMainWindow):
def addPicks(self, station, picks): def addPicks(self, station, picks):
stat_picks = self.getPicksOnStation(station) stat_picks = self.getPicksOnStation(station)
rval = False
if not stat_picks: if not stat_picks:
stat_picks = picks stat_picks = picks
else: else:
@ -684,11 +689,13 @@ class MainWindow(QMainWindow):
ret = msgBox.exec_() ret = msgBox.exec_()
if ret == QMessageBox.Save: if ret == QMessageBox.Save:
stat_picks = picks stat_picks = picks
rval = True
elif ret == QMessageBox.Cancel: elif ret == QMessageBox.Cancel:
pass pass
else: else:
raise Exception('FATAL: Should never occur!') raise Exception('FATAL: Should never occur!')
self.getPicks()[station] = stat_picks self.getPicks()[station] = stat_picks
return rval
def updatePicks(self): def updatePicks(self):
evt = self.getData().getEvtData() evt = self.getData().getEvtData()

1
icons.qrc
View File

@ -3,6 +3,7 @@
<file>icons/pylot.ico</file> <file>icons/pylot.ico</file>
<file>icons/pylot.png</file> <file>icons/pylot.png</file>
<file>icons/printer.png</file> <file>icons/printer.png</file>
<file>icons/delete.png</file>
<file>icons/key_E.png</file> <file>icons/key_E.png</file>
<file>icons/key_N.png</file> <file>icons/key_N.png</file>
<file>icons/key_P.png</file> <file>icons/key_P.png</file>

BIN
icons/delete.png Executable file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 5.0 KiB

8
icons_rc.py
View File

File diff suppressed because one or more lines are too long

58
pylot/core/analysis/magnitude.py
View File

@ -25,7 +25,8 @@ class Magnitude(object):
:type: float :type: float
:param: pwin, pick window [To To+pwin] to get maximum :param: pwin, pick window [To To+pwin] to get maximum
peak-to-peak amplitude peak-to-peak amplitude (WApp) or to calculate
source spectrum (DCfc)
:type: float :type: float
:param: iplot, no. of figure window for plotting interims results :param: iplot, no. of figure window for plotting interims results
@ -40,6 +41,7 @@ class Magnitude(object):
self.setpwin(pwin) self.setpwin(pwin)
self.setiplot(iplot) self.setiplot(iplot)
self.calcwapp() self.calcwapp()
self.calcsourcespec()
def getwfstream(self): def getwfstream(self):
@ -68,18 +70,23 @@ class Magnitude(object):
def getwapp(self): def getwapp(self):
return self.wapp return self.wapp
def getw0(self):
return self.w0
def getfc(self):
return self.fc
def calcwapp(self): def calcwapp(self):
self.wapp = None self.wapp = None
def calcsourcespec(self): def calcsourcespec(self):
self.sourcespek = None self.sourcespek = None
class WApp(Magnitude): class WApp(Magnitude):
''' '''
Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- Method to derive peak-to-peak amplitude as seen on a Wood-Anderson-
seismograph. Has to be derived from corrected traces! seismograph. Has to be derived from instrument corrected traces!
''' '''
def calcwapp(self): def calcwapp(self):
@ -110,6 +117,7 @@ class WApp(Magnitude):
iwin = getsignalwin(th, self.getTo(), self.getpwin()) iwin = getsignalwin(th, self.getTo(), self.getpwin())
self.wapp = np.max(sqH[iwin]) self.wapp = np.max(sqH[iwin])
print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp
if self.getiplot() > 1: if self.getiplot() > 1:
stream.plot() stream.plot()
f = plt.figure(2) f = plt.figure(2)
@ -128,10 +136,50 @@ class WApp(Magnitude):
class DCfc(Magnitude): class DCfc(Magnitude):
''' '''
Method to calculate the source spectrum and to derive from that the plateau Method to calculate the source spectrum and to derive from that the plateau
(the so-called DC-value) and the corner frequency assuming Aki's omega-square (so-called DC-value) and the corner frequency assuming Aki's omega-square
source model. Has to be derived from corrected traces! source model. Has to be derived from instrument corrected displacement traces!
''' '''
def calcsourcespec(self): def calcsourcespec(self):
print ("Calculating source spectrum ....") print ("Calculating source spectrum ....")
self.w0 = None # DC-value
self.fc = None # corner frequency
stream = self.getwfstream()
tr = stream[0]
# get time array
t = np.arange(0, len(tr) * tr.stats.delta, tr.stats.delta)
iwin = getsignalwin(t, self.getTo(), self.getpwin())
xdat = tr.data[iwin]
# fft
fny = tr.stats.sampling_rate / 2
l = len(xdat) / tr.stats.sampling_rate
n = tr.stats.sampling_rate * l # number of fft bins after Bath
# find next power of 2 of data length
m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
N = int(np.power(m, 2))
y = tr.stats.delta * np.fft.fft(xdat, N)
Y = abs(y[: N/2])
L = (N - 1) / tr.stats.sampling_rate
f = np.arange(0, fny, 1/L)
if self.getiplot() > 1:
f1 = plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t, np.multiply(tr, 1000), 'k') # show displacement in mm
plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') # show displacement in mm
plt.title('Seismogram and P pulse, station %s' % tr.stats.station)
plt.xlabel('Time since %s' % tr.stats.starttime)
plt.ylabel('Displacement [mm]')
plt.subplot(2,1,2)
plt.semilogy(f, Y.real)
plt.title('Source Spectrum from P Pulse')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [m/Hz]')
plt.show()
raw_input()
plt.close(f1)

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

@ -147,7 +147,7 @@ class AICPicker(AutoPicking):
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
@ -155,7 +155,7 @@ class AICPicker(AutoPicking):
#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))
@ -163,15 +163,15 @@ class AICPicker(AutoPicking):
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[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
@ -180,7 +180,7 @@ class AICPicker(AutoPicking):
#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))
@ -189,104 +189,104 @@ class AICPicker(AutoPicking):
#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], self.TSNR[1]) inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], 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 = 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 <= min([self.Pick + tslope, len(self.Data[0].data)])) \
& (self.Tcf >= self.Pick)) & (self.Tcf >= self.Pick))
# find maximum within slope determination window # find maximum within slope determination window
# 'cause slope should be calculated up to first local minimum only! # 'cause slope should be calculated up to first local minimum only!
imax = np.argmax(self.Data[0].data[islope]) imax = np.argmax(self.Data[0].data[islope])
if imax == 0: if imax == 0:
print 'AICPicker: Maximum for slope determination right at the beginning of the window!' print('AICPicker: Maximum for slope determination right at the beginning of the window!')
print 'Choose longer slope determination window!' print('Choose longer slope determination window!')
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')
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
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)
return return
islope = islope[0][0 :imax] islope = islope[0][0 :imax]
dataslope = self.Data[0].data[islope] dataslope = self.Data[0].data[islope]
# calculate slope as polynomal fit of order 1 # calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope), 1) xslope = np.arange(0, len(dataslope), 1)
P = np.polyfit(xslope, dataslope, 1) P = np.polyfit(xslope, dataslope, 1)
datafit = np.polyval(P, xslope) datafit = np.polyval(P, xslope)
if datafit[0] >= datafit[len(datafit) - 1]: if datafit[0] >= datafit[len(datafit) - 1]:
print 'AICPicker: Negative slope, bad onset skipped!' print('AICPicker: Negative slope, bad onset skipped!')
return 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', linewidth=2) p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2)
plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick']) plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick'])
else: else:
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF']) plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
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)
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], '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], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'],
loc='best') loc='best')
plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \ plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope)) self.SNR, self.slope))
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.ylabel('Counts') plt.ylabel('Counts')
plt.yticks([]) plt.yticks([])
plt.show() plt.show()
raw_input() raw_input()
plt.close(p) plt.close(p)
if self.Pick == None: if self.Pick == None:
print 'AICPicker: Could not find minimum, picking window too short?' print('AICPicker: Could not find minimum, picking window too short?')
class PragPicker(AutoPicking): class PragPicker(AutoPicking):
@ -297,102 +297,102 @@ 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
pickflag = 0 pickflag = 0
#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[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 >= self.getpick1() - self.PickWindow / 2) \ ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \
& (self.Tcf <= self.getpick1() + self.PickWindow / 2)) & (self.Tcf <= self.getpick1() + self.PickWindow / 2))
cfipick = self.cf[ipick] - np.mean(self.cf[ipick]) cfipick = self.cf[ipick] - np.mean(self.cf[ipick])
Tcfpick = self.Tcf[ipick] Tcfpick = self.Tcf[ipick]
cfsmoothipick = cfsmooth[ipick]- np.mean(self.cf[ipick]) cfsmoothipick = cfsmooth[ipick]- np.mean(self.cf[ipick])
ipick1 = np.argmin(abs(self.Tcf - self.getpick1())) ipick1 = np.argmin(abs(self.Tcf - self.getpick1()))
cfpick1 = 2 * self.cf[ipick1] 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
cfpick_r = 0 cfpick_r = 0
cfpick_l = 0 cfpick_l = 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)), 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
pickflag = 1
elif flagpick_l > 0 and flagpick_r > 0 and cfpick_l >= cfpick_r:
self.Pick = pick_r
pickflag = 1
elif flagpick_l == 0 and flagpick_r > 0 and cfpick_l >= cfpick_r:
self.Pick = pick_l self.Pick = pick_l
pickflag = 1 pickflag = 1
else: elif flagpick_l > 0 and flagpick_r > 0 and cfpick_l >= cfpick_r:
print 'PragPicker: Could not find reliable onset!' self.Pick = pick_r
self.Pick = None pickflag = 1
pickflag = 0 elif flagpick_l == 0 and flagpick_r > 0 and cfpick_l >= cfpick_r:
self.Pick = pick_l
pickflag = 1
else:
print('PragPicker: Could not find reliable onset!')
self.Pick = None
pickflag = 0
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')
if pickflag > 0: if pickflag > 0:
p3, = plt.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', linewidth=2) p3, = plt.plot([self.Pick, self.Pick], [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:
print 'PragPicker: No initial onset time given! Check input!' print('PragPicker: No initial onset time given! Check input!')
self.Pick = None self.Pick = None
return return

25
pylot/core/pick/autopick.py
View File

@ -11,12 +11,13 @@ function conglomerate utils.
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from scipy import integrate
from pylot.core.pick.Picker import AICPicker, PragPicker from pylot.core.pick.Picker import AICPicker, PragPicker
from pylot.core.pick.CharFuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf from pylot.core.pick.CharFuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker,\ from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker,\
getSNR, fmpicker, checkPonsets, wadaticheck getSNR, fmpicker, checkPonsets, wadaticheck, crossings_nonzero_all
from pylot.core.read.data import Data from pylot.core.read.data import Data
from pylot.core.analysis.magnitude import WApp from pylot.core.analysis.magnitude import WApp, DCfc
def autopickevent(data, param): def autopickevent(data, param):
stations = [] stations = []
@ -309,6 +310,26 @@ def autopickstation(wfstream, pickparam):
else: else:
FM = 'N' FM = 'N'
##############################################################
# get DC value (w0) and corner frequency (fc) of source spectrum
# from P pulse
# restitute streams
# initialize Data object
data = Data()
[corzdat, restflag] = data.restituteWFData(invdir, zdat)
if restflag == 1:
# integrate to displacement
corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta)
# class needs stream object => build it
z_copy = zdat.copy()
z_copy[0].data = corintzdat
# calculate source spectrum and get w0 and fc
calcwin = 1 / bpz2[0] # largest detectable period == window length
# around P pulse for calculating source spectrum
specpara = DCfc(z_copy, mpickP, calcwin, iplot)
w0 = specpara.getw0()
fc = specpara.getfc()
print 'autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \ print 'autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, ' \
'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM) 'Polarity: %s' % (Pweight, SNRP, SNRPdB, FM)
Sflag = 1 Sflag = 1

18
pylot/core/read/inputs.py
View File

@ -45,7 +45,7 @@ class AutoPickParameter(object):
self.__filename = fnin self.__filename = fnin
parFileCont = {} parFileCont = {}
# read from parsed arguments alternatively # read from parsed arguments alternatively
for key, val in kwargs.iteritems(): for key, val in kwargs.items():
parFileCont[key] = val parFileCont[key] = val
if self.__filename is not None: if self.__filename is not None:
@ -57,7 +57,7 @@ class AutoPickParameter(object):
for line in lines: for line in lines:
parspl = line.split('\t')[:2] parspl = line.split('\t')[:2]
parFileCont[parspl[0].strip()] = parspl[1] parFileCont[parspl[0].strip()] = parspl[1]
except Exception, e: except Exception as e:
self._printParameterError(e) self._printParameterError(e)
inputFile.seek(0) inputFile.seek(0)
lines = inputFile.readlines() lines = inputFile.readlines()
@ -65,7 +65,7 @@ class AutoPickParameter(object):
if not line.startswith(('#', '%', '\n', ' ')): if not line.startswith(('#', '%', '\n', ' ')):
parspl = line.split('#')[:2] parspl = line.split('#')[:2]
parFileCont[parspl[1].strip()] = parspl[0].strip() parFileCont[parspl[1].strip()] = parspl[0].strip()
for key, value in parFileCont.iteritems(): for key, value in parFileCont.items():
try: try:
val = int(value) val = int(value)
except: except:
@ -121,7 +121,7 @@ class AutoPickParameter(object):
return len(self.__parameter.keys()) return len(self.__parameter.keys())
def iteritems(self): def iteritems(self):
for key, value in self.__parameter.iteritems(): for key, value in self.__parameter.items():
yield key, value yield key, value
def hasParam(self, parameter): def hasParam(self, parameter):
@ -134,22 +134,22 @@ class AutoPickParameter(object):
for param in args: for param in args:
try: try:
return self.__getitem__(param) return self.__getitem__(param)
except KeyError, e: except KeyError as e:
self._printParameterError(e) self._printParameterError(e)
except TypeError: except TypeError:
try: try:
return self.__getitem__(args) return self.__getitem__(args)
except KeyError, e: except KeyError as e:
self._printParameterError(e) self._printParameterError(e)
def setParam(self, **kwargs): def setParam(self, **kwargs):
for param, value in kwargs.iteritems(): for param, value in kwargs.items():
self.__setitem__(param, value) self.__setitem__(param, value)
print self print(self)
@staticmethod @staticmethod
def _printParameterError(errmsg): def _printParameterError(errmsg):
print 'ParameterError:\n non-existent parameter %s' % errmsg print('ParameterError:\n non-existent parameter %s' % errmsg)
def export2File(self, fnout): def export2File(self, fnout):
fid_out = open(fnout, 'w') fid_out = open(fnout, 'w')

39
pylot/core/util/widgets.py
View File

@ -187,7 +187,7 @@ class PickDlg(QDialog):
try: try:
data = parent.getData().getWFData().copy() data = parent.getData().getWFData().copy()
self.data = data.select(station=station) self.data = data.select(station=station)
except AttributeError, e: except AttributeError as e:
errmsg = 'You either have to put in a data or an appropriate ' \ errmsg = 'You either have to put in a data or an appropriate ' \
'parent (PyLoT MainWindow) object: {0}'.format(e) 'parent (PyLoT MainWindow) object: {0}'.format(e)
raise Exception(errmsg) raise Exception(errmsg)
@ -239,6 +239,8 @@ class PickDlg(QDialog):
zoom_icon.addPixmap(QPixmap(':/icons/zoom_in.png')) zoom_icon.addPixmap(QPixmap(':/icons/zoom_in.png'))
home_icon = QIcon() home_icon = QIcon()
home_icon.addPixmap(QPixmap(':/icons/zoom_0.png')) home_icon.addPixmap(QPixmap(':/icons/zoom_0.png'))
del_icon = QIcon()
del_icon.addPixmap(QPixmap(':/icons/delete.png'))
# create actions # create actions
self.filterAction = createAction(parent=self, text='Filter', self.filterAction = createAction(parent=self, text='Filter',
@ -251,9 +253,12 @@ class PickDlg(QDialog):
slot=self.zoom, icon=zoom_icon, slot=self.zoom, icon=zoom_icon,
tip='Zoom into waveform', tip='Zoom into waveform',
checkable=True) checkable=True)
self.resetAction = createAction(parent=self, text='Home', self.resetZoomAction = createAction(parent=self, text='Home',
slot=self.resetZoom, icon=home_icon, slot=self.resetZoom, icon=home_icon,
tip='Reset zoom to original limits') tip='Reset zoom to original limits')
self.resetPicksAction = createAction(parent=self, text='Delete Picks',
slot=self.delPicks, icon=del_icon,
tip='Delete current picks.')
# create other widget elements # create other widget elements
self.selectPhase = QComboBox() self.selectPhase = QComboBox()
@ -269,7 +274,9 @@ class PickDlg(QDialog):
_dialtoolbar.addWidget(self.selectPhase) _dialtoolbar.addWidget(self.selectPhase)
_dialtoolbar.addAction(self.zoomAction) _dialtoolbar.addAction(self.zoomAction)
_dialtoolbar.addSeparator() _dialtoolbar.addSeparator()
_dialtoolbar.addAction(self.resetAction) _dialtoolbar.addAction(self.resetZoomAction)
_dialtoolbar.addSeparator()
_dialtoolbar.addAction(self.resetPicksAction)
# layout the innermost widget # layout the innermost widget
_innerlayout = QVBoxLayout() _innerlayout = QVBoxLayout()
@ -371,7 +378,7 @@ class PickDlg(QDialog):
traceIDs = [] traceIDs = []
for channel in channels: for channel in channels:
channel = channel.upper() channel = channel.upper()
for traceID, channelID in plotDict.iteritems(): for traceID, channelID in plotDict.items():
if channelID[1].upper().endswith(channel): if channelID[1].upper().endswith(channel):
traceIDs.append(traceID) traceIDs.append(traceID)
return traceIDs return traceIDs
@ -422,6 +429,13 @@ class PickDlg(QDialog):
def getPicks(self): def getPicks(self):
return self.picks return self.picks
def resetPicks(self):
self.picks = {}
def delPicks(self):
self.resetPicks()
self.resetPlot()
def setIniPick(self, gui_event): def setIniPick(self, gui_event):
trace_number = round(gui_event.ydata) trace_number = round(gui_event.ydata)
@ -672,8 +686,21 @@ class PickDlg(QDialog):
zoomx=self.getXLims(), zoomx=self.getXLims(),
zoomy=self.getYLims()) zoomy=self.getYLims())
self.setPlotLabels() self.setPlotLabels()
self.drawPicks()
self.draw() self.draw()
def resetPlot(self):
self.updateCurrentLimits()
data = self.getWFData().copy()
title = self.getPlotWidget().getAxes().get_title()
self.getPlotWidget().plotWFData(wfdata=data, title=title,
zoomx=self.getXLims(),
zoomy=self.getYLims())
self.setPlotLabels()
self.drawPicks()
self.draw()
def setPlotLabels(self): def setPlotLabels(self):
# get channel labels # get channel labels
@ -711,7 +738,7 @@ class PickDlg(QDialog):
else: else:
# deal with something that should never happen # deal with something that should never happen
scale_factor = 1 scale_factor = 1
print gui_event.button print(gui_event.button)
new_xlim = gui_event.xdata - \ new_xlim = gui_event.xdata - \
scale_factor * (gui_event.xdata - self.getXLims()) scale_factor * (gui_event.xdata - self.getXLims())
@ -742,7 +769,7 @@ class PickDlg(QDialog):
def apply(self): def apply(self):
picks = self.getPicks() picks = self.getPicks()
for pick in picks: for pick in picks:
print pick, picks[pick] print(pick, picks[pick])
def accept(self): def accept(self):
self.apply() self.apply()