Merge branch 'develop' of ariadne.geophysik.rub.de:/data/git/pylot into develop

Conflicts:
	pylot/core/pick/earllatepicker.py
	pylot/core/pick/fmpicker.py
	pylot/core/pick/getSNR.py
This commit is contained in:
Ludger Küperkoch 2015-04-13 09:36:22 +02:00
commit b42b87602b
15 changed files with 690 additions and 324 deletions

View File

@ -29,7 +29,7 @@ from PySide.QtCore import QCoreApplication, QSettings, Signal, QFile, \
QFileInfo, Qt
from PySide.QtGui import QMainWindow, QInputDialog, QIcon, QFileDialog, \
QWidget, QHBoxLayout, QStyle, QKeySequence, QLabel, QFrame, QAction, \
QDialog, QErrorMessage, QApplication
QDialog, QErrorMessage, QApplication, QPixmap
from obspy.core import UTCDateTime
from pylot.core.read import Data, FilterOptions
@ -38,12 +38,11 @@ from pylot.core.util import _getVersionString, FILTERDEFAULTS, fnConstructor, \
NewEventDlg, createEvent, MPLWidget, PropertiesDlg, HelpForm, \
DatastructureError, createAction, getLogin, createCreationInfo, PickDlg
from pylot.core.util.structure import DATASTRUCTURE
import qrc_resources
import icons_rc
# Version information
__version__ = _getVersionString()
class MainWindow(QMainWindow):
closing = Signal()
@ -101,8 +100,11 @@ class MainWindow(QMainWindow):
except:
self.startTime = UTCDateTime()
pylot_icon = QIcon()
pylot_icon.addPixmap(QPixmap(':/icons/pylot.ico'))
self.setWindowTitle("PyLoT - do seismic processing the python way")
self.setWindowIcon(QIcon(":/icon.ico"))
self.setWindowIcon(pylot_icon)
xlab = self.startTime.strftime('seconds since %Y/%m/%d %H:%M:%S (%Z)')
@ -123,7 +125,15 @@ class MainWindow(QMainWindow):
saveIcon = self.style().standardIcon(QStyle.SP_DriveHDIcon)
helpIcon = self.style().standardIcon(QStyle.SP_DialogHelpButton)
newIcon = self.style().standardIcon(QStyle.SP_FileIcon)
pickIcon = QIcon(':/pick.png')
# create resource icons
p_icon = QIcon()
p_icon.addPixmap(QPixmap(':/icons/picon.png'))
s_icon = QIcon()
s_icon.addPixmap(QPixmap(':/icons/sicon.png'))
print_icon = QIcon()
print_icon.addPixmap(QPixmap(':/icons/printer.png'))
newEventAction = self.createAction(self, "&New event ...",
self.createNewEvent,
QKeySequence.New, newIcon,
@ -140,10 +150,6 @@ class MainWindow(QMainWindow):
"Ctrl+W", QIcon(":/wfIcon.png"),
"""Open waveform data (event will
be closed).""")
selectStation = self.createAction(self, "Select station",
self.pickOnStation, "Alt+P", pickIcon,
"Select a station from overview "
"plot for picking")
prefsEventAction = self.createAction(self, "Preferences",
self.PyLoTprefs,
QKeySequence.Preferences,
@ -163,14 +169,14 @@ class MainWindow(QMainWindow):
"Alt+F", QIcon(None),
"""Adjust filter parameters.""")
self.selectPAction = self.createAction(self, "&P", self.alterPhase, "Alt+P",
QIcon(":/picon.png"),
p_icon,
"Toggle P phase.", True)
self.selectSAction = self.createAction(self, "&S", self.alterPhase, "Alt+S",
QIcon(":/sicon.png"),
s_icon,
"Toggle S phase", True)
printAction = self.createAction(self, "&Print event ...",
self.printEvent, QKeySequence.Print,
QIcon(":/printer.png"),
print_icon,
"Print waveform overview.")
helpAction = self.createAction(self, "&Help ...", self.helpHelp,
QKeySequence.HelpContents, helpIcon,
@ -204,10 +210,10 @@ class MainWindow(QMainWindow):
phaseToolBar.setObjectName("PhaseTools")
self.addActions(phaseToolBar, phaseToolActions)
pickToolBar = self.addToolBar("PickTools")
pickToolActions = (selectStation, )
pickToolBar.setObjectName("PickTools")
self.addActions(pickToolBar, pickToolActions)
# pickToolBar = self.addToolBar("PickTools")
# pickToolActions = (selectStation, )
# pickToolBar.setObjectName("PickTools")
# self.addActions(pickToolBar, pickToolActions)
self.eventLabel = QLabel()
self.eventLabel.setFrameStyle(QFrame.StyledPanel | QFrame.Sunken)
@ -336,9 +342,9 @@ class MainWindow(QMainWindow):
def getPlotWidget(self):
return self.DataPlot
def getWFID(self, event):
def getWFID(self, gui_event):
ycoord = event.ydata
ycoord = gui_event.ydata
statID = int(round(ycoord))
@ -370,6 +376,10 @@ class MainWindow(QMainWindow):
title = 'overview: {0} components'.format(zne_text[comp])
wfst = self.getData().getWFData().select(component=comp)
self.getPlotWidget().plotWFData(wfdata=wfst, title=title)
self.getPlotWidget().draw()
pos = self.getPlotWidget().getPlotDict().keys()
labels = [int(act) for act in pos]
self.getPlotWidget().setYTickLabels(pos, labels)
def filterWaveformData(self):
if self.getData():
@ -444,7 +454,7 @@ class MainWindow(QMainWindow):
return self.seismicPhase
def getStationName(self, wfID):
return self.getPlotWidget().getPlotDict()[wfID]
return self.getPlotWidget().getPlotDict()[wfID][0]
def alterPhase(self):
pass
@ -454,9 +464,9 @@ class MainWindow(QMainWindow):
self.updateStatus('Seismic phase changed to '
'{0}'.format(self.getSeismicPhase()))
def pickOnStation(self, event):
def pickOnStation(self, gui_event):
wfID = self.getWFID(event)
wfID = self.getWFID(gui_event)
station = self.getStationName(wfID)
print 'picking on station {0}'.format(station)
@ -517,13 +527,16 @@ class MainWindow(QMainWindow):
def main():
# create the Qt application
pylot_app = QApplication(sys.argv[0])
pylot_app = QApplication(sys.argv)
app_icon = QIcon()
app_icon.addPixmap(QPixmap(':/icons/pick.png'))
# set Application Information
pylot_app.setOrganizationName("Ruhr-University Bochum / MAGS2")
pylot_app.setOrganizationDomain("rub.de")
pylot_app.setApplicationName("PyLoT")
pylot_app.setWindowIcon(QIcon(":/icon.ico"))
pylot_app.setWindowIcon(app_icon)
# create the main window
pylot_form = MainWindow()

View File

@ -5,6 +5,8 @@
<file>icons/picon.png</file>
<file>icons/sicon.png</file>
<file>icons/pick.png</file>
<file>icons/filter.png</file>
<file>icons/zoom.png</file>
</qresource>
<qresource prefix="/help">
<file>help/index.html</file>

BIN
icons/filter.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 4.6 KiB

BIN
icons/zoom.png Executable file

Binary file not shown.

After

Width:  |  Height:  |  Size: 7.5 KiB

BIN
icons/zoom_minus.png Executable file

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.9 KiB

21
icons_rc.py Normal file

File diff suppressed because one or more lines are too long

View File

@ -0,0 +1,27 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Created Mar 2015
Transcription of the rezipe of Diehl et al. (2009) for consistent phase
picking. For a given inital (the most likely) pick, the corresponding earliest
and latest possible pick is calculated based on noise measurements in front of
the most likely pick and signal wavelength derived from zero crossings.
:author: Ludger Kueperkoch / MAGS2 EP3 working group
"""
import argparse
import obspy
from pylot.core.pick.utils import earllatepicker
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--X', type=~obspy.core.stream.Stream, help='time series (seismogram) read with obspy module read')
parser.add_argument('--nfac', type=int, help='(noise factor), nfac times noise level to calculate latest possible pick')
parser.add_argument('--TSNR', type=tuple, help='length of time windows around pick used to determine SNR \
[s] (Tnoise, Tgap, Tsignal)')
parser.add_argument('--Pick1', type=float, help='Onset time of most likely pick')
parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs')
args = parser.parse_args()
earllatepicker(args.X, args.nfac, args.TSNR, args.Pick1, args.iplot)

23
pylot/core/pick/fmpicker.py Executable file
View File

@ -0,0 +1,23 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Created Mar 2015
Function to derive first motion (polarity) for given phase onset based on zero crossings.
:author: MAGS2 EP3 working group / Ludger Kueperkoch
"""
import argparse
import obspy
from pylot.core.pick.utils import fmpicker
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--Xraw', type=~obspy.core.stream.Stream, help='unfiltered time series (seismogram) read with obspy module read')
parser.add_argument('--Xfilt', type=~obspy.core.stream.Stream, help='filtered time series (seismogram) read with obspy module read')
parser.add_argument('--pickwin', type=float, help='length of pick window [s] for first motion determination')
parser.add_argument('--Pick', type=float, help='Onset time of most likely pick')
parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs')
args = parser.parse_args()
fmpicker(args.Xraw, args.Xfilt, args.pickwin, args.Pick, args.iplot)

30
pylot/core/pick/getSNR.py Normal file
View File

@ -0,0 +1,30 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Created Mar/Apr 2015
Function to calculate SNR of certain part of seismogram relative
to given time. Returns SNR and SNR [dB].
:author: Ludger Kueperkoch /MAGS EP3 working group
"""
import argparse
import obspy
from pylot.core.pick.utils import getSNR
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--data', '-d', type=~obspy.core.stream.Stream,
help='time series (seismogram) read with obspy module '
'read',
dest='data')
parser.add_argument('--tsnr', '-s', type=tuple,
help='length of time windows around pick used to '
'determine SNR [s] (Tnoise, Tgap, Tsignal)',
dest='tsnr')
parser.add_argument('--time', '-t', type=float,
help='initial time from which noise and signal windows '
'are calculated',
dest='time')
args = parser.parse_args()
print getSNR(args.data, args.tsnr, args.time)

View File

@ -0,0 +1,15 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import numpy
from pylot.core.pick.utils import getnoisewin
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--t', type=~numpy.array, help='numpy array of time stamps')
parser.add_argument('--t1', type=float, help='time from which relativ to it noise window is extracted')
parser.add_argument('--tnoise', type=float, help='length of time window [s] for noise part extraction')
parser.add_argument('--tgap', type=float, help='safety gap between signal (t1=onset) and noise')
args = parser.parse_args()
getnoisewin(args.t, args.t1, args.tnoise, args.tgap)

View File

@ -0,0 +1,14 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import numpy
from pylot.core.pick.utils import getsignalwin
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--t', type=~numpy.array, help='numpy array of time stamps')
parser.add_argument('--t1', type=float, help='time from which relativ to it signal window is extracted')
parser.add_argument('--tsignal', type=float, help='length of time window [s] for signal part extraction')
args = parser.parse_args()
getsignalwin(args.t, args.t1, args.tsignal)

View File

@ -11,7 +11,7 @@
import numpy as np
import matplotlib.pyplot as plt
from obspy.core import Stream
import argparse
def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
'''
@ -45,8 +45,9 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...'
x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta)
#get latest possible pick
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
X[0].stats.delta)
# get latest possible pick
#get noise window
inoise = getnoisewin(t, Pick1, TSNR[0], TSNR[1])
#get signal window
@ -57,8 +58,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
ilup = np.where(x[isignal] > nlevel)
ildown = np.where(x[isignal] < -nlevel)
if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
print 'earllatepicker: Signal lower than noise level, misspick?'
return
print 'earllatepicker: Signal lower than noise level, misspick?'
return
il = min([ilup[0][0], ildown[0][0]])
LPick = t[isignal][il]
@ -68,62 +69,57 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
zc = []
zc.append(Pick1)
i = 0
for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
i = i+ 1
if x[j-1] <= 0 and x[j] >= 0:
zc.append(t[isignal][i])
elif x[j-1] > 0 and x[j] <= 0:
zc.append(t[isignal][i])
for j in range(isignal[0][1], isignal[0][len(t[isignal]) - 1]):
i = i + 1
if x[j - 1] <= 0 and x[j] >= 0:
zc.append(t[isignal][i])
elif x[j - 1] > 0 and x[j] <= 0:
zc.append(t[isignal][i])
if len(zc) == 3:
break
#calculate maximum period T0 of signal out of zero crossings
T0 = max(np.diff(zc)) #this is half wave length!
T0 = max(np.diff(zc)) #this is half wave length!
#T0/4 is assumed as time difference between most likely and earliest possible pick!
EPick = Pick1 - T0/2
EPick = Pick1 - T0 / 2
#get symmetric pick error as mean from earliest and latest possible pick
#by weighting latest possible pick two times earliest possible pick
diffti_tl = LPick - Pick1
diffti_te = Pick1 - EPick
PickError = (diffti_te + 2 * diffti_tl) / 3
PickError = (diffti_te + 2 * diffti_tl) / 3
if iplot is not None:
plt.figure(iplot)
p1, = plt.plot(t, x, 'k')
p2, = plt.plot(t[inoise], x[inoise])
p3, = plt.plot(t[isignal], x[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
p5, = plt.plot(zc, [0, 0, 0], '*g', markersize=14)
plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
plt.figure(iplot)
p1, = plt.plot(t, x, 'k')
p2, = plt.plot(t[inoise], x[inoise])
p3, = plt.plot(t[isignal], x[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
p5, = plt.plot(zc, [0, 0, 0], '*g', markersize=14)
plt.legend([p1, p2, p3, p4, p5],
['Data', 'Noise Window', 'Signal Window', 'Noise Level',
'Zero Crossings'], \
loc='best')
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2)
plt.plot([LPick, LPick], [max(x)/2, -max(x)/2], '--k')
plt.plot([EPick, EPick], [max(x)/2, -max(x)/2], '--k')
plt.plot([Pick1 + PickError, Pick1 + PickError], [max(x)/2, -max(x)/2], 'r--')
plt.plot([Pick1 - PickError, Pick1 - PickError], [max(x)/2, -max(x)/2], 'r--')
plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
plt.yticks([])
ax = plt.gca()
ax.set_xlim([t[inoise[0][0]] - 2, t[isignal[0][len(isignal) - 1]] + 3])
plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % X[0].stats.station)
plt.show()
raw_input()
plt.close(iplot)
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2)
plt.plot([LPick, LPick], [max(x) / 2, -max(x) / 2], '--k')
plt.plot([EPick, EPick], [max(x) / 2, -max(x) / 2], '--k')
plt.plot([Pick1 + PickError, Pick1 + PickError],
[max(x) / 2, -max(x) / 2], 'r--')
plt.plot([Pick1 - PickError, Pick1 - PickError],
[max(x) / 2, -max(x) / 2], 'r--')
plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
plt.yticks([])
ax = plt.gca()
ax.set_xlim([t[inoise[0][0]] - 2, t[isignal[0][len(isignal) - 1]] + 3])
plt.title(
'Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' %
X[0].stats.station)
plt.show()
raw_input()
plt.close(iplot)
return EPick, LPick, PickError
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--X', type=~obspy.core.stream.Stream, help='time series (seismogram) read with obspy module read')
parser.add_argument('--nfac', type=int, help='(noise factor), nfac times noise level to calculate latest possible pick')
parser.add_argument('--TSNR', type=tuple, help='length of time windows around pick used to determine SNR \
[s] (Tnoise, Tgap, Tsignal)')
parser.add_argument('--Pick1', type=float, help='Onset time of most likely pick')
parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs')
args = parser.parse_args()
earllatepicker(args.X, args.nfac, args.TSNR, args.Pick1, args.iplot)
def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
'''
@ -152,149 +148,146 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
FM = None
if Pick is not None:
print 'fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...'
print 'fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...'
xraw = Xraw[0].data
xfilt = Xfilt[0].data
t = np.arange(0, Xraw[0].stats.npts / Xraw[0].stats.sampling_rate, Xraw[0].stats.delta)
#get pick window
ipick = np.where((t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick))
#remove mean
xraw[ipick] = xraw[ipick] - np.mean(xraw[ipick])
xfilt[ipick] = xfilt[ipick] - np.mean(xfilt[ipick])
xraw = Xraw[0].data
xfilt = Xfilt[0].data
t = np.arange(0, Xraw[0].stats.npts / Xraw[0].stats.sampling_rate,
Xraw[0].stats.delta)
# get pick window
ipick = np.where(
(t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick))
#remove mean
xraw[ipick] = xraw[ipick] - np.mean(xraw[ipick])
xfilt[ipick] = xfilt[ipick] - np.mean(xfilt[ipick])
#get next zero crossing after most likely pick
#initial onset is assumed to be the first zero crossing
#first from unfiltered trace
zc1 = []
zc1.append(Pick)
index1 = []
i = 0
for j in range(ipick[0][1],ipick[0][len(t[ipick]) - 1]):
i = i+ 1
if xraw[j-1] <= 0 and xraw[j] >= 0:
zc1.append(t[ipick][i])
index1.append(i)
elif xraw[j-1] > 0 and xraw[j] <= 0:
zc1.append(t[ipick][i])
index1.append(i)
if len(zc1) == 3:
break
#get next zero crossing after most likely pick
#initial onset is assumed to be the first zero crossing
#first from unfiltered trace
zc1 = []
zc1.append(Pick)
index1 = []
i = 0
for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]):
i = i + 1
if xraw[j - 1] <= 0 and xraw[j] >= 0:
zc1.append(t[ipick][i])
index1.append(i)
elif xraw[j - 1] > 0 and xraw[j] <= 0:
zc1.append(t[ipick][i])
index1.append(i)
if len(zc1) == 3:
break
#if time difference betweeen 1st and 2cnd zero crossing
#is too short, get time difference between 1st and 3rd
#to derive maximum
if zc1[1] - zc1[0] <= Xraw[0].stats.delta:
li1 = index1[1]
else:
li1 = index1[0]
if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 0:
print 'earllatepicker: Onset on unfiltered trace too emergent for first motion determination!'
P1 = None
else:
imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]]))
islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1]))
#calculate slope as polynomal fit of order 1
xslope1 = np.arange(0, len(xraw[islope1]), 1)
P1 = np.polyfit(xslope1, xraw[islope1], 1)
datafit1 = np.polyval(P1, xslope1)
#if time difference betweeen 1st and 2cnd zero crossing
#is too short, get time difference between 1st and 3rd
#to derive maximum
if zc1[1] - zc1[0] <= Xraw[0].stats.delta:
li1 = index1[1]
else:
li1 = index1[0]
if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 0:
print 'earllatepicker: Onset on unfiltered trace too emergent for first motion determination!'
P1 = None
else:
imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]]))
islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1]))
#calculate slope as polynomal fit of order 1
xslope1 = np.arange(0, len(xraw[islope1]), 1)
P1 = np.polyfit(xslope1, xraw[islope1], 1)
datafit1 = np.polyval(P1, xslope1)
#now using filterd trace
#next zero crossing after most likely pick
zc2 = []
zc2.append(Pick)
index2 = []
i = 0
for j in range(ipick[0][1],ipick[0][len(t[ipick]) - 1]):
i = i+ 1
if xfilt[j-1] <= 0 and xfilt[j] >= 0:
zc2.append(t[ipick][i])
index2.append(i)
elif xfilt[j-1] > 0 and xfilt[j] <= 0:
zc2.append(t[ipick][i])
index2.append(i)
if len(zc2) == 3:
break
#now using filterd trace
#next zero crossing after most likely pick
zc2 = []
zc2.append(Pick)
index2 = []
i = 0
for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]):
i = i + 1
if xfilt[j - 1] <= 0 and xfilt[j] >= 0:
zc2.append(t[ipick][i])
index2.append(i)
elif xfilt[j - 1] > 0 and xfilt[j] <= 0:
zc2.append(t[ipick][i])
index2.append(i)
if len(zc2) == 3:
break
#if time difference betweeen 1st and 2cnd zero crossing
#is too short, get time difference between 1st and 3rd
#to derive maximum
if zc2[1] - zc2[0] <= Xfilt[0].stats.delta:
li2 = index2[1]
else:
li2 = index2[0]
if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 0:
print 'earllatepicker: Onset on filtered trace too emergent for first motion determination!'
P2 = None
else:
imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]]))
islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2]))
#calculate slope as polynomal fit of order 1
xslope2 = np.arange(0, len(xfilt[islope2]), 1)
P2 = np.polyfit(xslope2, xfilt[islope2], 1)
datafit2 = np.polyval(P2, xslope2)
#if time difference betweeen 1st and 2cnd zero crossing
#is too short, get time difference between 1st and 3rd
#to derive maximum
if zc2[1] - zc2[0] <= Xfilt[0].stats.delta:
li2 = index2[1]
else:
li2 = index2[0]
if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 0:
print 'earllatepicker: Onset on filtered trace too emergent for first motion determination!'
P2 = None
else:
imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]]))
islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2]))
#calculate slope as polynomal fit of order 1
xslope2 = np.arange(0, len(xfilt[islope2]), 1)
P2 = np.polyfit(xslope2, xfilt[islope2], 1)
datafit2 = np.polyval(P2, xslope2)
#compare results
if P1 is not None and P2 is not None:
if P1[0] < 0 and P2[0] < 0:
FM = 'D'
elif P1[0] >= 0 and P2[0] < 0:
FM = '-'
elif P1[0] < 0 and P2[0]>= 0:
FM = '-'
elif P1[0] > 0 and P2[0] > 0:
FM = 'U'
elif P1[0] <= 0 and P2[0] > 0:
FM = '+'
elif P1[0] > 0 and P2[0] <= 0:
FM = '+'
#compare results
if P1 is not None and P2 is not None:
if P1[0] < 0 and P2[0] < 0:
FM = 'D'
elif P1[0] >= 0 and P2[0] < 0:
FM = '-'
elif P1[0] < 0 and P2[0] >= 0:
FM = '-'
elif P1[0] > 0 and P2[0] > 0:
FM = 'U'
elif P1[0] <= 0 and P2[0] > 0:
FM = '+'
elif P1[0] > 0 and P2[0] <= 0:
FM = '+'
if iplot is not None:
plt.figure(iplot)
plt.subplot(2,1,1)
plt.plot(t, xraw, 'k')
p1, = plt.plot([Pick, Pick], [max(xraw), -max(xraw)], 'b', linewidth=2)
if P1 is not None:
p2, = plt.plot(t[islope1], xraw[islope1])
p3, = plt.plot(zc1, np.zeros(len(zc1)), '*g', markersize=14)
p4, = plt.plot(t[islope1], datafit1, '--g', linewidth=2)
plt.legend([p1, p2, p3, p4], ['Pick', 'Slope Window', 'Zero Crossings', 'Slope'], \
loc='best')
plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax = plt.gca()
ax.set_xlim([t[islope1[0][0]] - 0.1, t[islope1[0][len(islope1) - 1]] + 0.3])
plt.yticks([])
plt.title('First-Motion Determination, %s, Unfiltered Data' % Xraw[0].stats.station)
plt.figure(iplot)
plt.subplot(2, 1, 1)
plt.plot(t, xraw, 'k')
p1, = plt.plot([Pick, Pick], [max(xraw), -max(xraw)], 'b', linewidth=2)
if P1 is not None:
p2, = plt.plot(t[islope1], xraw[islope1])
p3, = plt.plot(zc1, np.zeros(len(zc1)), '*g', markersize=14)
p4, = plt.plot(t[islope1], datafit1, '--g', linewidth=2)
plt.legend([p1, p2, p3, p4],
['Pick', 'Slope Window', 'Zero Crossings', 'Slope'], \
loc='best')
plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax = plt.gca()
ax.set_xlim(
[t[islope1[0][0]] - 0.1, t[islope1[0][len(islope1) - 1]] + 0.3])
plt.yticks([])
plt.title('First-Motion Determination, %s, Unfiltered Data' % Xraw[
0].stats.station)
plt.subplot(2,1,2)
plt.title('First-Motion Determination, Filtered Data')
plt.plot(t, xfilt, 'k')
p1, = plt.plot([Pick, Pick], [max(xfilt), -max(xfilt)], 'b', linewidth=2)
if P2 is not None:
p2, = plt.plot(t[islope2], xfilt[islope2])
p3, = plt.plot(zc2, np.zeros(len(zc2)), '*g', markersize=14)
p4, = plt.plot(t[islope2], datafit2, '--g', linewidth=2)
plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax = plt.gca()
ax.set_xlim([t[islope2[0][0]] - 0.1, t[islope2[0][len(islope2) - 1]] + 0.3])
plt.xlabel('Time [s] since %s' % Xraw[0].stats.starttime)
plt.yticks([])
plt.show()
raw_input()
plt.close(iplot)
plt.subplot(2, 1, 2)
plt.title('First-Motion Determination, Filtered Data')
plt.plot(t, xfilt, 'k')
p1, = plt.plot([Pick, Pick], [max(xfilt), -max(xfilt)], 'b',
linewidth=2)
if P2 is not None:
p2, = plt.plot(t[islope2], xfilt[islope2])
p3, = plt.plot(zc2, np.zeros(len(zc2)), '*g', markersize=14)
p4, = plt.plot(t[islope2], datafit2, '--g', linewidth=2)
plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax = plt.gca()
ax.set_xlim(
[t[islope2[0][0]] - 0.1, t[islope2[0][len(islope2) - 1]] + 0.3])
plt.xlabel('Time [s] since %s' % Xraw[0].stats.starttime)
plt.yticks([])
plt.show()
raw_input()
plt.close(iplot)
return FM
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--Xraw', type=~obspy.core.stream.Stream, help='unfiltered time series (seismogram) read with obspy module read')
parser.add_argument('--Xfilt', type=~obspy.core.stream.Stream, help='filtered time series (seismogram) read with obspy module read')
parser.add_argument('--pickwin', type=float, help='length of pick window [s] for first motion determination')
parser.add_argument('--Pick', type=float, help='Onset time of most likely pick')
parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs')
args = parser.parse_args()
earllatepicker(args.Xraw, args.Xfilt, args.pickwin, args.Pick, args.iplot)
def getSNR(X, TSNR, t1):
'''
@ -314,35 +307,28 @@ def getSNR(X, TSNR, t1):
assert isinstance(X, Stream), "%s is not a stream object" % str(X)
SNR = None
SNRdB = None
x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta)
#get noise window
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
X[0].stats.delta)
# get noise window
inoise = getnoisewin(t, t1, TSNR[0], TSNR[1])
#get signal window
isignal = getsignalwin(t, t1, TSNR[2])
if np.size(inoise) < 1:
print 'getSNR: Empty array inoise, check noise window!'
return
print 'getSNR: Empty array inoise, check noise window!'
return
elif np.size(isignal) < 1:
print 'getSNR: Empty array isignal, check signal window!'
return
print 'getSNR: Empty array isignal, check signal window!'
return
#calculate ratios
SNR = max(abs(x[isignal])) / np.mean(abs(x[inoise]))
noiselevel = np.mean(abs(x[inoise]))
SNR = max(abs(x[isignal])) / noiselevel
SNRdB = 20 * np.log10(SNR)
return SNR, SNRdB
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--X', type=~obspy.core.stream.Stream, help='time series (seismogram) read with obspy module read')
parser.add_argument('--TSNR', type=tuple, help='length of time windows around pick used to determine SNR \
[s] (Tnoise, Tgap, Tsignal)')
parser.add_argument('--t1', type=float, help='initial time from which noise and signal windows are calculated')
args = parser.parse_args()
getSNR(args.X, args.TSNR, args.t1)
return SNR, SNRdB, noiselevel
def getnoisewin(t, t1, tnoise, tgap):
@ -365,22 +351,14 @@ def getnoisewin(t, t1, tnoise, tgap):
'''
inoise = None
#get noise window
# get noise window
inoise = np.where((t <= max([t1 - tgap, 0])) \
& (t >= max([t1 - tnoise - tgap, 0])))
& (t >= max([t1 - tnoise - tgap, 0])))
if np.size(inoise) < 1:
print 'getnoisewin: Empty array inoise, check noise window!'
print 'getnoisewin: Empty array inoise, check noise window!'
return inoise
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--t', type=array, help='numpy array of time stamps')
parser.add_argument('--t1', type=float, help='time from which relativ to it noise window is extracted')
parser.add_argument('--tnoise', type=float, help='length of time window [s] for noise part extraction')
parser.add_argument('--tgap', type=float, help='safety gap between signal (t1=onset) and noise')
args = parser.parse_args()
getnoisewin(args.t, args.t1, args.tnoise, args.tgap)
def getsignalwin(t, t1, tsignal):
'''
@ -398,18 +376,12 @@ def getsignalwin(t, t1, tsignal):
'''
inoise = None
#get signal window
# get signal window
isignal = np.where((t <= min([t1 + tsignal, len(t)])) \
& (t >= t1))
& (t >= t1))
if np.size(isignal) < 1:
print 'getsignalwin: Empty array isignal, check signal window!'
print 'getsignalwin: Empty array isignal, check signal window!'
return isignal
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument('--t', type=array, help='numpy array of time stamps')
parser.add_argument('--t1', type=float, help='time from which relativ to it signal window is extracted')
parser.add_argument('--tsignal', type=float, help='length of time window [s] for signal part extraction')
args = parser.parse_args()
getsignalwin(args.t, args.t1, args.tsignal)

View File

@ -14,34 +14,18 @@ matplotlib.rcParams['backend.qt4'] = 'PySide'
from matplotlib.figure import Figure
from matplotlib.backends.backend_qt4agg import FigureCanvas
from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg
from matplotlib.widgets import MultiCursor
from PySide.QtGui import (QAction,
QApplication,
QComboBox,
QDateTimeEdit,
QDialog,
QDialogButtonBox,
QDoubleSpinBox,
QGroupBox,
QGridLayout,
QHBoxLayout,
QIcon,
QKeySequence,
QLabel,
QLineEdit,
QMessageBox,
QSpinBox,
QTabWidget,
QToolBar,
QVBoxLayout,
QWidget)
from PySide.QtCore import (QSettings,
Qt,
QUrl,
Signal,
Slot)
from PySide.QtGui import QAction, QApplication,QComboBox, QDateTimeEdit,\
QDialog, QDialogButtonBox, QDoubleSpinBox, QGroupBox, QGridLayout,\
QIcon, QKeySequence, QLabel, QLineEdit, QMessageBox, QPixmap, QSpinBox,\
QTabWidget, QToolBar, QVBoxLayout, QWidget
from PySide.QtCore import QSettings, Qt, QUrl, Signal, Slot
from PySide.QtWebKit import QWebView
from obspy import Stream, UTCDateTime
from obspy.core.event import Pick, WaveformStreamID
from pylot.core.read import FilterOptions
from pylot.core.pick.utils import getSNR
from pylot.core.util.defaults import OUTPUTFORMATS
from pylot.core.util import prepTimeAxis, getGlobalTimes
@ -71,6 +55,7 @@ class MPLWidget(FigureCanvas):
self._parent = None
self.setParent(parent)
self.figure = Figure()
self.figure.set_facecolor((.92, .92, .92))
# attribute plotdict is an dictionary connecting position and a name
self.plotdict = dict()
# create axes
@ -91,39 +76,58 @@ class MPLWidget(FigureCanvas):
def setPlotDict(self, key, value):
self.plotdict[key] = value
def clearPlotDict(self):
self.plotdict = dict()
def getParent(self):
return self._parent
def setParent(self, parent):
self._parent = parent
def plotWFData(self, wfdata, title = None):
self.axes.lines = []
def plotWFData(self, wfdata, title=None, zoomx=None, zoomy=None):
self.axes.cla()
self.clearPlotDict()
wfstart = getGlobalTimes(wfdata)[0]
for n, trace in enumerate(wfdata):
channel = trace.stats.channel
station = trace.stats.station
print('plotting station: %s' % station)
msg = 'plotting %s channel of station %s' % (channel, station)
print(msg)
stime = trace.stats.starttime - wfstart
time_ax = prepTimeAxis(stime, trace)
trace.detrend()
trace.detrend('demean')
trace.normalize(trace.data.max() * 2)
self.axes.plot(time_ax, trace.data + n, 'k')
self.axes.hold(True)
xlabel = 'seconds since {0}'.format(wfstart)
ylabel = ''
self.updateWidget(xlabel, ylabel, title)
self.setPlotDict(n, station)
self.setPlotDict(n, (station, channel))
self.axes.autoscale(tight=True)
if zoomx:
self.axes.set_xlim(zoomx)
if zoomy:
self.axes.set_ylim(zoomy)
self.draw()
def setYTickLabels(self, pos, labels):
self.axes.set_yticks(pos)
self.axes.set_yticklabels(labels)
self.draw()
def updateXLabel(self, text):
self.axes.set_xlabel(text)
self.draw()
def updateYLabel(self, text):
self.axes.set_ylabel(text)
self.draw()
def updateTitle(self, text):
self.axes.set_title(text)
self.draw()
def updateWidget(self, xlabel, ylabel, title):
self.updateXLabel(xlabel)
@ -237,6 +241,12 @@ class PickDlg(QDialog):
self.station = station
self.rotate = rotate
self.components = 'ZNE'
self.picks = {}
# initialize panning attributes
self.press = None
self.xpress = None
self.ypress = None
# set attribute holding data
if data is None:
@ -259,17 +269,50 @@ class PickDlg(QDialog):
# plot data
self.getPlotWidget().plotWFData(wfdata=self.getWFData(),
title=self.getStation())
self.limits = {'xlims' : self.getPlotWidget().axes.get_xlim(),
'ylims' : self.getPlotWidget().axes.get_ylim()}
self.apd = self.getWFData()
# set plot labels
self.setPlotLabels()
# connect button press event to an action
self.cidpress = self.connectPressEvent(self.panPress)
self.cidmotion = self.connectMotionEvent()
self.cidrelease = self.connectReleaseEvent()
self.cidscroll = self.connectScrollEvent()
def setupUi(self):
# create matplotlib toolbar to inherit functionality
self.figToolBar = NavigationToolbar2QTAgg(self.getPlotWidget(), self)
self.figToolBar.hide()
# create icons
filter_icon = QIcon()
filter_icon.addPixmap(QPixmap(':/icons/filter.png'))
zoom_icon = QIcon()
zoom_icon.addPixmap(QPixmap(':/icons/zoom.png'))
# create actions
self.filterAction = createAction(parent=self, text='Filter',
slot=self.filterWFData,
icon=QIcon(':/filter.png'),
tip='Filter waveforms',
icon=filter_icon,
tip='Toggle filtered/original'
' waveforms',
checkable=True)
self.selectPhase = QComboBox()
self.selectPhase.addItems(['Pn', 'Pg', 'P1', 'P2'])
self.selectPhase.addItems([None, 'Pn', 'Pg', 'P1', 'P2'])
self.zoomAction = createAction(parent=self, text='Zoom',
slot=self.zoom, icon=zoom_icon,
tip='Zoom into waveform',
checkable=True)
# layout the outermost appearance of the Pick Dialog
_outerlayout = QVBoxLayout()
@ -280,19 +323,64 @@ class PickDlg(QDialog):
_dialtoolbar.addAction(self.filterAction)
_dialtoolbar.addWidget(self.selectPhase)
_innerlayout = QHBoxLayout()
_innerlayout = QVBoxLayout()
_toolslayout = QVBoxLayout()
_toolslabel = QLabel('Place for Tools')
_toolslayout.addWidget(_toolslabel)
_innerlayout.addLayout(_toolslayout)
_innerlayout.addWidget(self.multicompfig)
_buttonbox = QDialogButtonBox(QDialogButtonBox.Apply |
QDialogButtonBox.Ok |
QDialogButtonBox.Cancel)
_innerlayout.addWidget(_buttonbox)
_outerlayout.addWidget(_dialtoolbar)
_outerlayout.addLayout(_innerlayout)
self.selectPhase.currentIndexChanged.connect(self.verifyPhaseSelection)
self.setLayout(_outerlayout)
def disconnectPressEvent(self):
self.getPlotWidget().mpl_disconnect(self.cidpress)
def connectPressEvent(self, slot):
widget = self.getPlotWidget()
return widget.mpl_connect('button_press_event', slot)
def reconnectPressEvent(self, slot):
self.disconnectPressEvent()
return self.connectPressEvent(slot)
def disconnectScrollEvent(self):
widget = self.getPlotWidget()
widget.mpl_disconnect(self.cidscroll)
def connectScrollEvent(self):
widget = self.getPlotWidget()
return widget.mpl_connect('scroll_event', self.scrollZoom)
def disconnectMotionEvent(self):
widget = self.getPlotWidget()
widget.mpl_disconnect(self.cidmotion)
def connectMotionEvent(self):
widget = self.getPlotWidget()
return widget.mpl_connect('motion_notify_event', self.panMotion)
def disconnectReleaseEvent(self):
widget = self.getPlotWidget()
widget.mpl_disconnect(self.cidrelease)
def connectReleaseEvent(self):
widget = self.getPlotWidget()
return widget.mpl_connect('button_release_event', self.panRelease)
def verifyPhaseSelection(self):
phase = self.selectPhase.currentText()
if phase:
self.disconnectReleaseEvent()
self.disconnectScrollEvent()
self.disconnectMotionEvent()
self.reconnectPressEvent(self.setIniPick)
def getComponents(self):
return self.components
@ -303,14 +391,196 @@ class PickDlg(QDialog):
def getPlotWidget(self):
return self.multicompfig
def getChannelID(self, key):
return self.getPlotWidget().getPlotDict()[int(key)][1]
def getWFData(self):
return self.data
def filterWFData(self):
data = self.getWFData().copy().filter(type='bandpass', freqmin=.5, freqmax=15.)
title = self.getStation() + ' (filtered)'
self.getPlotWidget().plotWFData(wfdata=data, title=title)
def selectWFData(self, channel):
component = channel[-1].upper()
wfdata = Stream()
def selectTrace(trace, components):
if trace.stats.channel[-1].upper() in components:
return trace
if component == 'E' or component == 'N':
for trace in self.getWFData():
trace = selectTrace(trace, 'NE')
if trace:
wfdata.append(trace)
elif component == 'Z':
wfdata = self.getWFData().select(component=component)
return wfdata
def getPicks(self):
return self.picks
def getAPD(self):
return self.apd
def updateAPD(self, wfdata):
self.apd = wfdata
def setIniPick(self, gui_event):
channel = self.getChannelID(round(gui_event.ydata))
wfdata = self.selectWFData(channel)
self.disconnectScrollEvent()
self.cidpress = self.reconnectPressEvent(self.setPick)
ini_pick = gui_event.xdata
# calculate the resolution window width from SNR
# SNR >= 3 -> 2 sec HRW
# 3 > SNR >= 2 -> 5 sec MRW
# 2 > SNR >= 1.5 -> 10 sec LRW
# 1.5 > SNR -> 15 sec VLRW
# see also Diehl et al. 2009
res_wins = {
'HRW' : 2.,
'MRW' : 5.,
'LRW' : 10.,
'VLRW' : 15.
}
result = getSNR(wfdata, (10., 2., 1.5), ini_pick)
snr = result[0]
noiselevel = result[2] * 1.5
if snr < 1.5:
x_res = res_wins['VLRW']
elif snr < 2.:
x_res = res_wins['LRW']
elif snr < 3.:
x_res = res_wins['MRW']
else:
x_res = res_wins['HRW']
x_res /= 2
zoomx = [ini_pick - x_res, ini_pick + x_res]
zoomy = [noiselevel * 1.5, -noiselevel * 1.5]
self.getPlotWidget().plotWFData(wfdata=wfdata,
title=self.getStation() +
' picking mode',
zoomx=zoomx,
zoomy=zoomy)
self.updateAPD(wfdata)
# reset labels
self.setPlotLabels()
def setPick(self, gui_event):
pick = gui_event.xdata
ax = self.getPlotWidget().axes
ylims = ax.get_ylim()
ax.plot([pick, pick], ylims, 'r--')
self.getPlotWidget().draw()
def panPress(self, gui_event):
ax = self.getPlotWidget().axes
if gui_event.inaxes != ax: return
self.cur_xlim = ax.get_xlim()
self.cur_ylim = ax.get_ylim()
self.press = gui_event.xdata, gui_event.ydata
self.xpress, self.ypress = self.press
def panRelease(self, gui_event):
ax = self.getPlotWidget().axes
self.press = None
ax.figure.canvas.draw()
def panMotion(self, gui_event):
ax = self.getPlotWidget().axes
if self.press is None: return
if gui_event.inaxes != ax: return
dx = gui_event.xdata - self.xpress
dy = gui_event.ydata - self.ypress
self.cur_xlim -= dx
self.cur_ylim -= dy
ax.set_xlim(self.cur_xlim)
ax.set_ylim(self.cur_ylim)
ax.figure.canvas.draw()
def filterWFData(self):
ax = self.getPlotWidget().axes
ylims = ax.get_ylim()
xlims = ax.get_xlim()
if self.filterAction.isChecked():
data = self.getAPD().copy()
data.filter(type='bandpass', freqmin=.5, freqmax=15.)
title = self.getStation() + ' (filtered)'
else:
data = self.getAPD().copy()
title = self.getStation()
self.getPlotWidget().plotWFData(wfdata=data, title=title, zoomx=xlims,
zoomy=ylims)
self.setPlotLabels()
def setPlotLabels(self):
# get channel labels
pos = self.getPlotWidget().getPlotDict().keys()
labels = [self.getPlotWidget().getPlotDict()[key][1] for key in pos]
# set channel labels
self.getPlotWidget().setYTickLabels(pos, labels)
def zoom(self):
if self.zoomAction.isChecked():
self.disconnectPressEvent()
self.figToolBar.zoom()
else:
self.connectPressEvent(self.setIniPick)
def scrollZoom(self, gui_event, factor=2.):
widget = self.getPlotWidget()
curr_xlim = widget.axes.get_xlim()
curr_ylim = widget.axes.get_ylim()
if gui_event.button == 'up':
scale_factor = 1/factor
elif gui_event.button == 'down':
# deal with zoom out
scale_factor = factor
else:
# deal with something that should never happen
scale_factor = 1
print gui_event.button
new_xlim = gui_event.xdata - scale_factor * (gui_event.xdata - curr_xlim)
new_ylim = gui_event.ydata - scale_factor * (gui_event.ydata - curr_ylim)
new_xlim.sort()
new_xlim[0] = max(new_xlim[0], self.limits['xlims'][0])
new_xlim[1] = min(new_xlim[1], self.limits['xlims'][1])
new_ylim.sort()
new_ylim[0] = max(new_ylim[0], self.limits['ylims'][0])
new_ylim[1] = min(new_ylim[1], self.limits['ylims'][1])
widget.axes.set_xlim(new_xlim)
widget.axes.set_ylim(new_ylim)
widget.draw()
def apply(self):
picks = self.getPicks()
for pick in picks:
print pick
def reject(self):
QDialog.reject(self)
def accept(self):
self.apply()
QDialog.accept(self)
class PropertiesDlg(QDialog):

File diff suppressed because one or more lines are too long

0
testPickDlg.py Normal file → Executable file
View File