[add] for global seismology CF pick windows will now be calculated

relative to estimated tt from TauPy, metadata and source location (in XML
file) needed
[TO DO]: automatic export of XML (esp. source loc) before autopicking

QtPyLoT.py

@@ -67,14 +67,13 @@ from pylot.core.pick.compare import Comparison
 from pylot.core.pick.utils import symmetrize_error
 from pylot.core.io.phases import picksdict_from_picks
 import pylot.core.loc.nll as nll
-from pylot.core.util.defaults import FILTERDEFAULTS, SetChannelComponents, \
+from pylot.core.util.defaults import FILTERDEFAULTS, SetChannelComponents
-    readFilterInformation
 from pylot.core.util.errors import FormatError, DatastructureError, \
     OverwriteError
 from pylot.core.util.connection import checkurl
 from pylot.core.util.dataprocessing import read_metadata, restitute_data
 from pylot.core.util.utils import fnConstructor, getLogin, \
-    full_range
+    full_range, readFilterInformation
 from pylot.core.util.event import Event
 from pylot.core.io.location import create_creation_info, create_event
 from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \
@@ -722,7 +721,7 @@ class MainWindow(QMainWindow):
         fext = '.xml'
         for event in events:
             path = event.path
-            eventname = path.split('/')[-1]
+            eventname = path.split('/')[-1] #or event.pylot_id
             filename = os.path.join(path, 'PyLoT_' + eventname + fext)
             if os.path.isfile(filename):
                 self.load_data(filename, draw=False, event=event, overwrite=True)

autoPyLoT.py

@@ -111,8 +111,6 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             print('Parameters set and input file given. Choose either of both.')
             return
 
-    data = Data()
-
     evt = None
 
     # reading parameter file
@@ -192,16 +190,24 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             return
 
         # transform system path separator to '/'
-        for index, event in enumerate(events):
+        for index, eventpath in enumerate(events):
-            event = event.replace(SEPARATOR, '/')
+            eventpath = eventpath.replace(SEPARATOR, '/')
-            events[index] = event
+            events[index] = eventpath
 
-        for event in events:
+        for eventpath in events:
-            pylot_event = Event(event)  # event should be path to event directory
+            evID = os.path.split(eventpath)[-1]
-            data.setEvtData(pylot_event)
+            fext = '.xml'
+            filename = os.path.join(eventpath, 'PyLoT_' + evID + fext)
+            try:
+                data = Data(evtdata=filename)
+                print('Reading event data from filename {}...'.format(filename))
+            except Exception as e:
+                print('Could not read event from file {}: {}'.format(filename, e))
+                data = Data()
+                pylot_event = Event(eventpath)  # event should be path to event directory
+                data.setEvtData(pylot_event)
             if fnames == 'None':
-                data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
+                data.setWFData(glob.glob(os.path.join(datapath, eventpath, '*')))
-                evID = os.path.split(event)[-1]
                 # the following is necessary because within
                 # multiple event processing no event ID is provided
                 # in autopylot.in
@@ -218,7 +224,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             else:
                 data.setWFData(fnames)
 
-                event = events[0]
+                eventpath = events[0]
                 # now = datetime.datetime.now()
                 # evID = '%d%02d%02d%02d%02d' % (now.year,
                 #                               now.month,
@@ -238,16 +244,18 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             corr_dat = restitute_data(wfdat.copy(), *metadata, ncores=ncores)
             if not corr_dat and locflag:
                 locflag = 2
-            print('Working on event %s. Stations: %s' % (event, station))
+            print('Working on event %s. Stations: %s' % (eventpath, station))
             print(wfdat)
             ##########################################################
             # !automated picking starts here!
             if input_dict:
                 if 'fig_dict' in input_dict:
                     fig_dict = input_dict['fig_dict']
-                    picks = autopickevent(wfdat, parameter, iplot=iplot, fig_dict=fig_dict, ncores=ncores)
+                    picks = autopickevent(wfdat, parameter, iplot=iplot, fig_dict=fig_dict,
+                                          ncores=ncores, metadata=metadata, origin=data.get_evt_data().origins)
             else:
-                picks = autopickevent(wfdat, parameter, iplot=iplot, ncores=ncores)
+                picks = autopickevent(wfdat, parameter, iplot=iplot,
+                                      ncores=ncores, metadata=metadata, origin=data.get_evt_data().origins)
             ##########################################################
             # locating
             if locflag > 0:
@@ -393,29 +401,29 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             # ObsPy event object
             data.applyEVTData(picks)
             if evt is not None:
-                event_id = event.split('/')[-1]
+                event_id = eventpath.split('/')[-1]
                 evt.resource_id = ResourceIdentifier('smi:local/' + event_id)
                 data.applyEVTData(evt, 'event')
-            fnqml = '%s/PyLoT_%s' % (event, evID)
+            fnqml = '%s/PyLoT_%s' % (eventpath, evID)
             data.exportEvent(fnqml, fnext='.xml', fcheck=['auto', 'magnitude', 'origin'])
             if locflag == 1:
                 # HYPO71
-                hypo71file = '%s/PyLoT_%s_HYPO71_phases' % (event, evID)
+                hypo71file = '%s/PyLoT_%s_HYPO71_phases' % (eventpath, evID)
                 hypo71.export(picks, hypo71file, parameter)
                 # HYPOSAT
-                hyposatfile = '%s/PyLoT_%s_HYPOSAT_phases' % (event, evID)
+                hyposatfile = '%s/PyLoT_%s_HYPOSAT_phases' % (eventpath, evID)
                 hyposat.export(picks, hyposatfile, parameter)
                 # VELEST
-                velestfile = '%s/PyLoT_%s_VELEST_phases.cnv' % (event, evID)
+                velestfile = '%s/PyLoT_%s_VELEST_phases.cnv' % (eventpath, evID)
                 velest.export(picks, velestfile, parameter, evt)
                 # hypoDD
-                hypoddfile = '%s/PyLoT_%s_hypoDD_phases.pha' % (event, evID)
+                hypoddfile = '%s/PyLoT_%s_hypoDD_phases.pha' % (eventpath, evID)
                 hypodd.export(picks, hypoddfile, parameter, evt)
                 # FOCMEC
-                focmecfile = '%s/PyLoT_%s_FOCMEC.in' % (event, evID)
+                focmecfile = '%s/PyLoT_%s_FOCMEC.in' % (eventpath, evID)
                 focmec.export(picks, focmecfile, parameter, evt)
                 # HASH
-                hashfile = '%s/PyLoT_%s_HASH' % (event, evID)
+                hashfile = '%s/PyLoT_%s_HASH' % (eventpath, evID)
                 hash.export(picks, hashfile, parameter, evt)
 
             endsplash = '''------------------------------------------\n'

pylot/core/io/default_parameters.py

@@ -87,12 +87,14 @@ defaults = {'rootpath': {'type': str,
                    'namestring': ('Quality factor', 'Qp1', 'Qp2')},
 
             'pstart': {'type': float,
-                       'tooltip': 'start time [s] for calculating CF for P-picking',
+                       'tooltip': 'start time [s] for calculating CF for P-picking (if TauPy:'
+                                  ' seconds relative to estimated onset)',
                        'value': 15.0,
                        'namestring': 'P start'},
 
             'pstop': {'type': float,
-                      'tooltip': 'end time [s] for calculating CF for P-picking',
+                      'tooltip': 'end time [s] for calculating CF for P-picking (if TauPy:'
+                                 ' seconds relative to estimated onset)',
                       'value': 60.0,
                       'namestring': 'P stop'},
 
@@ -381,12 +383,12 @@ defaults = {'rootpath': {'type': str,
             'use_taup': {'type': bool,
                          'tooltip': 'use estimated traveltimes from TauPy for calculating windows for CF',
                          'value': True,
-                         'namestring': 'Use Taupy'},
+                         'namestring': 'Use TauPy'},
 
             'taup_model': {'type': str,
                            'tooltip': 'define TauPy model for traveltime estimation',
                            'value': 'iasp91',
-                           'namestring': 'Taupy model'}
+                           'namestring': 'TauPy model'}
             }
 
 settings_main = {

pylot/core/pick/autopick.py

@@ -18,10 +18,13 @@ from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
 from pylot.core.pick.picker import AICPicker, PragPicker
 from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
     getSNR, fmpicker, checkPonsets, wadaticheck
-from pylot.core.util.utils import getPatternLine, gen_Pool
+from pylot.core.util.utils import getPatternLine, gen_Pool, identifyPhase, loopIdentifyPhase, \
+    full_range
+
+from obspy.taup import TauPyModel
 
 
-def autopickevent(data, param, iplot=0, fig_dict=None, ncores=0):
+def autopickevent(data, param, iplot=0, fig_dict=None, ncores=0, metadata=None, origin=None):
     stations = []
     all_onsets = {}
     input_tuples = []
@@ -42,9 +45,11 @@ def autopickevent(data, param, iplot=0, fig_dict=None, ncores=0):
         topick = data.select(station=station)
 
         if not iplot:
-            input_tuples.append((topick, param, apverbose))
+            input_tuples.append((topick, param, apverbose, metadata, origin))
         if iplot > 0:
-            all_onsets[station] = autopickstation(topick, param, verbose=apverbose, iplot=iplot, fig_dict=fig_dict)
+            all_onsets[station] = autopickstation(topick, param, verbose=apverbose,
+                                                  iplot=iplot, fig_dict=fig_dict,
+                                                  metadata=metadata, origin=origin)
 
     if iplot > 0:
         print('iPlot Flag active: NO MULTIPROCESSING possible.')
@@ -69,12 +74,13 @@ def autopickevent(data, param, iplot=0, fig_dict=None, ncores=0):
 
 
 def call_autopickstation(input_tuple):
-    wfstream, pickparam, verbose = input_tuple
+    wfstream, pickparam, verbose, metadata, origin = input_tuple
     # multiprocessing not possible with interactive plotting
-    return autopickstation(wfstream, pickparam, verbose, iplot=0)
+    return autopickstation(wfstream, pickparam, verbose, iplot=0, metadata=metadata, origin=origin)
 
 
-def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None):
+def autopickstation(wfstream, pickparam, verbose=False,
+                    iplot=0, fig_dict=None, metadata=None, origin=None):
     """
     :param wfstream: `~obspy.core.stream.Stream`  containing waveform
     :type wfstream: obspy.core.stream.Stream
@@ -117,6 +123,8 @@ def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None):
     algoS = pickparam.get('algoS')
     sstart = pickparam.get('sstart')
     sstop = pickparam.get('sstop')
+    use_taup = pickparam.get('use_taup')
+    taup_model = pickparam.get('taup_model')
     bph1 = pickparam.get('bph1')
     bph2 = pickparam.get('bph2')
     tsnrh = pickparam.get('tsnrh')
@@ -182,6 +190,8 @@ def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None):
     if len(ndat) == 0:  # check for other components
         ndat = wfstream.select(component="1")
 
+    wfstart, wfend = full_range(wfstream)
+
     if algoP == 'HOS' or algoP == 'ARZ' and zdat is not None:
         msg = '##################################################\nautopickstation:' \
               ' Working on P onset of station {station}\nFiltering vertical ' \
@@ -197,7 +207,47 @@ def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None):
         z_copy[0].data = tr_filt.data
         ##############################################################
         # check length of waveform and compare with cut times
-        Lc = pstop - pstart
+
+        # for global seismology: use tau-p method for estimating travel times (needs source and station coords.)
+        # if not given: sets Lc to infinity to use full stream
+        if use_taup:
+            Lc = np.inf
+            print('autopickstation: use_taup flag active.')
+            if not metadata[1]:
+                print('Warning: Could not use TauPy to estimate onsets as there are no metadata given.')
+            else:
+                if origin:
+                    source_origin = origin[0]
+                    station_id = wfstream[0].get_id()
+                    parser = metadata[1]
+                    station_coords = parser.get_coordinates(station_id)
+                    model = TauPyModel(taup_model)
+                    arrivals = model.get_travel_times_geo(
+                        source_origin.depth,
+                        source_origin.latitude,
+                        source_origin.longitude,
+                        station_coords['latitude'],
+                        station_coords['longitude']
+                    )
+                    phases = {'P': [],
+                              'S': []}
+                    for arr in arrivals:
+                        phases[identifyPhase(loopIdentifyPhase(arr.phase.name))].append(arr)
+
+                    # get first P and S onsets from arrivals list
+                    arrP, estFirstP = min([(arr, arr.time) for arr in phases['P']], key = lambda t: t[1])
+                    arrS, estFirstS = min([(arr, arr.time) for arr in phases['S']], key = lambda t: t[1])
+                    print('autopick: estimated first arrivals for P: {}, S:{} using TauPy'.format(estFirstP, estFirstS))
+
+                    # modifiy pstart and pstop relative to estimated first P arrival (relative to station time axis)
+                    pstart += (source_origin.time + estFirstP) - wfstart
+                    pstop += (source_origin.time + estFirstP) - wfstart
+                    Lc = pstop - pstart
+                else:
+                    print('No source origins given!')
+
+        else:
+            Lc = pstop - pstart
         Lwf = zdat[0].stats.endtime - zdat[0].stats.starttime
         Ldiff = Lwf - Lc
         if Ldiff < 0:
@@ -238,6 +288,12 @@ def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None):
         else:
             fig = None
         aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP, fig=fig)
+        # add pstart and pstop to aic plot
+        if fig.axes:
+            for ax in fig.axes:
+                ax.vlines(pstart, ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed', label='P start')
+                ax.vlines(pstop, ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed', label='P stop')
+                ax.legend()
         ##############################################################
         if aicpick.getpick() is not None:
             # check signal length to detect spuriously picked noise peaks

pylot/core/util/defaults.py

@@ -9,7 +9,7 @@ Created on Wed Feb 26 12:31:25 2014
 import os
 import platform
 
-from pylot.core.io.inputs import PylotParameter
+from pylot.core.util.utils import readDefaultFilterInformation
 from pylot.core.loc import hypo71
 from pylot.core.loc import hypodd
 from pylot.core.loc import hyposat
@@ -17,23 +17,6 @@ from pylot.core.loc import nll
 from pylot.core.loc import velest
 
 
-def readDefaultFilterInformation(fname):
-    pparam = PylotParameter(fname)
-    return readFilterInformation(pparam)
-
-
-def readFilterInformation(pylot_parameter):
-    p_filter = {'filtertype': pylot_parameter['filter_type'][0],
-                'freq': [pylot_parameter['minfreq'][0], pylot_parameter['maxfreq'][0]],
-                'order': int(pylot_parameter['filter_order'][0])}
-    s_filter = {'filtertype': pylot_parameter['filter_type'][1],
-                'freq': [pylot_parameter['minfreq'][1], pylot_parameter['maxfreq'][1]],
-                'order': int(pylot_parameter['filter_order'][1])}
-    filter_information = {'P': p_filter,
-                          'S': s_filter}
-    return filter_information
-
-
 # determine system dependent path separator
 system_name = platform.system()
 if system_name in ["Linux", "Darwin"]:

pylot/core/util/utils.py

@@ -10,8 +10,8 @@ import subprocess
 import numpy as np
 from obspy import UTCDateTime, read
 from pylot.core.io.inputs import PylotParameter
-from scipy.interpolate import splrep, splev
 
+from scipy.interpolate import splrep, splev
 
 def _pickle_method(m):
     if m.im_self is None:
@@ -20,6 +20,23 @@ def _pickle_method(m):
         return getattr, (m.im_self, m.im_func.func_name)
 
 
+def readDefaultFilterInformation(fname):
+    pparam = PylotParameter(fname)
+    return readFilterInformation(pparam)
+
+
+def readFilterInformation(pylot_parameter):
+    p_filter = {'filtertype': pylot_parameter['filter_type'][0],
+                'freq': [pylot_parameter['minfreq'][0], pylot_parameter['maxfreq'][0]],
+                'order': int(pylot_parameter['filter_order'][0])}
+    s_filter = {'filtertype': pylot_parameter['filter_type'][1],
+                'freq': [pylot_parameter['minfreq'][1], pylot_parameter['maxfreq'][1]],
+                'order': int(pylot_parameter['filter_order'][1])}
+    filter_information = {'P': p_filter,
+                          'S': s_filter}
+    return filter_information
+
+
 def fit_curve(x, y):
     return splev, splrep(x, y)
 
@@ -551,6 +568,51 @@ def which(program, infile=None):
     return None
 
 
+def loopIdentifyPhase(phase):
+    '''
+    Loop through phase string and try to recognize its type (P or S wave).
+    Global variable ALTSUFFIX gives alternative suffix for phases if they do not end with P, p or S, s.
+    If ALTSUFFIX is not given, the function will cut the last letter of the phase string until string ends
+    with P or S.
+    :param phase: phase name (str)
+    :return:
+    '''
+    from pylot.core.util.defaults import ALTSUFFIX
+
+    phase_copy = phase
+    while not identifyPhase(phase_copy):
+        identified = False
+        for alt_suf in ALTSUFFIX:
+            if phase_copy.endswith(alt_suf):
+                phase_copy = phase_copy.split(alt_suf)[0]
+                identified = True
+        if not identified:
+            phase_copy = phase_copy[:-1]
+        if len(phase_copy) < 1:
+            print('Warning: Could not identify phase {}!'.format(phase))
+            return
+    return phase_copy
+
+
+def identifyPhase(phase):
+    '''
+    Returns capital P or S if phase string is identified by last letter. Else returns False.
+    :param phase: phase name (str)
+    :return: 'P', 'S' or False
+    '''
+    # common phase suffix for P and S
+    common_P = ['P', 'p']
+    common_S = ['S', 's']
+    if phase[-1] in common_P:
+        return 'P'
+    if phase[-1] in common_S:
+        return 'S'
+    else:
+        return False
+
+
+
+
 if __name__ == "__main__":
     import doctest
 

pylot/core/util/widgets.py

@@ -21,7 +21,7 @@ except:
     pg = None
 
 from matplotlib.figure import Figure
-from pylot.core.util.utils import find_horizontals
+from pylot.core.util.utils import find_horizontals, identifyPhase, loopIdentifyPhase
 
 try:
     from matplotlib.backends.backend_qt4agg import FigureCanvas
@@ -46,7 +46,7 @@ from pylot.core.io.inputs import FilterOptions, PylotParameter
 from pylot.core.pick.utils import getSNR, earllatepicker, getnoisewin, \
     getResolutionWindow
 from pylot.core.pick.compare import Comparison
-from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS, ALTSUFFIX, \
+from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS, \
     SetChannelComponents
 from pylot.core.util.utils import prepTimeAxis, full_range, scaleWFData, \
     demeanTrace, isSorted, findComboBoxIndex, clims
@@ -118,34 +118,6 @@ def createAction(parent, text, slot=None, shortcut=None, icon=None,
     return action
 
 
-def loopIdentifyPhase(phase):
-    phase_copy = phase
-    while not identifyPhase(phase_copy):
-        identified = False
-        for alt_suf in ALTSUFFIX:
-            if phase_copy.endswith(alt_suf):
-                phase_copy = phase_copy.split(alt_suf)[0]
-                identified = True
-        if not identified:
-            phase_copy = phase_copy[:-1]
-        if len(phase_copy) < 1:
-            print('Warning: Could not identify phase {}!'.format(phase))
-            return
-    return phase_copy
-
-
-def identifyPhase(phase):
-    # common phase suffix for P and S
-    common_P = ['P', 'p']
-    common_S = ['S', 's']
-    if phase[-1] in common_P:
-        return 'P'
-    if phase[-1] in common_S:
-        return 'S'
-    else:
-        return False
-
-
 class ComparisonDialog(QDialog):
     def __init__(self, c, parent=None):
         self._data = c