Pulling latest release from server ariadne.

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

autoPyLoT.py

@@ -0,0 +1,90 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+
+import os
+import argparse
+
+from pylot.core.util import _getVersionString
+from pylot.core.read import Data, AutoPickParameter
+from pylot.core.pick.CharFuns import HOScf, AICcf
+from pylot.core.util.structure import DATASTRUCTURE
+
+
+__version__ = _getVersionString()
+
+METHOD = {'HOS':HOScf, 'AIC':AICcf}
+
+def autoPyLoT(fnames, inputfile):
+    '''
+    Determine phase onsets automatically utilizing the automatic picking
+    algorithm by Kueperkoch et al. 2011.
+
+    :param fnames: list of strings containing the paths or urls to the
+    waveform data to be picked
+    :type fnames: list
+    :param inputfile: path to the input file containing all parameter
+    information for automatic picking (for formatting details, see.
+    `~pylot.core.read.input.AutoPickParameter`
+    :type inputfile: str
+    :return:
+
+    .. rubric:: Example
+
+    '''
+    parameter = AutoPickParameter(inputfile)
+
+    data = Data()
+
+    meth = parameter.getParam('algoP')
+    tsnr1 = parameter.getParam('tsnr1')
+    tsnr2 = parameter.getParam('tsnr2')
+    tnoise = parameter.getParam('pnoiselen')
+    tsignal = parameter.getParam('tlim')
+    order = parameter.getParam('hosorder')
+    thosmw = parameter.getParam('tlta')
+
+    if parameter.hasParam('datastructure'):
+        datastructure = DATASTRUCTURE[parameter.getParam('datastructure')]()
+        dsfields = {'root':parameter.getParam('rootpath')} #see TODO: in data.py
+        datastructure.modifyFields(dsfields)
+
+        def expandSDS(datastructure):
+            return datastructure.expandDataPath()
+
+        def expandPDS(datastructure):
+            return os.path.join(datastructure.expandDataPath(),'*')
+
+        dsem = {'PDS':expandPDS, 'SDS':expandSDS}
+
+        expandMethod = dsem[datastructure.getName()]
+
+        data.setWFData(expandMethod())
+    else:
+        data.setWFData()
+
+    cfP = METHOD[meth](data.getWFData(), (tnoise, tsignal), thosmw, order)
+
+if __name__ == "__main__":
+    # parse arguments
+    parser = argparse.ArgumentParser(
+        description='''This program visualizes Probabilistic Power Spectral Densities.''')
+
+    parser.add_argument('-d', '-D', '--data', type=str, action='store',
+                        help='''path or url to the data to be picked''',
+                        default='http://examples.obspy.org/BW.KW1..EHZ.D.2011.037'
+                        )
+    parser.add_argument('-i', '-I', '--inputfile', type=str,
+                        action='store',
+                        help='''full path to the file containing the input
+                        parameters for autoPyLoT''',
+                        default=os.path.join(os.path.expanduser('~'), '.pylot',
+                                             'autoPyLoT.in')
+                        )
+    parser.add_argument('-v', '-V', '--version', action='version',
+                        version='autoPyLoT ' + __version__,
+                        help='show version information and exit')
+
+    cla = parser.parse_args()
+
+    autoPyLoT(cla.data, str(cla.inputfile))

pylot/RELEASE-VERSION

@@ -1 +1 @@
-043c-dirty
+694a-dirty

pylot/core/pick/utils.py

@@ -53,32 +53,22 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
     #get signal window
     isignal = getsignalwin(t, Pick1, TSNR[2])
     #calculate noise level
-    nlevel = max(abs(x[inoise])) * nfac
+    nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac
     #get time where signal exceeds nlevel
-    ilup = np.where(x[isignal] > nlevel)
+    ilup, = np.where(x[isignal] > nlevel)
-    ildown = np.where(x[isignal] < -nlevel)
+    ildown, = np.where(x[isignal] < -nlevel)
-    if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
+    if not ilup.size and not ildown.size:
-        print 'earllatepicker: Signal lower than noise level, misspick?'
+        raise ValueError('earllatepicker: Signal lower than noise level')
-        return
+    il = min(np.min(ilup) if ilup.size else float('inf'),
-    il = min([ilup[0][0], ildown[0][0]])
+             np.min(ildown) if ildown.size else float('inf'))
     LPick = t[isignal][il]
 
     #get earliest possible pick
-    #get next 2 zero crossings after most likely pick
+
-    #initial onset is assumed to be the first zero crossing
+    #determine all zero crossings in signal window
-    zc = []
+    zc = crossings_nonzero_all(x[isignal])
-    zc.append(Pick1)
+    #calculate mean half period T0 of signal as the average of the
-    i = 0
+    T0 = np.mean(np.diff(zc)) * X[0].stats.delta  #this is half wave length!
-    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/4 is assumed as time difference between most likely and earliest possible pick!
     EPick = Pick1 - T0 / 2
 
@@ -288,6 +278,10 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
 
     return FM
 
+def crossings_nonzero_all(data):
+    pos = data > 0
+    npos = ~pos
+    return ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0]
 
 def getSNR(X, TSNR, t1):
     '''

pylot/core/read/data.py

@@ -3,15 +3,13 @@
 
 import os
 
-import numpy as np
 from obspy.core import (read, Stream, UTCDateTime)
 from obspy import readEvents
 from obspy.core.event import (Event, Catalog)
 
 from pylot.core.read import readPILOTEvent
-
+from pylot.core.util import fnConstructor, FormatError, \
-from pylot.core.util import fnConstructor, createEvent, FormatError, \
+    getGlobalTimes
-    prepTimeAxis, getGlobalTimes
 
 
 class Data(object):
@@ -129,8 +127,8 @@ class Data(object):
             except TypeError:
                 try:
                     self.wfdata += read(fname, format='GSE2')
-                except Exception:
+                except Exception, e:
-                    warnmsg += '{0}\n'.format(fname)
+                    warnmsg += '{0}\n{1}\n'.format(fname, e)
         if warnmsg:
             warnmsg = 'WARNING: unable to read\n' + warnmsg
             print warnmsg
@@ -148,6 +146,7 @@ class Data(object):
     def getEvtData(self):
         return self.evtdata
 
+#TODO: write superclass DataStructure instead of three fully featured individual classes
 
 class GenericDataStructure(object):
     '''
@@ -157,6 +156,7 @@ class GenericDataStructure(object):
 
     def __init__(self, structexp='$R/$D/$E', folderdepth=2, **kwargs):
         structureOptions = ('$R', '$D', '$E')
+        self.__name = 'GDS'
         structExpression = []
         depth = 0
         while structexp is not os.path.sep:
@@ -182,6 +182,9 @@ class GenericDataStructure(object):
             key = str(key).upper()
             self.__gdsFields[key] = value
 
+    def getName(self):
+        return self.__name
+
     def getFields(self):
         return self.__gdsFields
 
@@ -198,7 +201,9 @@ class PilotDataStructure(object):
     def __init__(self, dataformat='GSE2', fsuffix='gse',
                  root='/data/Egelados/EVENT_DATA/LOCAL/', database='2006.01',
                  **kwargs):
-        self.dataType = dataformat
+        self.__name = 'PDS'
+        self.dataFormat = dataformat
+        self.dataType = 'waveform'
         self.__pdsFields = {'ROOT': root,
                             'DATABASE': database,
                             'SUFFIX': fsuffix
@@ -233,6 +238,9 @@ class PilotDataStructure(object):
     def getFields(self):
         return self.__pdsFields
 
+    def getName(self):
+        return self.__name
+
     def expandDataPath(self):
         datapath = os.path.join(self.getFields()['ROOT'],
                                 self.getFields()['DATABASE'])
@@ -260,8 +268,8 @@ class SeiscompDataStructure(object):
     }
 
     def __init__(self, dataType='waveform', sdate=None, edate=None, **kwargs):
-        # imports
+
-        from obspy.core import UTCDateTime
+        self.__name = 'SDS'
 
         def checkDate(date):
             if not isinstance(date, UTCDateTime):
@@ -337,6 +345,9 @@ class SeiscompDataStructure(object):
     def getFields(self):
         return self.__sdsFields
 
+    def getName(self):
+        return self.__name
+
     def expandDataPath(self):
         fullChan = '{0}.{1}'.format(self.getFields()['CHAN'], self.getType())
         dataPath = os.path.join(self.getFields()['SDSdir'],

pylot/core/read/inputs.py

@@ -5,7 +5,7 @@
 class AutoPickParameter(object):
     '''
     AutoPickParameters is a parameter type object capable to read and/or write
-    parameter ASCII and binary files.
+    parameter ASCII.
 
     :param fn str: Filename of the input file
 
@@ -34,51 +34,64 @@ class AutoPickParameter(object):
     ==========  ==========  =======================================
     '''
 
-    def __init__(self, fn=None):
+    def __init__(self, fnin=None, fnout=None, **kwargs):
         '''
         Initialize parameter object:
 
         read content of an ASCII file an form a type consistent dictionary
         contain all parameters.
         '''
-        self.__filename = fn
+
+        self.__filename = fnin
         parFileCont = {}
-        try:
+        # read from parsed arguments alternatively
-            if self.__filename is not None:
+        for key, val in kwargs.iteritems():
+            parFileCont[key] = val
+
+        if self.__filename is not None:
                 inputFile = open(self.__filename, 'r')
-                lines = inputFile.readlines()
+        else:
-                for line in lines:
+            return
-                    parspl = line.split('\t')[:2]
+        try:
-                    parFileCont[parspl[0].strip()] = parspl[1]
+            lines = inputFile.readlines()
-                for key, value in parFileCont.iteritems():
+            for line in lines:
-                    try:
+                parspl = line.split('\t')[:2]
-                        val = int(value)
+                parFileCont[parspl[0].strip()] = parspl[1]
-                    except:
-                        try:
-                            val = float(value)
-                        except:
-                            if value.find(';') > 0:
-                                vallist = value.strip().split(';')
-                                val = []
-                                for val0 in vallist:
-                                    val0 = float(val0)
-                                    val.append(val0)
-                            else:
-                                val = str(value.strip())
-                    parFileCont[key] = val
-            else:
-                parFileCont = {}
         except Exception, e:
             self._printParameterError(e)
-            parFileCont = {}
+            inputFile.seek(0)
+            lines = inputFile.readlines()
+            for line in lines:
+                if not line.startswith(('#', '%', '\n', ' ')):
+                    parspl = line.split('#')[:2]
+                    parFileCont[parspl[1].strip()] = parspl[0].strip()
+        for key, value in parFileCont.iteritems():
+            try:
+                val = int(value)
+            except:
+                try:
+                    val = float(value)
+                except:
+                    if len(value.split(' ')) > 1:
+                        vallist = value.strip().split(' ')
+                        val = []
+                        for val0 in vallist:
+                            val0 = float(val0)
+                            val.append(val0)
+                    else:
+                        val = str(value.strip())
+            parFileCont[key] = val
         self.__parameter = parFileCont
 
+        if fnout:
+            self.export2File(fnout)
+
     # Human-readable string representation of the object
     def __str__(self):
         string = ''
         string += 'Automated picking parameter:\n\n'
         if self.__parameter:
-            for key, value in self.__parameter.iteritems():
+            for key, value in self.iteritems():
                 string += '%s:\t\t%s\n' % (key, value)
         else:
             string += 'Empty parameter dictionary.'
@@ -107,6 +120,25 @@ class AutoPickParameter(object):
     def __len__(self):
         return len(self.__parameter.keys())
 
+    def iteritems(self):
+        for key, value in self.__parameter.iteritems():
+            yield key, value
+
+    def hasParam(self, *args):
+
+        def test(param):
+            try:
+                self.__parameter[param]
+                return True
+            except KeyError:
+                return False
+
+        try:
+            for param in args:
+                test(param)
+        except TypeError:
+            test(args)
+
     def getParam(self, *args):
         try:
             for param in args:
@@ -122,15 +154,19 @@ class AutoPickParameter(object):
 
     def setParam(self, **kwargs):
         for param, value in kwargs.iteritems():
-            try:
+            self.__setitem__(param, value)
-                self.__setitem__(param, value)
-            except KeyError, e:
-                self._printParameterError(e)
         print self
 
     def _printParameterError(self, errmsg):
         print 'ParameterError:\n non-existent parameter %s' % errmsg
 
+    def export2File(self, fnout):
+        fid_out = open(fnout, 'w')
+        lines = []
+        for key, value in self.iteritems():
+            lines.append('{key}\t{value}'.format(key=key, value=value))
+        fid_out.writelines(lines)
+
 
 class FilterOptions(object):
     '''

pylot/core/util/widgets.py

@@ -285,9 +285,9 @@ class PickDlg(QDialog):
 
         # connect button press event to an action
         self.cidpress = self.connectPressEvent(self.panPress)
-        self.cidmotion = self.connectMotionEvent()
+        self.cidmotion = self.connectMotionEvent(self.panMotion)
-        self.cidrelease = self.connectReleaseEvent()
+        self.cidrelease = self.connectReleaseEvent(self.panRelease)
-        self.cidscroll = self.connectScrollEvent()
+        self.cidscroll = self.connectScrollEvent(self.scrollZoom)
 
     def setupUi(self):
 
@@ -328,6 +328,7 @@ class PickDlg(QDialog):
 
         _dialtoolbar.addAction(self.filterAction)
         _dialtoolbar.addWidget(self.selectPhase)
+        _dialtoolbar.addAction(self.zoomAction)
 
         _innerlayout = QVBoxLayout()
 
@@ -346,39 +347,40 @@ class PickDlg(QDialog):
         self.setLayout(_outerlayout)
 
     def disconnectPressEvent(self):
-        self.getPlotWidget().mpl_disconnect(self.cidpress)
+        widget = self.getPlotWidget()
+        widget.mpl_disconnect(self.cidpress)
+        self.cidpress = None
 
     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)
+        self.cidscroll = None
 
-    def connectScrollEvent(self):
+    def connectScrollEvent(self, slot):
         widget = self.getPlotWidget()
-        return widget.mpl_connect('scroll_event', self.scrollZoom)
+        return widget.mpl_connect('scroll_event', slot)
 
     def disconnectMotionEvent(self):
         widget = self.getPlotWidget()
         widget.mpl_disconnect(self.cidmotion)
+        self.cidmotion = None
 
-    def connectMotionEvent(self):
+    def connectMotionEvent(self, slot):
         widget = self.getPlotWidget()
-        return widget.mpl_connect('motion_notify_event', self.panMotion)
+        return widget.mpl_connect('motion_notify_event', slot)
 
     def disconnectReleaseEvent(self):
         widget = self.getPlotWidget()
         widget.mpl_disconnect(self.cidrelease)
+        self.cidrelease = None
 
-    def connectReleaseEvent(self):
+    def connectReleaseEvent(self, slot):
         widget = self.getPlotWidget()
-        return widget.mpl_connect('button_release_event', self.panRelease)
+        return widget.mpl_connect('button_release_event', slot)
 
     def verifyPhaseSelection(self):
         phase = self.selectPhase.currentText()
@@ -386,12 +388,13 @@ class PickDlg(QDialog):
             self.disconnectReleaseEvent()
             self.disconnectScrollEvent()
             self.disconnectMotionEvent()
-            self.reconnectPressEvent(self.setIniPick)
+            self.disconnectPressEvent()
+            self.cidpress = self.connectPressEvent(self.setIniPick)
         else:
             self.cidpress = self.connectPressEvent(self.panPress)
-            self.cidmotion = self.connectMotionEvent()
+            self.cidmotion = self.connectMotionEvent(self.panMotion)
-            self.cidrelease = self.connectReleaseEvent()
+            self.cidrelease = self.connectReleaseEvent(self.panRelease)
-            self.cidscroll = self.connectScrollEvent()
+            self.cidscroll = self.connectScrollEvent(self.scrollZoom)
 
     def getComponents(self):
         return self.components
@@ -438,8 +441,10 @@ class PickDlg(QDialog):
         wfdata = self.selectWFData(channel)
 
         self.disconnectScrollEvent()
-
+        self.disconnectPressEvent()
-        self.cidpress = self.reconnectPressEvent(self.setPick)
+        self.disconnectReleaseEvent()
+        self.disconnectMotionEvent()
+        self.cidpress = self.connectPressEvent(self.setPick)
 
         ini_pick = gui_event.xdata
 
@@ -457,7 +462,7 @@ class PickDlg(QDialog):
             'VLRW' : 15.
         }
 
-        result = getSNR(wfdata, (5., .5, 1.), ini_pick)
+        result = getSNR(wfdata, (5., .5, 2.), ini_pick)
 
         snr = result[0]
         noiselevel = result[2] * 1.5
@@ -473,7 +478,7 @@ class PickDlg(QDialog):
         x_res /= 2
 
         zoomx = [ini_pick - x_res, ini_pick + x_res]
-        zoomy = [noiselevel * 1.5, -noiselevel * 1.5]
+        zoomy = [-noiselevel * 1.5, noiselevel * 1.5]
         self.getPlotWidget().plotWFData(wfdata=wfdata,
                                         title=self.getStation() +
                                               ' picking mode',
@@ -493,7 +498,7 @@ class PickDlg(QDialog):
 
         wfdata = self.getAPD().copy().select(channel=channel)
         # get earliest and latest possible pick
-        [epp, lpp, pickerror] = earllatepicker(wfdata, 1.5, (5., .5, 1.), pick)
+        [epp, lpp, pickerror] = earllatepicker(wfdata, 1.5, (5., .5, 2.), pick)
 
         # plotting picks
         ax = self.getPlotWidget().axes
@@ -559,9 +564,11 @@ class PickDlg(QDialog):
     def zoom(self):
         if self.zoomAction.isChecked():
             self.disconnectPressEvent()
+            self.disconnectMotionEvent()
+            self.disconnectReleaseEvent()
             self.figToolBar.zoom()
         else:
-            self.connectPressEvent(self.setIniPick)
+            self.cidpress = self.connectPressEvent(self.setIniPick)
 
     def scrollZoom(self, gui_event, factor=2.):