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

QtPyLoT.py

@@ -43,7 +43,7 @@ from obspy import UTCDateTime
 from pylot.core.io.data import Data
 from pylot.core.io.inputs import FilterOptions, AutoPickParameter
 from pylot.core.pick.autopick import autopickevent
-from pylot.core.io.phases import picks_to_dict
+from pylot.core.io.phases import picksdict_from_picks
 from pylot.core.loc.nll import locate as locateNll
 from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP,\
     AUTOMATIC_DEFAULTS
@@ -735,7 +735,7 @@ class MainWindow(QMainWindow):
         return rval
 
     def updatePicks(self, type='manual'):
-        picks = picks_to_dict(evt=self.getData().getEvtData())
+        picks = picksdict_from_picks(evt=self.getData().getEvtData())
         if type == 'manual':
             self.picks.update(picks)
         elif type == 'auto':

pylot/core/io/data.py

@@ -9,7 +9,7 @@ from obspy.core import read, Stream, UTCDateTime
 from obspy import read_events, read_inventory
 from obspy.core.event import Event, ResourceIdentifier, Pick, WaveformStreamID
 
-from pylot.core.io.phases import readPILOTEvent, picks_from_dict
+from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict
 from pylot.core.util.utils import fnConstructor, getGlobalTimes
 from pylot.core.util.errors import FormatError, OverwriteError
 
@@ -411,19 +411,18 @@ class Data(object):
             :raise OverwriteError: raises an OverwriteError if the picks list is
              not empty. The GUI will then ask for a decision.
             """
-            firstonset = None
+
+            #firstonset = find_firstonset(picks)
             if self.getEvtData().picks:
                 raise OverwriteError('Actual picks would be overwritten!')
-            picks, firstonset = picks_from_dict(picks)
+            picks = picks_from_picksdict(picks)
             self.getEvtData().picks = picks
+            # if 'smi:local' in self.getID() and firstonset:
+            #     fonset_str = firstonset.strftime('%Y_%m_%d_%H_%M_%S')
+            #     ID = ResourceIdentifier('event/' + fonset_str)
+            #     ID.convertIDToQuakeMLURI(authority_id=authority_id)
+            #     self.getEvtData().resource_id = ID
 
-            if 'smi:local' in self.getID() and firstonset:
-                fonset_str = firstonset.strftime('%Y_%m_%d_%H_%M_%S')
-                ID = ResourceIdentifier('event/' + fonset_str)
-                ID.convertIDToQuakeMLURI(authority_id=authority_id)
-                self.getEvtData().resource_id = ID
-            else:
-                print('No picks to apply!')
 
         def applyArrivals(arrivals):
             """

pylot/core/io/phases.py

@@ -8,9 +8,9 @@ import scipy.io as sio
 import obspy.core.event as ope
 from obspy.core import UTCDateTime
 
+from pylot.core.pick.utils import select_for_phase
 from pylot.core.util.utils import getOwner, createPick, createArrival, \
-    createEvent, createOrigin, createMagnitude
+    createEvent, createOrigin, createMagnitude, getGlobalTimes
-
 
 def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs):
     """
@@ -192,14 +192,14 @@ def picksdict_from_obs(fn):
     return picks
 
 
-def picks_to_dict(evt):
+def picksdict_from_picks(evt):
-    '''
+    """
     Takes an Event object and return the pick dictionary commonly used within
     PyLoT
     :param evt: Event object contain all available information
     :type evt: `~obspy.core.event.Event`
     :return: pick dictionary
-    '''
+    """
     picks = {}
     for pick in evt.picks:
         phase = {}
@@ -229,12 +229,12 @@ def picks_to_dict(evt):
         picks[station] = onsets.copy()
     return picks
 
-def picks_from_dict(picks):
+def picks_from_picksdict(picks):
-    firstonset = None
+    picks_list = list()
     for station, onsets in picks.items():
         print('Reading picks on station %s' % station)
         for label, phase in onsets.items():
-            if not isinstance(phase, dict):
+            if not isinstance(phase, dict) or len(phase) < 3:
                 continue
             onset = phase['mpp']
             epp = phase['epp']
@@ -262,52 +262,90 @@ def picks_from_dict(picks):
                 else:
                     pick.polarity = 'undecidable'
             except KeyError as e:
-                print('No polarity information found for %s' % phase)
+                print(e.message, 'No polarity information found for %s' % phase)
-            if firstonset is None or firstonset > onset:
+            picks_list.append(pick)
-                firstonset = onset
+    return picks_list
 
 
-def reassess_pilot_event(root_dir, event_id):
+def reassess_pilot_db(root_dir, out_dir=None, fn_param=None):
+    import glob
+
+    evt_list = glob.glob1(root_dir,'e????.???.??')
+
+    for evt in evt_list:
+        reassess_pilot_event(root_dir, evt, out_dir, fn_param)
+
+
+def reassess_pilot_event(root_dir, event_id, out_dir=None, fn_param=None, verbosity=0):
     from obspy import read
-    from pylot.core.util.defaults import AUTOMATIC_DEFAULTS
+
     from pylot.core.io.inputs import AutoPickParameter
     from pylot.core.pick.utils import earllatepicker
 
-    default = AutoPickParameter(AUTOMATIC_DEFAULTS)
+    if fn_param is None:
+        import pylot.core.util.defaults as defaults
+        fn_param = defaults.AUTOMATIC_DEFAULTS
+
+    default = AutoPickParameter(fn_param)
 
     search_base = os.path.join(root_dir, event_id)
     phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat'))
-    picks_dict = picks_from_pilot(phases_file)
+    if not phases_file:
+        return
+    print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0]))
+    picks_dict = picks_from_pilot(phases_file[0])
+    print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict))
     for station in picks_dict.keys():
         fn_pattern = os.path.join(search_base, '{0}*'.format(station))
         try:
-            st = read(fn_pattern)
+            st = read(fn_pattern, format='GSE2')
         except TypeError as e:
             print(e.message)
-            st = read(fn_pattern, format='GSE2')
+            st = read(fn_pattern)
-        if not st:
-            raise RuntimeError('no waveform data found for station {station}'.format(station=station))
         for phase in picks_dict[station].keys():
             try:
                 mpp = picks_dict[station][phase]['mpp']
             except KeyError as e:
                 print(e.message, station)
                 continue
-            epp, lpp, spe = earllatepicker(st,
+            sel_st = select_for_phase(st, phase)
+            if not sel_st:
+                raise warnings.formatwarning(
+                    'no waveform data found for station {station}'.format(
+                        station=station), category=RuntimeWarning)
+            print(sel_st)
+            stime, etime = getGlobalTimes(sel_st)
+            rel_pick = mpp - stime
+            epp, lpp, spe = earllatepicker(sel_st,
                                            default.get('nfac{0}'.format(phase)),
                                            default.get('tsnrz' if phase == 'P' else 'tsnrh'),
-                                           mpp
+                                           Pick1=rel_pick,
+                                           iplot=None,
                                            )
+            if epp is None or lpp is None:
+                continue
+            epp = stime + epp
+            lpp = stime + lpp
+            min_diff = 3 * st[0].stats.delta
+            if lpp - mpp < min_diff:
+                lpp = mpp + min_diff
+            if mpp - epp < min_diff:
+                epp = mpp - min_diff
             picks_dict[station][phase] = dict(epp=epp, mpp=mpp, lpp=lpp, spe=spe)
     # create Event object for export
     evt = ope.Event(resource_id=event_id)
-    evt.picks = picks_from_dict(picks_dict)
+    evt.picks = picks_from_picksdict(picks_dict)
     # write phase information to file
-    evt.write('{0}.xml'.format(event_id), format='QUAKEML')
+    if not out_dir:
+        fnout_prefix = os.path.join(root_dir, event_id, '{0}.'.format(event_id))
+    else:
+        fnout_prefix = os.path.join(out_dir, '{0}.'.format(event_id))
+    evt.write(fnout_prefix + 'xml', format='QUAKEML')
+    #evt.write(fnout_prefix + 'cnv', format='VELEST')
 
 
 def writephases(arrivals, fformat, filename):
-    '''
+    """
     Function of methods to write phases to the following standard file
     formats used for locating earthquakes:
 
@@ -325,7 +363,7 @@ def writephases(arrivals, fformat, filename):
 
     :param: filename, full path and name of phase file
     :type: string
-    '''
+    """
 
     if fformat == 'NLLoc':
         print ("Writing phases to %s for NLLoc" % filename)
@@ -387,7 +425,6 @@ def writephases(arrivals, fformat, filename):
                                                                                            sweight))
 
         fid.close()
-
     elif fformat == 'HYPO71':
         print ("Writing phases to %s for HYPO71" % filename)
         fid = open("%s" % filename, 'w')

pylot/core/pick/compare.py

@@ -9,7 +9,7 @@ import matplotlib.pyplot as plt
 
 from obspy import read_events
 
-from pylot.core.io.phases import picks_to_dict
+from pylot.core.io.phases import picksdict_from_picks
 from pylot.core.util.pdf import ProbabilityDensityFunction
 from pylot.core.util.version import get_git_version as _getVersionString
 
@@ -86,8 +86,8 @@ class Comparison(object):
         """
         compare_pdfs = dict()
 
-        pdf_a = self.get(self.names[0]).pdf_data(type)
+        pdf_a = self.get(self.names[0]).generate_pdf_data(type)
-        pdf_b = self.get(self.names[1]).pdf_data(type)
+        pdf_b = self.get(self.names[1]).generate_pdf_data(type)
 
         for station, phases in pdf_a.items():
             if station in pdf_b.keys():
@@ -206,6 +206,7 @@ class PDFDictionary(object):
 
     def __init__(self, data):
         self._pickdata = data
+        self._pdfdata = self.generate_pdf_data()
 
     def __nonzero__(self):
         if len(self.pick_data) < 1:
@@ -213,6 +214,14 @@ class PDFDictionary(object):
         else:
             return True
 
+    @property
+    def pdf_data(self):
+        return self._pdfdata
+
+    @pdf_data.setter
+    def pdf_data(self, data):
+        self._pdfdata = data
+
     @property
     def pick_data(self):
         return self._pickdata
@@ -221,15 +230,23 @@ class PDFDictionary(object):
     def pick_data(self, data):
         self._pickdata = data
 
+    @property
+    def stations(self):
+        return self.pick_data.keys()
+
+    @property
+    def nstations(self):
+        return len(self.stations)
+
     @classmethod
     def from_quakeml(self, fn):
         cat = read_events(fn)
         if len(cat) > 1:
             raise NotImplementedError('reading more than one event at the same '
                                       'time is not implemented yet! Sorry!')
-        return PDFDictionary(picks_to_dict(cat[0]))
+        return PDFDictionary(picksdict_from_picks(cat[0]))
 
-    def pdf_data(self, type='exp'):
+    def generate_pdf_data(self, type='exp'):
         """
         Returns probabiliy density function dictionary containing the
         representation of the actual pick_data.
@@ -251,3 +268,57 @@ class PDFDictionary(object):
 
         return pdf_picks
 
+    def plot(self, stations=None):
+        '''
+        plots the all probability density function for either desired STATIONS
+        or all available date
+        :param stations: list of stations to be plotted
+        :type stations: list
+        :return: matplotlib figure object containing the plot
+        '''
+        assert stations is not None or not isinstance(stations, list), \
+            'parameter stations should be a list not {0}'.format(type(stations))
+        if not stations:
+            nstations = self.nstations
+            stations = self.stations
+        else:
+            nstations = len(stations)
+
+        istations = range(nstations)
+        fig, axarr = plt.subplots(nstations, 2, sharex='col', sharey='row')
+        hide_labels = True
+
+        for n in istations:
+            station = stations[n]
+            pdfs = self.pdf_data[station]
+            for l, phase in enumerate(pdfs.keys()):
+                try:
+                    axarr[n, l].plot(pdfs[phase].axis, pdfs[phase].data)
+                    if n is 0:
+                        axarr[n, l].set_title(phase)
+                    if l is 0:
+                        axann = axarr[n, l].annotate(station, xy=(.05, .5),
+                                                     xycoords='axes fraction')
+                        bbox_props = dict(boxstyle='round', facecolor='lightgrey',
+                                          alpha=.7)
+                        axann.set_bbox(bbox_props)
+                    if n == int(np.median(istations)) and l is 0:
+                        label = 'probability density (qualitative)'
+                        axarr[n, l].set_ylabel(label)
+                except IndexError as e:
+                    print('trying aligned plotting\n{0}'.format(e))
+                    hide_labels = False
+                    axarr[l].plot(pdfs[phase].axis, pdfs[phase].data)
+                    axarr[l].set_title(phase)
+                    if l is 0:
+                        axann = axarr[l].annotate(station, xy=(.05, .5),
+                                                     xycoords='axes fraction')
+                        bbox_props = dict(boxstyle='round', facecolor='lightgrey',
+                                          alpha=.7)
+                        axann.set_bbox(bbox_props)
+        if hide_labels:
+            plt.setp([a.get_xticklabels() for a in axarr[0, :]], visible=False)
+            plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False)
+            plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False)
+
+        return fig

pylot/core/pick/utils.py

@@ -464,6 +464,37 @@ def getResolutionWindow(snr):
     return time_resolution / 2
 
 
+def select_for_phase(st, phase):
+    '''
+    takes a STream object and a phase name and returns that particular component
+    which presumably shows the chosen PHASE best
+
+    :param st: stream object containing one or more component[s]
+    :type st: `~obspy.core.stream.Stream`
+    :param phase: label of the phase for which the stream selection is carried
+        out; 'P' or 'S'
+    :type phase: str
+    :return:
+    '''
+    from pylot.core.util.defaults import COMPNAME_MAP
+
+    sel_st = Stream()
+    if phase.upper() == 'P':
+        comp = 'Z'
+        alter_comp = COMPNAME_MAP[comp]
+        sel_st += st.select(component=comp)
+        sel_st += st.select(component=alter_comp)
+    elif phase.upper() == 'S':
+        comps = 'NE'
+        for comp in comps:
+            alter_comp = COMPNAME_MAP[comp]
+            sel_st += st.select(component=comp)
+            sel_st += st.select(component=alter_comp)
+    else:
+        raise TypeError('Unknown phase label: {0}'.format(phase))
+    return sel_st
+
+
 def wadaticheck(pickdic, dttolerance, iplot):
     '''
     Function to calculate Wadati-diagram from given P and S onsets in order
@@ -938,10 +969,6 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
     return returnflag
 
 
-def reassess_pilot_event():
-    pass
-
-
 if __name__ == '__main__':
     import doctest
 

pylot/core/util/dataprocessing.py

@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import os
+import glob
+from obspy import UTCDateTime
+
+def gse_time_header(lines):
+    '''
+    takes a path FILE to a GSE data file and returns the time header of the file
+    :param file: path to GSE data file
+    :type file: str
+    :return: time header from FILE
+    :rtype: str
+
+    >>> gse_time_header('test.gse')
+    "WID2 2005/10/09 20:17:25.000 ATH   SHZ  NSP CM6     9000   50.000000   0.10E+01   1.000    NSP  -1.0  0.0"
+    '''
+
+    return lines[1]
+
+def time_from_header(header):
+    timeline = header.split(' ')
+    time = timeline[1].split('/') + timeline[2].split(':')
+    time = time[:-1] + time[-1].split('.')
+    time[-1] += '000'
+    return [int(t) for t in time]
+
+def check_time(time):
+    try:
+        UTCDateTime(time)
+        return True
+    except ValueError:
+        return False

pylot/core/util/pdf.py

@@ -256,9 +256,10 @@ class ProbabilityDensityFunction(object):
         '''
 
         rval = 0
-        for n, x in enumerate(self.axis):
+        axis = self.axis - self.x0
+        for n, x in enumerate(axis):
             rval += x * self.data[n]
-        return rval * self.incr
+        return rval * self.incr + self.x0
 
     def standard_deviation(self):
         mu = self.expectation()

scripts/pylot-reasses-pilot-db.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import argparse
+
+from pylot.core.util.version import get_git_version as _getVersionString
+from pylot.core.io.phases import reassess_pilot_db
+
+__version__ = _getVersionString()
+__author__ = 'S. Wehling-Benatelli'
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(
+        description='reassess old PILOT event data base in terms of consistent '
+                    'automatic uncertainty estimation',
+        epilog='Script written by {author} belonging to PyLoT version'
+               ' {version}\n'.format(author=__author__,
+                                     version=__version__)
+    )
+
+    parser.add_argument(
+        'dbroot', type=str, help='specifies the root directory (in '
+                                            'most cases PILOT database folder)'
+    )
+    parser.add_argument(
+        '--output', '-o', type=str, help='path to the output directory', dest='output'
+    )
+    parser.add_argument(
+        '--parameterfile', '-p', type=str, help='full path to the parameterfile', dest='parfile'
+    )
+
+    args = parser.parse_args()
+    reassess_pilot_db(args.dbroot, args.output, args.parfile)

scripts/pylot-reasses-pilot-event.py

@@ -3,23 +3,34 @@
 
 import argparse
 
-from pylot.core.pick.utils import reassess_pilot_event
 from pylot.core.util.version import get_git_version as _getVersionString
 from pylot.core.io.phases import reassess_pilot_event
 
 __version__ = _getVersionString()
-__author__ = 'sebastianw'
+__author__ = 'S. Wehling-Benatelli'
 
 if __name__ == '__main__':
-    parser = argparse.ArgumentParser()
+    parser = argparse.ArgumentParser(
+        description='reassess old PILOT event data in terms of consistent '
+                    'automatic uncertainty estimation',
+        epilog='Script written by {author} belonging to PyLoT version'
+               ' {version}\n'.format(author=__author__,
+                                     version=__version__)
+    )
 
     parser.add_argument(
-        '--directory', '-d', type=str, help='specifies the root directory (in '
+        'dbroot', type=str, help='specifies the root directory (in '
                                             'most cases PILOT database folder)'
     )
     parser.add_argument(
-        '--eventid', '-i', type=str, help='PILOT event identifier'
+        'id', type=str, help='PILOT event identifier'
+    )
+    parser.add_argument(
+        '--output', '-o', type=str, help='path to the output directory', dest='output'
+    )
+    parser.add_argument(
+        '--parameterfile', '-p', type=str, help='full path to the parameterfile', dest='parfile'
     )
 
     args = parser.parse_args()
-    reassess_pilot_event(args.dir, args.id)
+    reassess_pilot_event(args.dbroot, args.id, args.output, args.parfile)