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

QtPyLoT.py

@@ -235,10 +235,10 @@ class MainWindow(QMainWindow):
         fileToolBar.setObjectName("FileTools")
         self.addActions(fileToolBar, fileToolActions)
 
-        phaseToolBar = self.addToolBar("PhaseTools")
+        # phaseToolBar = self.addToolBar("PhaseTools")
-        phaseToolActions = (self.selectPAction, self.selectSAction)
+        # phaseToolActions = (self.selectPAction, self.selectSAction)
-        phaseToolBar.setObjectName("PhaseTools")
+        # phaseToolBar.setObjectName("PhaseTools")
-        self.addActions(phaseToolBar, phaseToolActions)
+        # self.addActions(phaseToolBar, phaseToolActions)
 
         # create button group for component selection
 
@@ -536,39 +536,31 @@ class MainWindow(QMainWindow):
         self.plotWaveformData()
         self.drawPicks()
 
+    def pushFilterWF(self, param_args):
+        self.getData().filterWFData(param_args)
+
     def filterWaveformData(self):
         if self.getData():
-            def hasfreq(kwdict):
+            if self.getFilterOptions() and self.filterAction.isChecked():
-                for key in kwdict.keys():
+                kwargs = self.getFilterOptions().parseFilterOptions()
-                    if not key.startswith('freq'):
+                self.pushFilterWF(kwargs)
-                        return True
+            elif self.filterAction.isChecked():
-                return False
+                self.adjustFilterOptions()
-
-            if self.filterAction.isChecked():
-                kwargs = {}
-                freq = self.getFilterOptions().getFreq()
-                if freq is not None and len(freq) > 1:
-                    kwargs['freqmin'] = freq[0]
-                    kwargs['freqmax'] = freq[1]
-                elif freq is not None and len(freq) == 1:
-                    kwargs['freq'] = freq
-                if hasfreq(kwargs):
-                    kwargs['type'] = self.getFilterOptions().getFilterType()
-                    kwargs['corners'] = self.getFilterOptions().getOrder()
-                    self.getData().filterWFData(kwargs)
             else:
                 self.getData().resetWFData()
         self.plotWaveformData()
 
     def adjustFilterOptions(self):
-        filteroptions = self.getFilterOptions()
         fstring = "Filter Options ({0})".format(self.getSeismicPhase())
         filterDlg = FilterOptionsDialog(titleString=fstring,
-                                        parent=self,
+                                        parent=self)
-                                        filterOptions=filteroptions)
         if filterDlg.exec_():
             filteroptions = filterDlg.getFilterOptions()
             self.setFilterOptions(filteroptions)
+            if self.filterAction.isChecked():
+                kwargs = self.getFilterOptions().parseFilterOptions()
+                self.pushFilterWF(kwargs)
+                self.plotWaveformData()
 
     def getFilterOptions(self):
         try:

autoPyLoT.py

@@ -5,20 +5,21 @@
 import os
 import argparse
 import glob
+import subprocess
 
-import matplotlib.pyplot as plt
 from obspy.core import read
 from pylot.core.read.data import Data
 from pylot.core.read.inputs import AutoPickParameter
 from pylot.core.util.structure import DATASTRUCTURE
 from pylot.core.pick.autopick import autopickevent
+from pylot.core.pick.utils import writephases
 from pylot.core.util.version import get_git_version as _getVersionString
 
 __version__ = _getVersionString()
 
 
 def autoPyLoT(inputfile):
-    '''
+    """
     Determine phase onsets automatically utilizing the automatic picking
     algorithms by Kueperkoch et al. 2010/2012.
 
@@ -30,7 +31,7 @@ def autoPyLoT(inputfile):
 
     .. rubric:: Example
 
-    '''
+    """
     print '************************************'
     print '*********autoPyLoT starting*********'
     print 'The Python picking and Location Tool'
@@ -55,9 +56,9 @@ def autoPyLoT(inputfile):
 
     if parameter.hasParam('datastructure'):
         datastructure = DATASTRUCTURE[parameter.getParam('datastructure')]()
-        dsfields = {'root':parameter.getParam('rootpath'),
+        dsfields = {'root' :parameter.getParam('rootpath'),
-                    'dpath':parameter.getParam('datapath'),
+                    'dpath' :parameter.getParam('datapath'),
-                    'dbase':parameter.getParam('database')}
+                    'dbase' :parameter.getParam('database')}
 
         exf = ['root', 'dpath', 'dbase']
 
@@ -68,8 +69,29 @@ def autoPyLoT(inputfile):
         datastructure.modifyFields(**dsfields)
         datastructure.setExpandFields(exf)
 
-        # get path to inventory or dataless-seed file with station meta data
+        # check if default location routine NLLoc is available
-        invdir = parameter.getParam('invdir')
+        if parameter.hasParam('nllocbin'):
+            locflag = 1
+            # get NLLoc-root path
+            nllocroot = parameter.getParam('nllocroot')
+            # get path to NLLoc executable
+            nllocbin = parameter.getParam('nllocbin')
+            nlloccall = '%s/NLLoc' % nllocbin
+            # get name of phase file  
+            phasef = parameter.getParam('phasefile')
+            phasefile = '%s/obs/%s' % (nllocroot, phasef)
+            # get name of NLLoc-control file
+            locf = parameter.getParam('locfile')
+            locfile = '%s/run/%s' % (nllocroot, locf)
+            # patter of NLLoc ttimes from location grid
+            ttpat = parameter.getParam('ttpatter')
+            ttpatter = '%s/time/%s' % (nllocroot, ttpat)
+            # patter of NLLoc-output file
+            nllocoutpatter = parameter.getParam('outpatter')
+        else:
+            locflag = 0
+            print ("!!No location routine available, autoPyLoT just picks the events without locating them!!")
+
 
         # multiple event processing
         # read each event in database
@@ -77,14 +99,41 @@ def autoPyLoT(inputfile):
         if not parameter.hasParam('eventID'):
             for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]:
                 data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
-                print 'Working on event %s' %event
+                print 'Working on event %s' % event
                 print data
 
-                wfdat = data.getWFData() # all available streams
+                wfdat = data.getWFData()  # all available streams
                 ##########################################################
                 # !automated picking starts here!
                 picks = autopickevent(wfdat, parameter)
 
+                ##########################################################
+                # locating
+                if locflag == 1:
+                    # write phases to NLLoc-phase file
+                    writephases(picks, 'NLLoc', phasefile)
+
+                    # For locating the events we have to modify the NLLoc-control file!
+                    # create comment line for NLLoc-control file
+                    # NLLoc-output file
+                    nllocout = '%s/loc/%s_%s' % (nllocroot, event, nllocoutpatter)
+                    locfiles = 'LOCFILES %s NLLOC_OBS %s %s 0' % (phasefile, ttpatter, nllocout)
+                    print ("Modifying  NLLoc-control file %s ..." % locfile)
+                    # modification of NLLoc-control file
+                    filedata = None
+                    nllfile = open(locfile, 'r')
+                    filedata = nllfile.read()
+                    if filedata.find(locfiles) < 0:
+                        # replace old command
+                        filedata = filedata.replace('LOCFILES', locfiles)
+                        nllfile = open(locfile, 'w')
+                        nllfile.write(filedata)
+                        nllfile.close()
+
+                    # locate the event
+                    subprocess.call([nlloccall, locfile])
+                ##########################################################
+
                 print '------------------------------------------'
                 print '-----Finished event %s!-----' % event
                 print '------------------------------------------'
@@ -95,10 +144,38 @@ def autoPyLoT(inputfile):
             print 'Working on event ', parameter.getParam('eventID')
             print data
 
-            wfdat = data.getWFData() # all available streams
+            wfdat = data.getWFData()  # all available streams
             ##########################################################
             # !automated picking starts here!
             picks = autopickevent(wfdat, parameter)
+            
+            ##########################################################
+            # locating
+            if locflag == 1:
+                # write phases to NLLoc-phase file
+                writephases(picks, 'NLLoc', phasefile)
+
+                # For locating the event we have to modify the NLLoc-control file!
+                # create comment line for NLLoc-control file NLLoc-output file
+                nllocout = '%s/loc/%s_%s' % (nllocroot, parameter.getParam('eventID'), nllocoutpatter)
+                locfiles = 'LOCFILES %s NLLOC_OBS %s %s 0' % (phasefile, ttpatter, nllocout)
+                print ("Modifying  NLLoc-control file %s ..." % locfile)
+                # modification of NLLoc-control file
+                filedata = None
+                nllfile = open(locfile, 'r')
+                filedata = nllfile.read()
+                if filedata.find(locfiles) < 0:
+                    # replace old command
+            	    filedata = filedata.replace('LOCFILES', locfiles)
+            	    nllfile = open(locfile, 'w')
+            	    nllfile.write(filedata)
+            	    nllfile.close()
+
+                # locate the event
+                subprocess.call([nlloccall, locfile])
+            ##########################################################
+                   
+
 
             print '------------------------------------------'
             print '-------Finished event %s!-------' % parameter.getParam('eventID')

autoPyLoT_local.in

@@ -3,7 +3,6 @@
 %and picking are to be set here!
 %Parameters are optimized for local data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 #main settings#
 /DATA/Insheim                       #rootpath#     %project path
 EVENT_DATA/LOCAL		    #datapath#     %data path
@@ -12,10 +11,20 @@ e0019.048.13                        #eventID#      %event ID for single event pr
 /DATA/Insheim/STAT_INFO             #invdir#       %full path to inventory or dataless-seed file
 PILOT				    #datastructure#  %choose data structure
 0                                   #iplot#        %flag for plotting: 0 none, 1, partly, >1 everything
-AUTOPHASES_AIC_HOS4_ARH             #phasefile#    %name of autoPILOT output phase file 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-AUTOLOC_AIC_HOS4_ARH                #locfile#      %name of autoPILOT output location file 
+#NLLoc settings
+/home/ludger/NLLOC                  #nllocbin#     %path to NLLoc executable
+/home/ludger/NLLOC/Insheim          #nllocroot#    %root of NLLoc-processing directory
+AUTOPHASES.obs                      #phasefile#    %name of autoPyLoT-output phase file for NLLoc
+                                                   %(in nllocroot/obs)
+Insheim_min1d2015.in                #locfile#      %name of autoPyLoT-output control file for NLLoc
+                                                   %(in nllocroot/run)
+ttime                               #ttpatter#     %pattern of NLLoc ttimes from grid
+                                                   %(in nllocroot/times)
+AUTOLOC_nlloc                       #outpatter#    %pattern of NLLoc-output file
+                                                   %(returns 'eventID_outpatter')
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 AUTOFOCMEC_AIC_HOS4_ARH.in          #focmecin#     %name of focmec input file containing polarities
-HYPOSAT                             #locrt#        %location routine used ("HYPOINVERSE" or "HYPOSAT")
 6                                   #pmin#         %minimum required P picks for location
 4                                   #p0min#        %minimum required P picks for location if at least
                                                    %3 excellent P picks are found

pylot/core/pick/utils.py

@@ -931,6 +931,81 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
 
     return returnflag
 
+
+def writephases(arrivals, fformat, filename):
+    '''
+    Function of methods to write phases to the following standard file
+    formats used for locating earthquakes:
+
+    HYPO71, NLLoc, VELEST, HYPOSAT, HYPOINVERSE and hypoDD
+
+    :param: arrivals
+    :type: dictionary containing all phase information including
+           station ID, phase, first motion, weight (uncertainty), 
+           ....
+
+    :param: fformat
+    :type:  string, chosen file format (location routine),
+            choose between NLLoc, HYPO71, HYPOSAT, VELEST, 
+            HYPOINVERSE, and hypoDD
+
+    :param: filename, full path and name of phase file
+    :type: string
+    '''
+
+
+    if fformat == 'NLLoc':
+    	print ("Writing phases to %s for NLLoc" % filename)
+        fid = open("%s" % filename, 'w')
+        # write header
+        fid.write('# EQEVENT:  Label: EQ001  Loc:  X 0.00  Y 0.00  Z 10.00  OT 0.00 \n')
+        for key in arrivals:
+            if arrivals[key]['P']['weight'] < 4:
+                fm = arrivals[key]['P']['fm']
+                onset = arrivals[key]['P']['mpp']
+                year = onset.year
+                month = onset.month
+                day = onset.day
+                hh = onset.hour
+                mm = onset.minute
+                ss = onset.second
+                ms = onset.microsecond
+                ss_ms = ss + ms / 1000000.0
+                fid.write('%s ? ? ? P   %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 1 \n' % (key,
+                 fm,
+                 year,
+                 month,
+                 day,
+                 hh,
+                 mm,
+                 ss_ms))
+            if arrivals[key]['S']['weight'] < 4:
+                fm = '?'
+                onset = arrivals[key]['S']['mpp']
+                year = onset.year
+                month = onset.month
+                day = onset.day
+                hh = onset.hour
+                mm = onset.minute
+                ss = onset.second
+                ms = onset.microsecond
+                ss_ms = ss + ms / 1000000.0
+                fid.write('%s ? ? ? S   %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 1 \n' % (key,
+                 fm,
+                 year,
+                 month,
+                 day,
+                 hh,
+                 mm,
+                 ss_ms))
+
+
+        fid.close()
+
+
+
+        
+
 if __name__ == '__main__':
     import doctest
     doctest.testmod()

pylot/core/read/inputs.py

@@ -206,8 +206,11 @@ class FilterOptions(object):
                    order=self.getOrder())
         return hrs
 
+    def __nonzero__(self):
+        return bool(self.getFilterType())
+
     def parseFilterOptions(self):
-        if self.getFilterType():
+        if self:
             robject = {'type': self.getFilterType(), 'corners': self.getOrder()}
             if len(self.getFreq()) > 1:
                 robject['freqmin'] = self.getFreq()[0]

pylot/core/util/utils.py

@@ -32,7 +32,8 @@ def runProgram(cmd, parameter=None):
     output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd),
                                      shell = True)
 
-
+def isSorted(iterable):
+    return sorted(iterable) == iterable
 
 def fnConstructor(s):
     if type(s) is str:

pylot/core/util/widgets.py

@@ -27,7 +27,7 @@ from pylot.core.pick.utils import getSNR, earllatepicker, getnoisewin,\
     getResolutionWindow
 from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS
 from pylot.core.util.utils import prepTimeAxis, getGlobalTimes, scaleWFData, \
-    demeanTrace
+    demeanTrace, isSorted
 
 
 def createAction(parent, text, slot=None, shortcut=None, icon=None,
@@ -985,7 +985,7 @@ class FilterOptionsDialog(QDialog):
         """
         super(FilterOptionsDialog, self).__init__()
 
-        if parent is not None:
+        if parent is not None and parent.getFilterOptions():
             self.filterOptions = parent.getFilterOptions()
         elif filterOptions is not None:
             self.filterOptions = FilterOptions(filterOptions)
@@ -1021,8 +1021,8 @@ class FilterOptionsDialog(QDialog):
             try:
                 self.freqmaxSpinBox.setValue(self.getFilterOptions().getFreq())
                 self.freqminSpinBox.setValue(self.getFilterOptions().getFreq())
-            except TypeError, e:
+            except TypeError as e:
-                print e
+                print(e)
                 self.freqmaxSpinBox.setValue(1.)
                 self.freqminSpinBox.setValue(.1)
 
@@ -1037,6 +1037,7 @@ class FilterOptionsDialog(QDialog):
         self.selectTypeLabel.setText("Select filter type:")
         self.selectTypeCombo = QComboBox()
         self.selectTypeCombo.addItems(typeOptions)
+        self.selectTypeCombo.setCurrentIndex(typeOptions.index(self.getFilterOptions().getFilterType()))
         self.selectTypeLayout = QVBoxLayout()
         self.selectTypeLayout.addWidget(self.orderLabel)
         self.selectTypeLayout.addWidget(self.orderSpinBox)
@@ -1071,30 +1072,28 @@ class FilterOptionsDialog(QDialog):
         self.buttonBox.rejected.connect(self.reject)
 
     def updateUi(self):
-        _enable = False
+        type = self.selectTypeCombo.currentText()
-        if self.selectTypeCombo.currentText() not in ['bandpass', 'bandstop']:
+        _enable = type in ['bandpass', 'bandstop']
-            self.freqminLabel.setText("cutoff:")
+        freq = [self.freqminSpinBox.value(), self.freqmaxSpinBox.value()]
-            self.freqmaxSpinBox.setValue(self.freqminSpinBox.value())
-        else:
-            _enable = True
-            self.freqminLabel.setText("minimum:")
-
         self.freqmaxLabel.setEnabled(_enable)
         self.freqmaxSpinBox.setEnabled(_enable)
 
-        self.getFilterOptions().setFilterType(
+        if not _enable:
-            self.selectTypeCombo.currentText())
+            self.freqminLabel.setText("cutoff:")
-        freq = [self.freqminSpinBox.value()]
+            self.freqmaxSpinBox.setValue(freq[0])
-        if _enable:
+            freq.remove(freq[1])
-            if self.freqminSpinBox.value() > self.freqmaxSpinBox.value():
+        else:
+            self.freqminLabel.setText("minimum:")
+            if not isSorted(freq):
                 QMessageBox.warning(self, "Value error",
                                     "Maximum frequency must be at least the "
                                     "same value as minimum frequency (notch)!")
-                self.freqmaxSpinBox.setValue(self.freqminSpinBox.value())
+                self.freqmaxSpinBox.setValue(freq[0])
                 self.freqmaxSpinBox.selectAll()
                 self.freqmaxSpinBox.setFocus()
                 return
-            freq.append(self.freqmaxSpinBox.value())
+
+        self.getFilterOptions().setFilterType(type)
         self.getFilterOptions().setFreq(freq)
         self.getFilterOptions().setOrder(self.orderSpinBox.value())