Merge branch 'feature/magnitude4QtPyLoT' into develop

Conflicts:
	pylot/core/analysis/magnitude.py
	pylot/core/util/widgets.py

QtPyLoT.py

@@ -23,10 +23,11 @@ https://www.iconfinder.com/iconsets/flavour
     (http://www.gnu.org/copyleft/lesser.html)
 """
 
-import matplotlib
 import os
 import sys
 
+import matplotlib
+
 matplotlib.use('Qt4Agg')
 matplotlib.rcParams['backend.qt4'] = 'PySide'
 
@@ -37,24 +38,22 @@ from PySide.QtGui import QMainWindow, QInputDialog, QIcon, QFileDialog, \
     QDialog, QErrorMessage, QApplication, QPixmap, QMessageBox, QSplashScreen, \
     QActionGroup, QListWidget, QDockWidget, QLineEdit
 import numpy as np
-import subprocess
 from obspy import UTCDateTime
-from obspy.geodetics import degrees2kilometers
-from obspy.core.event import Magnitude
 
-from pylot.core.analysis.magnitude import calcsourcespec, calcMoMw
+from pylot.core.analysis.magnitude import RichterMagnitude, MomentMagnitude
 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.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, COMPNAME_MAP, \
     AUTOMATIC_DEFAULTS
 from pylot.core.util.errors import FormatError, DatastructureError, \
-    OverwriteError
+    OverwriteError, ProcessingError
 from pylot.core.util.connection import checkurl
-from pylot.core.util.dataprocessing import read_metadata
+from pylot.core.util.dataprocessing import read_metadata, restitute_data
 from pylot.core.util.utils import fnConstructor, getLogin, \
     full_range
 from pylot.core.io.location import create_creation_info, create_event
@@ -124,6 +123,11 @@ class MainWindow(QMainWindow):
         # load and display waveform data
         self.dirty = False
         self.load_data()
+        finv = settings.value("inventoryFile", None)
+        if finv is not None:
+            self._metadata = read_metadata(finv)
+        else:
+            self._metadata = None
         if self.loadWaveformData():
             self.updateFilterOptions()
         else:
@@ -368,6 +372,17 @@ class MainWindow(QMainWindow):
 
         self.setCentralWidget(_widget)
 
+
+    @property
+    def metadata(self):
+        return self._metadata
+
+
+    @metadata.setter
+    def metadata(self, value):
+        self._metadata = value
+
+
     def updateFileMenu(self):
 
         self.fileMenu.clear()
@@ -911,6 +926,8 @@ class MainWindow(QMainWindow):
             epp = picks['epp'] - stime
             lpp = picks['lpp'] - stime
             spe = picks['spe']
+            if not spe:
+                spe = symmetrize_error(mpp - epp, lpp - mpp)
 
             if picktype == 'manual':
                 ax.fill_between([epp, lpp], ylims[0], ylims[1],
@@ -984,53 +1001,55 @@ class MainWindow(QMainWindow):
             os.remove(phasepath)
 
         self.get_data().applyEVTData(lt.read_location(locpath), type='event')
-        self.get_data().get_evt_data().magnitudes.append(self.calc_magnitude())
+        self.get_data().applyEVTData(self.calc_magnitude(), type='event')
 
-    def calc_magnitude(self):
+
-        e = self.get_data().get_evt_data()
+    def calc_magnitude(self, type='ML'):
+        def set_inv(settings):
+            fninv, _ = QFileDialog.getOpenFileName(self, self.tr(
+                "Select inventory..."), self.tr("Select file"))
+            if not fninv:
+                return False
+            ans = QMessageBox.question(self, self.tr("Make default..."),
+                                       self.tr(
+                                           "New inventory filename set.\n" + \
+                                           "Do you want to make it the default value?"),
+                                       QMessageBox.Yes | QMessageBox.No,
+                                       QMessageBox.No)
+            if ans == QMessageBox.Yes:
+                settings.setValue("inventoryFile", fninv)
+                settings.sync()
+            self.metadata = read_metadata(fninv)
+            return True
+        
         settings = QSettings()
-        if e.origins:
+        fninv = settings.value("inventoryFile", None)
-            o = e.origins[0]
-            mags = dict()
-            fninv = settings.value("inventoryFile", None)
-            if fninv is None:
-                fninv, _ = QFileDialog.getOpenFileName(self, self.tr(
-                    "Select inventory..."), self.tr("Select file"))
-                ans = QMessageBox.question(self, self.tr("Make default..."),
-                                           self.tr(
-                                               "New inventory filename set.\n" + \
-                                               "Do you want to make it the default value?"),
-                                           QMessageBox.Yes | QMessageBox.No,
-                                           QMessageBox.No)
-                if ans == QMessageBox.Yes:
-                    settings.setValue("inventoryFile", fninv)
-                    settings.sync()
-            metadata = read_metadata(fninv)
-            for a in o.arrivals:
-                if a.phase in 'sS':
-                    continue
-                pick = a.pick_id.get_referred_object()
-                station = pick.waveform_id.station_code
-                wf = self.get_data().getWFData().select(station=station)
-                if not wf:
-                    continue
-                onset = pick.time
-                dist = degrees2kilometers(a.distance)
-                w0, fc = calcsourcespec(wf, onset, metadata, self.inputs.get('vp'), dist,
-                                        a.azimuth, a.takeoff_angle,
-                                        self.inputs.get('Qp'), 0)
-                if w0 is None or fc is None:
-                    continue
-                station_mag = calcMoMw(wf, w0, self.inputs.get('rho'),
-                                       self.inputs.get('vp'), dist)
-                mags[station] = station_mag
-            mag = np.median([M[1] for M in mags.values()])
-            # give some information on the processing
-            print('number of stations used: {0}\n'.format(len(mags.values())))
-            print('stations used:\n')
-            for s in mags.keys(): print('\t{0}'.format(s))
 
-            return Magnitude(mag=mag, magnitude_type='Mw')
+        if fninv is None and not self.metadata:
+            if not set_inv(settings):
+                return None
+        elif fninv is not None and not self.metadata:
+            ans = QMessageBox.question(self, self.tr("Use default..."),
+                                       self.tr(
+                                           "Do you want to use the default value?"),
+                                       QMessageBox.Yes | QMessageBox.No,
+                                       QMessageBox.Yes)
+            if ans == QMessageBox.No:
+                if not set_inv(settings):
+                    return None
+            else:
+                self.metadata = read_metadata(fninv)
+                
+        wf_copy = self.get_data().getWFData().copy()
+        [corr_wf, rest_flag] = restitute_data(wf_copy, *self.metadata)
+        if not rest_flag:
+            raise ProcessingError('Restitution of waveform data failed!')
+        if type == 'ML':
+            local_mag = RichterMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('sstop'), verbosity = True)
+            return local_mag.updated_event()
+        elif type == 'Mw':
+            moment_mag = MomentMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho'), verbosity = True)
+            return moment_mag.updated_event()
         else:
             return None
 

autoPyLoT.py

@@ -5,17 +5,18 @@ from __future__ import print_function
 
 import argparse
 import glob
-import string
 import os
 
-import numpy as np
+from obspy import read_events
 
-from pylot.core.analysis.magnitude import M0Mw
+import pylot.core.loc.hsat as hsat
+import pylot.core.loc.nll as nll
+from pylot.core.analysis.magnitude import MomentMagnitude, RichterMagnitude
 from pylot.core.io.data import Data
 from pylot.core.io.inputs import AutoPickParameter
-import pylot.core.loc.nll as nll
-import pylot.core.loc.hsat as hsat
 from pylot.core.pick.autopick import autopickevent, iteratepicker
+from pylot.core.util.dataprocessing import restitute_data, read_metadata, \
+    remove_underscores
 from pylot.core.util.structure import DATASTRUCTURE
 from pylot.core.util.version import get_git_version as _getVersionString
 
@@ -54,6 +55,7 @@ def autoPyLoT(inputfile):
 
     data = Data()
 
+    evt = None
     # getting information on data structure
 
     if parameter.hasParam('datastructure'):
@@ -97,139 +99,26 @@ def autoPyLoT(inputfile):
             print("!!No source parameter estimation possible!!")
             print("                 !!!              ")
 
-        # multiple event processing
-        # read each event in database
         datapath = datastructure.expandDataPath()
         if not parameter.hasParam('eventID'):
-            for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]:
+            # multiple event processing
-                data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
+            # read each event in database
-                print('Working on event %s' % event)
+            events = [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]
-                print(data)
-                wfdat = data.getWFData()  # all available streams
-                ##########################################################
-                # !automated picking starts here!
-                picks = autopickevent(wfdat, parameter)
-                finalpicks = picks
-                ##########################################################
-                # locating
-                if locflag == 1:
-                    # write phases to NLLoc-phase file
-                    nll.export(picks, phasefile)
-
-                    # For locating the event the NLLoc-control file has to be modified!
-                    evID = event[string.rfind(event, "/") + 1: len(events) - 1]
-                    nllocout = '%s_%s' % (evID, nllocoutpatter)
-                    # create comment line for NLLoc-control file
-                    nll.modify_inputs(ctrf, nllocroot, nllocout, phasef,
-                                      ttpat)
-
-                    # locate the event
-                    nll.locate(ctrfile)
-
-                    # !iterative picking if traces remained unpicked or occupied with bad picks!
-                    # get theoretical onset times for picks with weights >= 4
-                    # in order to reprocess them using smaller time windows around theoretical onset
-                    # get stations with bad onsets
-                    badpicks = []
-                    for key in picks:
-                        if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
-                            badpicks.append([key, picks[key]['P']['mpp']])
-
-                    if len(badpicks) == 0:
-                        print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!")
-                        # get NLLoc-location file
-                        locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
-                        if len(glob.glob(locsearch)) > 0:
-                            # get latest NLLoc-location file if several are available
-                            nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
-                            # calculating seismic moment Mo and moment magnitude Mw
-                            finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
-                                              nllocfile, picks, parameter.get('rho'), \
-                                              parameter.get('vp'), parameter.get('Qp'), \
-                                              parameter.get('invdir'))
-                        else:
-                            print("autoPyLoT: No NLLoc-location file available!")
-                            print("No source parameter estimation possible!")
-                    else:
-                        # get theoretical P-onset times from NLLoc-location file
-                        locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
-                        if len(glob.glob(locsearch)) > 0:
-                            # get latest file if several are available
-                            nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
-                            nlloccounter = 0
-                            while len(badpicks) > 0 and nlloccounter <= maxnumit:
-                                nlloccounter += 1
-                                if nlloccounter > maxnumit:
-                                    print("autoPyLoT: Number of maximum iterations reached, stop iterative picking!")
-                                    break
-                                print("autoPyLoT: Starting with iteration No. %d ..." % nlloccounter)
-                                picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
-                                # write phases to NLLoc-phase file
-                                nll.export(picks, phasefile)
-                                # remove actual NLLoc-location file to keep only the last
-                                os.remove(nllocfile)
-                                # locate the event
-                                nll.locate(ctrfile)
-                                print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
-                                # get updated NLLoc-location file
-                                nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
-                                # check for bad picks
-                                badpicks = []
-                                for key in picks:
-                                    if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
-                                        badpicks.append([key, picks[key]['P']['mpp']])
-                                print("autoPyLoT: After iteration No. %d: %d bad onsets found ..." % (nlloccounter, \
-                                                                                                      len(badpicks)))
-                                if len(badpicks) == 0:
-                                    print("autoPyLoT: No more bad onsets found, stop iterative picking!")
-                                    nlloccounter = maxnumit
-
-                            # calculating seismic moment Mo and moment magnitude Mw
-                            finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
-                                              nllocfile, picks, parameter.get('rho'), \
-                                              parameter.get('vp'), parameter.get('Qp'), \
-                                              parameter.get('invdir'))
-                            # get network moment magntiude
-                            netMw = []
-                            for key in finalpicks.getpicdic():
-                                if finalpicks.getpicdic()[key]['P']['Mw'] is not None:
-                                    netMw.append(finalpicks.getpicdic()[key]['P']['Mw'])
-                            netMw = np.median(netMw)
-                            print("Network moment magnitude: %4.1f" % netMw)
-                        else:
-                            print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
-                ##########################################################
-                # write phase files for various location routines
-                # HYPO71
-                hypo71file = '%s/autoPyLoT_HYPO71.pha' % event
-                if hasattr(finalpicks, 'getpicdic') and finalpicks.getpicdic() is not None:
-                    hsat.export(finalpicks.getpicdic(), hypo71file)
-                    data.applyEVTData(finalpicks.getpicdic())
-                else:
-                    hsat.export(picks, hypo71file)
-                    data.applyEVTData(picks)
-                fnqml = '%s/autoPyLoT' % event
-                data.exportEvent(fnqml)
-
-                endsplash = '''------------------------------------------\n'
-                               -----Finished event %s!-----\n'
-                               ------------------------------------------'''.format \
-                                (version=_getVersionString()) % evID
-                print(endsplash)
-                if locflag == 0:
-                    print("autoPyLoT was running in non-location mode!")
-
-        # single event processing
         else:
-            data.setWFData(glob.glob(os.path.join(datapath, parameter.get('eventID'), '*')))
+            # single event processing
-            print("Working on event {0}".format(parameter.get('eventID')))
+            events = glob.glob(os.path.join(datapath, parameter.get('eventID')))
+        for event in events:
+            data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
+            evID = os.path.split(event)[-1]
+            print('Working on event %s' % event)
             print(data)
-
             wfdat = data.getWFData()  # all available streams
+            wfdat = remove_underscores(wfdat)
+            metadata =  read_metadata(parameter.get('invdir'))
+            corr_dat, rest_flag = restitute_data(wfdat.copy(), *metadata)
             ##########################################################
             # !automated picking starts here!
             picks = autopickevent(wfdat, parameter)
-            finalpicks = picks
             ##########################################################
             # locating
             if locflag == 1:
@@ -237,12 +126,14 @@ def autoPyLoT(inputfile):
                 nll.export(picks, phasefile)
 
                 # For locating the event the NLLoc-control file has to be modified!
-                nllocout = '%s_%s' % (parameter.get('eventID'), nllocoutpatter)
+                nllocout = '%s_%s' % (evID, nllocoutpatter)
                 # create comment line for NLLoc-control file
-                nll.modify_inputs(ctrf, nllocroot, nllocout, phasef, ttpat)
+                nll.modify_inputs(ctrf, nllocroot, nllocout, phasef,
+                                  ttpat)
 
                 # locate the event
                 nll.locate(ctrfile)
+
                 # !iterative picking if traces remained unpicked or occupied with bad picks!
                 # get theoretical onset times for picks with weights >= 4
                 # in order to reprocess them using smaller time windows around theoretical onset
@@ -252,18 +143,30 @@ def autoPyLoT(inputfile):
                     if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
                         badpicks.append([key, picks[key]['P']['mpp']])
 
+                # TODO keep code DRY (Don't Repeat Yourself) the following part is written twice
+                # suggestion: delete block and modify the later similar block to work properly
+
                 if len(badpicks) == 0:
                     print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!")
                     # get NLLoc-location file
                     locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
                     if len(glob.glob(locsearch)) > 0:
-                        # get latest NLLOc-location file if several are available
+                        # get latest NLLoc-location file if several are available
                         nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
+                        evt = read_events(nllocfile)[0]
                         # calculating seismic moment Mo and moment magnitude Mw
-                        finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
+                        moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'),
-                                          nllocfile, picks, parameter.get('rho'), \
+                                                     parameter.get('Qp'),
-                                          parameter.get('vp'), parameter.get('Qp'), \
+                                                     parameter.get('rho'), True, 0)
-                                          parameter.get('invdir'))
+                        # update pick with moment property values (w0, fc, Mo)
+                        for station, props in moment_mag.moment_props.items():
+                            picks[station]['P'].update(props)
+                        evt = moment_mag.updated_event()
+                        local_mag = RichterMagnitude(corr_dat, evt,
+                                                     parameter.get('sstop'), True, 0)
+                        for station, amplitude in local_mag.amplitudes.items():
+                            picks[station]['S']['Ao'] = amplitude.generic_amplitude
+                        evt = local_mag.updated_event()
                     else:
                         print("autoPyLoT: No NLLoc-location file available!")
                         print("No source parameter estimation possible!")
@@ -300,38 +203,39 @@ def autoPyLoT(inputfile):
                             if len(badpicks) == 0:
                                 print("autoPyLoT: No more bad onsets found, stop iterative picking!")
                                 nlloccounter = maxnumit
-
+                        evt = read_events(nllocfile)[0]
                         # calculating seismic moment Mo and moment magnitude Mw
-                        finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
+                        moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'),
-                                          nllocfile, picks, parameter.get('rho'), \
+                                                     parameter.get('Qp'),
-                                          parameter.get('vp'), parameter.get('Qp'), \
+                                                     parameter.get('rho'), True, 0)
-                                          parameter.get('invdir'))
+                        # update pick with moment property values (w0, fc, Mo)
-                        # get network moment magntiude
+                        for station, props in moment_mag.moment_props.items():
-                        netMw = []
+                            picks[station]['P'].update(props)
-                        for key in finalpicks.getpicdic():
+                        evt = moment_mag.updated_event()
-                            if finalpicks.getpicdic()[key]['P']['Mw'] is not None:
+                        local_mag = RichterMagnitude(corr_dat, evt,
-                                netMw.append(finalpicks.getpicdic()[key]['P']['Mw'])
+                                                     parameter.get('sstop'), True, 0)
-                        netMw = np.median(netMw)
+                        for station, amplitude in local_mag.amplitudes.items():
-                        print("Network moment magnitude: %4.1f" % netMw)
+                            picks[station]['S']['Ao'] = amplitude.generic_amplitude
+                        evt = local_mag.updated_event()
+                        net_mw = moment_mag.net_magnitude()
+                        print("Network moment magnitude: %4.1f" % net_mw.mag)
                     else:
                         print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
             ##########################################################
             # write phase files for various location routines
             # HYPO71
-            hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.get('eventID'))
+            hypo71file = '%s/autoPyLoT_HYPO71.pha' % event
-            if hasattr(finalpicks, 'getpicdic') and finalpicks.getpicdic() is not None:
+            hsat.export(picks, hypo71file)
-                hsat.export(finalpicks.getpicdic(), hypo71file)
+            data.applyEVTData(picks)
-                data.applyEVTData(finalpicks.getpicdic())
+            if evt is not None:
-            else:
+                data.applyEVTData(evt, 'event')
-                hsat.export(picks, hypo71file)
+            fnqml = '%s/autoPyLoT' % event
-                data.applyEVTData(picks)
-            fnqml =  '%s/%s/autoPyLoT' % (datapath, parameter.get('eventID'))
             data.exportEvent(fnqml)
 
             endsplash = '''------------------------------------------\n'
                            -----Finished event %s!-----\n'
                            ------------------------------------------'''.format \
-                            (version=_getVersionString()) % parameter.get('eventID')
+                            (version=_getVersionString()) % evID
             print(endsplash)
             if locflag == 0:
                 print("autoPyLoT was running in non-location mode!")

pylot/RELEASE-VERSION

@@ -1 +1 @@
-55a58-dirty
+5f92-dirty

pylot/core/active/seismicArrayPreparation.py

@@ -3,6 +3,41 @@ import sys
 import numpy as np
 from scipy.interpolate import griddata
 
+def readMygridNlayers(filename):
+    infile = open(filename, 'r')
+    nlayers = len(infile.readlines()) / 2
+    infile.close()
+
+    return nlayers
+
+def readMygrid(filename):
+    ztop = [];
+    zbot = [];
+    vtop = [];
+    vbot = []
+    infile = open(filename, 'r')
+    nlayers = readMygridNlayers(filename)
+
+    print('\nreadMygrid: Reading file %s.' % filename)
+    for index in range(nlayers):
+        line1 = infile.readline()
+        line2 = infile.readline()
+        ztop.append(float(line1.split()[0]))
+        vtop.append(float(line1.split()[1]))
+        zbot.append(float(line2.split()[0]))
+        vbot.append(float(line2.split()[1]))
+        print('Layer %s:\n[Top: v = %s [km/s], z = %s [m]]'
+              '\n[Bot: v = %s [km/s], z = %s [m]]'
+              % (index + 1, vtop[index], ztop[index],
+                 vbot[index], zbot[index]))
+
+    if not ztop[0] == 0:
+        print('ERROR: there must be a velocity set for z = 0 in the file %s' % filename)
+        print('e.g.:\n0    0.33\n-5    1.0\netc.')
+
+    infile.close()
+    return ztop, zbot, vtop, vbot
+
 
 class SeisArray(object):
     '''
@@ -714,41 +749,6 @@ class SeisArray(object):
         #     rad = angle / 180 * PI
         #     return rad
 
-        def readMygridNlayers(filename):
-            infile = open(filename, 'r')
-            nlayers = len(infile.readlines()) / 2
-            infile.close()
-
-            return nlayers
-
-        def readMygrid(filename):
-            ztop = [];
-            zbot = [];
-            vtop = [];
-            vbot = []
-            infile = open(filename, 'r')
-            nlayers = readMygridNlayers(filename)
-
-            print('\nreadMygrid: Reading file %s.' % filename)
-            for index in range(nlayers):
-                line1 = infile.readline()
-                line2 = infile.readline()
-                ztop.append(float(line1.split()[0]))
-                vtop.append(float(line1.split()[1]))
-                zbot.append(float(line2.split()[0]))
-                vbot.append(float(line2.split()[1]))
-                print('Layer %s:\n[Top: v = %s [km/s], z = %s [m]]'
-                      '\n[Bot: v = %s [km/s], z = %s [m]]'
-                      % (index + 1, vtop[index], ztop[index],
-                         vbot[index], zbot[index]))
-
-            if not ztop[0] == 0:
-                print('ERROR: there must be a velocity set for z = 0 in the file %s' % filename)
-                print('e.g.:\n0    0.33\n-5    1.0\netc.')
-
-            infile.close()
-            return ztop, zbot, vtop, vbot
-
         R = 6371.
         vmin = 0.34
         decm = 0.3  # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)

pylot/core/analysis/magnitude.py

@@ -6,214 +6,261 @@ Created autumn/winter 2015.
 :author: Ludger Küperkoch / MAGS2 EP3 working group
 """
 import os
+
 import matplotlib.pyplot as plt
 import numpy as np
-from obspy.core import Stream, UTCDateTime
+import obspy.core.event as ope
-from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all
+from obspy.geodetics import degrees2kilometers
-from pylot.core.util.utils import getPatternLine
-from scipy.optimize import curve_fit
 from scipy import integrate, signal
-from pylot.core.io.data import Data
+from scipy.optimize import curve_fit
-from pylot.core.util.dataprocessing import restitute_data, read_metadata
+
+from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, \
+    select_for_phase
 from pylot.core.util.utils import common_range, fit_curve
 
+
 def richter_magnitude_scaling(delta):
     relation = np.loadtxt(os.path.join(os.path.expanduser('~'),
-                            '.pylot', 'richter_scaling.data'))
+                                       '.pylot', 'richter_scaling.data'))
     # prepare spline interpolation to calculate return value
-    func, params = fit_curve(relation[:,0], relation[:, 1])
+    func, params = fit_curve(relation[:, 0], relation[:, 1])
     return func(delta, params)
 
+
 class Magnitude(object):
-    '''
+    """
-    Superclass for calculating Wood-Anderson peak-to-peak
+    Base class object for Magnitude calculation within PyLoT.
-    amplitudes, local magnitudes, source spectra, seismic moments
+    """
-    and moment magnitudes.
-    '''
 
-    def __init__(self, wfstream, To, pwin, iplot, NLLocfile=None, \
+    def __init__(self, stream, event, verbosity=False, iplot=0):
-                 picks=None, rho=None, vp=None, Qp=None, invdir=None):
+        self._type = "M"
-        '''
+        self._plot_flag = iplot
-        :param: wfstream
+        self._verbosity = verbosity
-        :type: `~obspy.core.stream.Stream
+        self._event = event
+        self._stream = stream
+        self._magnitudes = dict()
 
-        :param: To, onset time, P- or S phase
+    def __str__(self):
-        :type:  float
+        print(
+        'number of stations used: {0}\n'.format(len(self.magnitudes.values())))
+        print('\tstation\tmagnitude')
+        for s, m in self.magnitudes.items(): print('\t{0}\t{1}'.format(s, m))
 
-        :param: pwin, pick window [To To+pwin] to get maximum
+    def __nonzero__(self):
-                peak-to-peak amplitude (WApp) or to calculate
+        return bool(self.magnitudes)
-                source spectrum (DCfc) around P onset
-        :type:  float
 
-        :param: iplot, no. of figure window for plotting interims results
+    @property
-        :type:  integer
+    def type(self):
+        return self._type
 
-        :param: NLLocfile, name and full path to NLLoc-location file
+    @property
-                needed when calling class MoMw
+    def plot_flag(self):
-        :type:  string
+        return self._plot_flag
 
-        :param: picks, dictionary containing picking results
+    @plot_flag.setter
-        :type:  dictionary
+    def plot_flag(self, value):
+        self._plot_flag = value
 
-        :param: rho [kg/m³], rock density, parameter from autoPyLoT.in
+    @property
-        :type:  integer
+    def verbose(self):
+        return self._verbosity
 
-        :param: vp [m/s], P-velocity
+    @verbose.setter
-        :param: integer
+    def verbose(self, value):
+        if not isinstance(value, bool):
+            print('WARNING: only boolean values accepted...\n')
+            value = bool(value)
+        self._verbosity = value
 
-        :param: invdir, name and path to inventory or dataless-SEED file
+    @property
-        :type:  string
+    def stream(self):
-        '''
+        return self._stream
 
-        assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream)
+    @stream.setter
+    def stream(self, value):
+        self._stream = value
 
-        self.setwfstream(wfstream)
+    @property
-        self.setTo(To)
+    def event(self):
-        self.setpwin(pwin)
+        return self._event
-        self.setiplot(iplot)
-        self.setNLLocfile(NLLocfile)
-        self.setrho(rho)
-        self.setpicks(picks)
-        self.setvp(vp)
-        self.setQp(Qp)
-        self.setinvdir(invdir)
-        self.calcwapp()
-        self.calcsourcespec()
-        self.run_calcMoMw()
 
-    def getwfstream(self):
+    @property
-        return self.wfstream
+    def origin_id(self):
+        return self._event.origins[0].resource_id
 
-    def setwfstream(self, wfstream):
+    @property
-        self.wfstream = wfstream
+    def arrivals(self):
+        return self._event.origins[0].arrivals
 
-    def getTo(self):
+    @property
-        return self.To
+    def magnitudes(self):
+        return self._magnitudes
 
-    def setTo(self, To):
+    @magnitudes.setter
-        self.To = To
+    def magnitudes(self, value):
+        """
+        takes a tuple and saves the key value pair to private
+        attribute _magnitudes
+        :param value: station, magnitude value pair
+        :type value: tuple or list
+        :return:
+        """
+        station, magnitude = value
+        self._magnitudes[station] = magnitude
 
-    def getpwin(self):
+    def calc(self):
-        return self.pwin
+        pass
 
-    def setpwin(self, pwin):
+    def updated_event(self):
-        self.pwin = pwin
+        self.event.magnitudes.append(self.net_magnitude())
+        return self.event
 
-    def getiplot(self):
+    def net_magnitude(self):
-        return self.iplot
+        if self:
-
+            # TODO if an average Magnitude instead of the median is calculated
-    def setiplot(self, iplot):
+            # StationMagnitudeContributions should be added to the returned
-        self.iplot = iplot
+            # Magnitude object
-
+            # mag_error => weights (magnitude error estimate from peak_to_peak, calcsourcespec?)
-    def setNLLocfile(self, NLLocfile):
+            # weights => StationMagnitdeContribution
-        self.NLLocfile = NLLocfile
+            mag = ope.Magnitude(
-
+                mag=np.median([M.mag for M in self.magnitudes.values()]),
-    def getNLLocfile(self):
+                magnitude_type=self.type,
-        return self.NLLocfile
+                origin_id=self.origin_id,
-
+                station_count=len(self.magnitudes),
-    def setrho(self, rho):
+                azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap)
-        self.rho = rho
+            return mag
-
+        return None
-    def getrho(self):
-        return self.rho
-
-    def setvp(self, vp):
-        self.vp = vp
-
-    def getvp(self):
-        return self.vp
-
-    def setQp(self, Qp):
-        self.Qp = Qp
-
-    def getQp(self):
-        return self.Qp
-
-    def setpicks(self, picks):
-        self.picks = picks
-
-    def getpicks(self):
-        return self.picks
-
-    def getwapp(self):
-        return self.wapp
-
-    def getw0(self):
-        return self.w0
-
-    def getfc(self):
-        return self.fc
-
-    def setinvdir(self, invdir):
-        self.invdir = invdir
-
-    def get_metadata(self):
-        return read_metadata(self.invdir)
-
-    def getpicdic(self):
-        return self.picdic
-
-    def calcwapp(self):
-        self.wapp = None
-
-    def calcsourcespec(self):
-        self.sourcespek = None
-
-    def run_calcMoMw(self):
-        self.pickdic = None
 
 
-class WApp(Magnitude):
+class RichterMagnitude(Magnitude):
-    '''
+    """
     Method to derive peak-to-peak amplitude as seen on a Wood-Anderson-
     seismograph. Has to be derived from instrument corrected traces!
-    '''
+    """
 
-    def calcwapp(self):
+    # poles, zeros and sensitivity of WA seismograph
-        print ("Getting Wood-Anderson peak-to-peak amplitude ...")
+    # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716)
-        print ("Simulating Wood-Anderson seismograph ...")
+    _paz = {
+        'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j],
+        'zeros': [0j, 0j],
+        'gain': 2080,
+        'sensitivity': 1
+    }
 
-        self.wapp = None
+    _amplitudes = dict()
-        stream = self.getwfstream()
 
-        # poles, zeros and sensitivity of WA seismograph
+    def __init__(self, stream, event, calc_win, verbosity=False, iplot=0):
-        # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716)
+        super(RichterMagnitude, self).__init__(stream, event, verbosity, iplot)
-        paz_wa = {
-            'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j],
-            'zeros': [0j, 0j],
-            'gain': 2080,
-            'sensitivity': 1}
 
-        stream.simulate(paz_remove=None, paz_simulate=paz_wa)
+        self._calc_win = calc_win
+        self._type = 'ML'
+        self.calc()
+
+    @property
+    def calc_win(self):
+        return self._calc_win
+
+    @calc_win.setter
+    def calc_win(self, value):
+        self._calc_win = value
+
+    @property
+    def amplitudes(self):
+        return self._amplitudes
+
+    @amplitudes.setter
+    def amplitudes(self, value):
+        station, a0 = value
+        self._amplitudes[station] = a0
+
+    def peak_to_peak(self, st, t0):
+
+        # simulate Wood-Anderson response
+        st.simulate(paz_remove=None, paz_simulate=self._paz)
+
+        # trim waveform to common range
+        stime, etime = common_range(st)
+        st.trim(stime, etime)
+
+        # get time delta from waveform data
+        dt = st[0].stats.delta
+
+        power = [np.power(tr.data, 2) for tr in st if tr.stats.channel[-1] not
+                 in 'Z3']
+        if len(power) != 2:
+            raise ValueError('Wood-Anderson amplitude defintion only valid for '
+                             'two horizontals: {0} given'.format(len(power)))
+        power_sum = power[0] + power[1]
+        #
+        sqH = np.sqrt(power_sum)
 
-        trH1 = stream[0].data
-        trH2 = stream[1].data
-        ilen = min([len(trH1), len(trH2)])
-        # get RMS of both horizontal components
-        sqH = np.sqrt(np.power(trH1[0:ilen], 2) + np.power(trH2[0:ilen], 2))
         # get time array
-        th = np.arange(0, len(sqH) * stream[0].stats.delta, stream[0].stats.delta)
+        th = np.arange(0, len(sqH) * dt, dt)
         # get maximum peak within pick window
-        iwin = getsignalwin(th, self.getTo(), self.getpwin())
+        iwin = getsignalwin(th, t0 - stime, self.calc_win)
-        self.wapp = np.max(sqH[iwin])
+        wapp = np.max(sqH[iwin])
-        print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp
+        if self.verbose:
+            print("Determined Wood-Anderson peak-to-peak amplitude: {0} "
+                  "mm".format(wapp))
 
-        if self.getiplot() > 1:
+        # check for plot flag (for debugging only)
-            stream.plot()
+        if self.plot_flag > 1:
+            st.plot()
             f = plt.figure(2)
             plt.plot(th, sqH)
             plt.plot(th[iwin], sqH[iwin], 'g')
-            plt.plot([self.getTo(), self.getTo()], [0, max(sqH)], 'r', linewidth=2)
+            plt.plot([t0, t0], [0, max(sqH)], 'r', linewidth=2)
-            plt.title('Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \
+            plt.title(
-                      % (stream[0].stats.station, self.wapp))
+                'Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \
+                % (st[0].stats.station, wapp))
             plt.xlabel('Time [s]')
             plt.ylabel('Displacement [mm]')
             plt.show()
             raw_input()
             plt.close(f)
 
+        return wapp
 
-class M0Mw(Magnitude):
+    def calc(self):
+        for a in self.arrivals:
+            if a.phase not in 'sS':
+                continue
+            pick = a.pick_id.get_referred_object()
+            station = pick.waveform_id.station_code
+            wf = select_for_phase(self.stream.select(
+                station=station), a.phase)
+            if not wf:
+                if self.verbose:
+                    print(
+                    'WARNING: no waveform data found for station {0}'.format(
+                        station))
+                continue
+            delta = degrees2kilometers(a.distance)
+            onset = pick.time
+            a0 = self.peak_to_peak(wf, onset)
+            amplitude = ope.Amplitude(generic_amplitude=a0 * 1e-3)
+            amplitude.unit = 'm'
+            amplitude.category = 'point'
+            amplitude.waveform_id = pick.waveform_id
+            amplitude.magnitude_hint = self.type
+            amplitude.pick_id = pick.resource_id
+            amplitude.type = 'AML'
+            self.event.amplitudes.append(amplitude)
+            self.amplitudes = (station, amplitude)
+            # using standard Gutenberg-Richter relation
+            # TODO make the ML calculation more flexible by allowing
+            # use of custom relation functions
+            magnitude = ope.StationMagnitude(
+                mag=np.log10(a0) + richter_magnitude_scaling(delta))
+            magnitude.origin_id = self.origin_id
+            magnitude.waveform_id = pick.waveform_id
+            magnitude.amplitude_id = amplitude.resource_id
+            magnitude.station_magnitude_type = self.type
+            self.event.station_magnitudes.append(magnitude)
+            self.magnitudes = (station, magnitude)
+
+
+class MomentMagnitude(Magnitude):
     '''
     Method to calculate seismic moment Mo and moment magnitude Mw.
     Requires results of class calcsourcespec for calculating plateau w0
@@ -223,55 +270,82 @@ class M0Mw(Magnitude):
     corresponding moment magntiude Mw.
     '''
 
-    def run_calcMoMw(self):
+    _props = dict()
 
-        picks = self.getpicks()
+    def __init__(self, stream, event, vp, Qp, density, verbosity=False,
-        nllocfile = self.getNLLocfile()
+                 iplot=False):
-        wfdat = self.getwfstream()
+        super(MomentMagnitude, self).__init__(stream, event, verbosity, iplot)
-        self.picdic = None
 
-        for key in picks:
+        self._vp = vp
-            if picks[key]['P']['weight'] < 4:
+        self._Qp = Qp
-                # select waveform
+        self._density = density
-                selwf = wfdat.select(station=key)
+        self._type = 'Mw'
-                if len(key) > 4:
+        self.calc()
-                    Ppattern = '%s  ?    ?    ? P' % key
-                elif len(key) == 4:
-                    Ppattern = '%s   ?    ?    ? P' % key
-                elif len(key) < 4:
-                    Ppattern = '%s    ?    ?    ? P' % key
-                nllocline = getPatternLine(nllocfile, Ppattern)
-                # get hypocentral distance, station azimuth and
-                # angle of incidence from NLLoc-location file
-                delta = float(nllocline.split(None)[21])
-                az = float(nllocline.split(None)[22])
-                inc = float(nllocline.split(None)[24])
-                # call subfunction to estimate source spectrum
-                # and to derive w0 and fc
-                [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \
-                                          self.get_metadata(), self.getvp(), delta, az, \
-                                          inc, self.getQp(), self.getiplot())
 
-                if w0 is not None:
+    @property
-                    # call subfunction to calculate Mo and Mw
+    def p_velocity(self):
-                    zdat = selwf.select(component="Z")
+        return self._vp
-                    if len(zdat) == 0:  # check for other components
-                        zdat = selwf.select(component="3")
-                    [Mo, Mw] = calcMoMw(zdat, w0, self.getrho(), self.getvp(),
-                                        delta)
-                else:
-                    Mo = None
-                    Mw = None
 
-                # add w0, fc, Mo and Mw to dictionary
+    @property
-                picks[key]['P']['w0'] = w0
+    def p_attenuation(self):
-                picks[key]['P']['fc'] = fc
+        return self._Qp
-                picks[key]['P']['Mo'] = Mo
+
-                picks[key]['P']['Mw'] = Mw
+    @property
-                self.picdic = picks
+    def rock_density(self):
+        return self._density
+
+    @property
+    def moment_props(self):
+        return self._props
+
+    @moment_props.setter
+    def moment_props(self, value):
+        station, props = value
+        self._props[station] = props
+
+    @property
+    def seismic_moment(self):
+        return self._m0
+
+    @seismic_moment.setter
+    def seismic_moment(self, value):
+        self._m0 = value
+
+    def calc(self):
+        for a in self.arrivals:
+            if a.phase not in 'pP':
+                continue
+            pick = a.pick_id.get_referred_object()
+            station = pick.waveform_id.station_code
+            wf = select_for_phase(self.stream.select(
+                station=station), a.phase)
+            if not wf:
+                continue
+            onset = pick.time
+            distance = degrees2kilometers(a.distance)
+            azimuth = a.azimuth
+            incidence = a.takeoff_angle
+            w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance,
+                                    azimuth,
+                                    incidence, self.p_attenuation,
+                                    self.plot_flag, self.verbose)
+            if w0 is None or fc is None:
+                if self.verbose:
+                    print("WARNING: insufficient frequency information")
+                continue
+            wf = select_for_phase(wf, "P")
+            m0, mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity,
+                              distance, self.verbose)
+            self.moment_props = (station, dict(w0=w0, fc=fc, Mo=m0))
+            magnitude = ope.StationMagnitude(mag=mw)
+            magnitude.origin_id = self.origin_id
+            magnitude.waveform_id = pick.waveform_id
+            magnitude.station_magnitude_type = self.type
+            self.event.station_magnitudes.append(magnitude)
+            self.magnitudes = (station, magnitude)
 
 
-def calcMoMw(wfstream, w0, rho, vp, delta):
+def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
     '''
     Subfunction of run_calcMoMw to calculate individual
     seismic moments and corresponding moment magnitudes.
@@ -295,11 +369,14 @@ def calcMoMw(wfstream, w0, rho, vp, delta):
     tr = wfstream[0]
     delta = delta * 1000  # hypocentral distance in [m]
 
-    print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \
+    if verbosity:
-          % tr.stats.station)
+        print(
+        "calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station {0} ...".format(
+            tr.stats.station))
 
     # additional common parameters for calculating Mo
-    rP = 2 / np.sqrt(15)  # average radiation pattern of P waves (Aki & Richards, 1980)
+    rP = 2 / np.sqrt(
+        15)  # average radiation pattern of P waves (Aki & Richards, 1980)
     freesurf = 2.0  # free surface correction, assuming vertical incidence
 
     Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf)
@@ -307,13 +384,16 @@ def calcMoMw(wfstream, w0, rho, vp, delta):
     # Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 # after Hanks & Kanamori (1979), defined for [dyn*cm]!
     Mw = np.log10(Mo) * 2 / 3 - 6.7  # for metric units
 
-    print("calcMoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (Mo, Mw))
+    if verbosity:
+        print(
+        "calcMoMw: Calculated seismic moment Mo = {0} Nm => Mw = {1:3.1f} ".format(
+            Mo, Mw))
 
     return Mo, Mw
 
 
-def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence,
+def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
-                   qp, iplot):
+                   qp, iplot=0, verbosity=False):
     '''
     Subfunction to calculate the source spectrum and to derive from that the plateau
     (usually called omega0) and the corner frequency assuming Aki's omega-square
@@ -321,16 +401,12 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence,
     thus restitution and integration necessary! Integrated traces are rotated
     into ray-coordinate system ZNE => LQT using Obspy's rotate modul!
 
-    :param: wfstream
+    :param: wfstream (corrected for instrument)
     :type:  `~obspy.core.stream.Stream`
 
     :param: onset, P-phase onset time
     :type: float
 
-    :param: metadata, tuple or list containing type of inventory and either
-    list of files or inventory object
-    :type:  tuple or list
-
     :param: vp, Vp-wave velocity
     :type:  float
 
@@ -349,169 +425,162 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence,
     :param: iplot, show results (iplot>1) or not (iplot<1)
     :type:  integer
     '''
-    print ("Calculating source spectrum ....")
+    if verbosity:
+        print ("Calculating source spectrum ....")
 
     # get Q value
     Q, A = qp
 
-    delta = delta * 1000  # hypocentral distance in [m]
+    dist = delta * 1000  # hypocentral distance in [m]
 
     fc = None
     w0 = None
-    data = Data()
-    wf_copy = wfstream.copy()
 
-    invtype, inventory = metadata
+    zdat = select_for_phase(wfstream, "P")
 
-    [cordat, restflag] = restitute_data(wf_copy, invtype, inventory)
+    dt = zdat[0].stats.delta
-    if restflag is True:
-        zdat = cordat.select(component="Z")
-        if len(zdat) == 0:
-            zdat = cordat.select(component="3")
-        cordat_copy = cordat.copy()
-        # get equal time stamps and lengths of traces
-        # necessary for rotation of traces
-        trstart, trend = common_range(cordat_copy)
-        cordat_copy.trim(trstart, trend)
-        try:
-            # rotate into LQT (ray-coordindate-) system using Obspy's rotate
-            # L: P-wave direction
-            # Q: SV-wave direction
-            # T: SH-wave direction
-            LQT = cordat_copy.rotate('ZNE->LQT', azimuth, incidence)
-            ldat = LQT.select(component="L")
-            if len(ldat) == 0:
-                # if horizontal channels are 2 and 3
-                # no azimuth information is available and thus no
-                # rotation is possible!
-                print("calcsourcespec: Azimuth information is missing, "
-                      "no rotation of components possible!")
-                ldat = LQT.select(component="Z")
 
-            # integrate to displacement
+    freq = zdat[0].stats.sampling_rate
-            # unrotated vertical component (for copmarison)
-            inttrz = signal.detrend(integrate.cumtrapz(zdat[0].data, None,
-                                                       zdat[0].stats.delta))
-            # rotated component Z => L
-            Ldat = signal.detrend(integrate.cumtrapz(ldat[0].data, None,
-                                                     ldat[0].stats.delta))
 
-            # get window after P pulse for
+    # trim traces to common range (for rotation)
-            # calculating source spectrum 
+    trstart, trend = common_range(wfstream)
-            tstart = UTCDateTime(zdat[0].stats.starttime)
+    wfstream.trim(trstart, trend)
-            tonset = onset.timestamp - tstart.timestamp
-            impickP = tonset * zdat[0].stats.sampling_rate
-            wfzc = Ldat[impickP: len(Ldat) - 1]
-            # get time array
-            t = np.arange(0, len(inttrz) * zdat[0].stats.delta, \
-                          zdat[0].stats.delta)
-            # calculate spectrum using only first cycles of
-            # waveform after P onset!
-            zc = crossings_nonzero_all(wfzc)
-            if np.size(zc) == 0 or len(zc) <= 3:
-                print ("calcsourcespec: Something is wrong with the waveform, "
-                       "no zero crossings derived!")
-                print ("No calculation of source spectrum possible!")
-                plotflag = 0
-            else:
-                plotflag = 1
-                index = min([3, len(zc) - 1])
-                calcwin = (zc[index] - zc[0]) * zdat[0].stats.delta
-                iwin = getsignalwin(t, tonset, calcwin)
-                xdat = Ldat[iwin]
 
-                # fft
+    # rotate into LQT (ray-coordindate-) system using Obspy's rotate
-                fny = zdat[0].stats.sampling_rate / 2
+    # L: P-wave direction
-                l = len(xdat) / zdat[0].stats.sampling_rate
+    # Q: SV-wave direction
-                # number of fft bins after Bath
+    # T: SH-wave direction
-                n = zdat[0].stats.sampling_rate * l
+    LQT = wfstream.rotate('ZNE->LQT', azimuth, incidence)
-                # find next power of 2 of data length
+    ldat = LQT.select(component="L")
-                m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
+    if len(ldat) == 0:
-                N = int(np.power(m, 2))
+        # if horizontal channels are 2 and 3
-                y = zdat[0].stats.delta * np.fft.fft(xdat, N)
+        # no azimuth information is available and thus no
-                Y = abs(y[: N / 2])
+        # rotation is possible!
-                L = (N - 1) / zdat[0].stats.sampling_rate
+        if verbosity:
-                f = np.arange(0, fny, 1 / L)
+            print("calcsourcespec: Azimuth information is missing, "
+                  "no rotation of components possible!")
+        ldat = LQT.select(component="Z")
 
-                # remove zero-frequency and frequencies above
+    # integrate to displacement
-                # corner frequency of seismometer (assumed
+    # unrotated vertical component (for comparison)
-                # to be 100 Hz)
+    inttrz = signal.detrend(integrate.cumtrapz(zdat[0].data, None, dt))
-                fi = np.where((f >= 1) & (f < 100))
-                F = f[fi]
-                YY = Y[fi]
 
-                # correction for attenuation
+    # rotated component Z => L
-                wa = 2 * np.pi * F  # angular frequency
+    Ldat = signal.detrend(integrate.cumtrapz(ldat[0].data, None, dt))
-                D = np.exp((wa * delta) / (2 * vp * Q * F ** A))
-                YYcor = YY.real * D
 
-                # get plateau (DC value) and corner frequency
+    # get window after P pulse for
-                # initial guess of plateau
+    # calculating source spectrum
-                w0in = np.mean(YYcor[0:100])
+    rel_onset = onset - trstart
-                # initial guess of corner frequency
+    impickP = int(rel_onset * freq)
-                # where spectral level reached 50% of flat level
+    wfzc = Ldat[impickP: len(Ldat) - 1]
-                iin = np.where(YYcor >= 0.5 * w0in)
+    # get time array
-                Fcin = F[iin[0][np.size(iin) - 1]]
+    t = np.arange(0, len(inttrz) * dt, dt)
+    # calculate spectrum using only first cycles of
+    # waveform after P onset!
+    zc = crossings_nonzero_all(wfzc)
+    if np.size(zc) == 0 or len(zc) <= 3:
+        if verbosity:
+            print ("calcsourcespec: Something is wrong with the waveform, "
+                   "no zero crossings derived!\n")
+            print ("No calculation of source spectrum possible!")
+        plotflag = 0
+    else:
+        plotflag = 1
+        index = min([3, len(zc) - 1])
+        calcwin = (zc[index] - zc[0]) * dt
+        iwin = getsignalwin(t, rel_onset, calcwin)
+        xdat = Ldat[iwin]
 
-                # use of implicit scipy otimization function
+        # fft
-                fit = synthsourcespec(F, w0in, Fcin)
+        fny = freq / 2
-                [optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in,
+        l = len(xdat) / freq
-                                                                        Fcin])
+        # number of fft bins after Bath
-                w01 = optspecfit[0]
+        n = freq * l
-                fc1 = optspecfit[1]
+        # find next power of 2 of data length
-                print ("calcsourcespec: Determined w0-value: %e m/Hz, \n"
+        m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
-                       "Determined corner frequency: %f Hz" % (w01, fc1))
+        N = int(np.power(m, 2))
+        y = dt * np.fft.fft(xdat, N)
+        Y = abs(y[: N / 2])
+        L = (N - 1) / freq
+        f = np.arange(0, fny, 1 / L)
 
-                # use of conventional fitting
+        # remove zero-frequency and frequencies above
-                [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot)
+        # corner frequency of seismometer (assumed
+        # to be 100 Hz)
+        fi = np.where((f >= 1) & (f < 100))
+        F = f[fi]
+        YY = Y[fi]
 
-                # get w0 and fc as median of both
+        # correction for attenuation
-                # source spectrum fits
+        wa = 2 * np.pi * F  # angular frequency
-                w0 = np.median([w01, w02])
+        D = np.exp((wa * dist) / (2 * vp * Q * F ** A))
-                fc = np.median([fc1, fc2])
+        YYcor = YY.real * D
-                print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc))
 
-        except TypeError as er:
+        # get plateau (DC value) and corner frequency
-            raise TypeError('''{0}'''.format(er))
+        # initial guess of plateau
+        w0in = np.mean(YYcor[0:100])
+        # initial guess of corner frequency
+        # where spectral level reached 50% of flat level
+        iin = np.where(YYcor >= 0.5 * w0in)
+        Fcin = F[iin[0][np.size(iin) - 1]]
 
-        if iplot > 1:
+        # use of implicit scipy otimization function
-            f1 = plt.figure()
+        fit = synthsourcespec(F, w0in, Fcin)
-            tLdat = np.arange(0, len(Ldat) * zdat[0].stats.delta, \
+        [optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, Fcin])
-                              zdat[0].stats.delta)
+        w01 = optspecfit[0]
-            plt.subplot(2, 1, 1)
+        fc1 = optspecfit[1]
-            # show displacement in mm
+        if verbosity:
-            p1, = plt.plot(t, np.multiply(inttrz, 1000), 'k')
+            print ("calcsourcespec: Determined w0-value: %e m/Hz, \n"
-            p2, = plt.plot(tLdat, np.multiply(Ldat, 1000))
+                   "Determined corner frequency: %f Hz" % (w01, fc1))
-            plt.legend([p1, p2], ['Displacement', 'Rotated Displacement'])
-            if plotflag == 1:
-                plt.plot(t[iwin], np.multiply(xdat, 1000), 'g')
-                plt.title('Seismogram and P Pulse, Station %s-%s' \
-                          % (zdat[0].stats.station, zdat[0].stats.channel))
-            else:
-                plt.title('Seismogram, Station %s-%s' \
-                          % (zdat[0].stats.station, zdat[0].stats.channel))
-            plt.xlabel('Time since %s' % zdat[0].stats.starttime)
-            plt.ylabel('Displacement [mm]')
 
-            if plotflag == 1:
+        # use of conventional fitting
-                plt.subplot(2, 1, 2)
+        [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot, verbosity)
-                p1, = plt.loglog(f, Y.real, 'k')
+
-                p2, = plt.loglog(F, YY.real)
+        # get w0 and fc as median of both
-                p3, = plt.loglog(F, YYcor, 'r')
+        # source spectrum fits
-                p4, = plt.loglog(F, fit, 'g')
+        w0 = np.median([w01, w02])
-                plt.loglog([fc, fc], [w0 / 100, w0], 'g')
+        fc = np.median([fc1, fc2])
-                plt.legend([p1, p2, p3, p4], ['Raw Spectrum', \
+        if verbosity:
-                                              'Used Raw Spectrum', \
+            print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (
-                                              'Q-Corrected Spectrum', \
+            w0, fc))
-                                              'Fit to Spectrum'])
+
-                plt.title('Source Spectrum from P Pulse, w0=%e m/Hz, fc=%6.2f Hz' \
+    if iplot > 1:
-                          % (w0, fc))
+        f1 = plt.figure()
-                plt.xlabel('Frequency [Hz]')
+        tLdat = np.arange(0, len(Ldat) * dt, dt)
-                plt.ylabel('Amplitude [m/Hz]')
+        plt.subplot(2, 1, 1)
-                plt.grid()
+        # show displacement in mm
-            plt.show()
+        p1, = plt.plot(t, np.multiply(inttrz, 1000), 'k')
-            raw_input()
+        p2, = plt.plot(tLdat, np.multiply(Ldat, 1000))
-            plt.close(f1)
+        plt.legend([p1, p2], ['Displacement', 'Rotated Displacement'])
+        if plotflag == 1:
+            plt.plot(t[iwin], np.multiply(xdat, 1000), 'g')
+            plt.title('Seismogram and P Pulse, Station %s-%s' \
+                      % (zdat[0].stats.station, zdat[0].stats.channel))
+        else:
+            plt.title('Seismogram, Station %s-%s' \
+                      % (zdat[0].stats.station, zdat[0].stats.channel))
+        plt.xlabel('Time since %s' % zdat[0].stats.starttime)
+        plt.ylabel('Displacement [mm]')
+
+        if plotflag == 1:
+            plt.subplot(2, 1, 2)
+            p1, = plt.loglog(f, Y.real, 'k')
+            p2, = plt.loglog(F, YY.real)
+            p3, = plt.loglog(F, YYcor, 'r')
+            p4, = plt.loglog(F, fit, 'g')
+            plt.loglog([fc, fc], [w0 / 100, w0], 'g')
+            plt.legend([p1, p2, p3, p4], ['Raw Spectrum', \
+                                          'Used Raw Spectrum', \
+                                          'Q-Corrected Spectrum', \
+                                          'Fit to Spectrum'])
+            plt.title('Source Spectrum from P Pulse, w0=%e m/Hz, fc=%6.2f Hz' \
+                      % (w0, fc))
+            plt.xlabel('Frequency [Hz]')
+            plt.ylabel('Amplitude [m/Hz]')
+            plt.grid()
+        plt.show()
+        raw_input()
+        plt.close(f1)
 
     return w0, fc
 
@@ -537,7 +606,7 @@ def synthsourcespec(f, omega0, fcorner):
     return ssp
 
 
-def fitSourceModel(f, S, fc0, iplot):
+def fitSourceModel(f, S, fc0, iplot, verbosity=False):
     '''
     Calculates synthetic source spectrum by varying corner frequency fc.
     Returns best approximated plateau omega0 and corner frequency, i.e. with least
@@ -591,9 +660,9 @@ def fitSourceModel(f, S, fc0, iplot):
     elif len(STD) == 0:
         fc = fc0
         w0 = max(S)
-
+    if verbosity:
-    print("fitSourceModel: best fc: %fHz, best w0: %e m/Hz" \
+        print(
-          % (fc, w0))
+        "fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0))
 
     if iplot > 1:
         plt.figure(iplot)

pylot/core/io/phases.py

@@ -15,6 +15,24 @@ from pylot.core.pick.utils import select_for_phase
 from pylot.core.util.utils import getOwner, full_range, four_digits
 
 
+def add_amplitudes(event, amplitudes):
+    amplitude_list = []
+    for pick in event.picks:
+        try:
+            a0 = amplitudes[pick.waveform_id.station_code]
+            amplitude = ope.Amplitude(generic_amplitude=a0 * 1e-3)
+            amplitude.unit = 'm'
+            amplitude.category = 'point'
+            amplitude.waveform_id = pick.waveform_id
+            amplitude.magnitude_hint = 'ML'
+            amplitude.pick_id = pick.resource_id
+            amplitude.type = 'AML'
+            amplitude_list.append(amplitude)
+        except KeyError:
+            continue
+    event.amplitudes = amplitude_list
+    return event
+
 def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
     """
     readPILOTEvent - function
@@ -549,4 +567,4 @@ def merge_picks(event, picks):
             if p.waveform_id.station_code == station and p.phase_hint == phase:
                 p.time, p.time_errors = time, err
         del time, err, phase, station
-    return event
+    return event

pylot/core/pick/autopick.py

@@ -17,10 +17,8 @@ from pylot.core.pick.charfuns import CharacteristicFunction
 from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
 from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
     getSNR, fmpicker, checkPonsets, wadaticheck
-from pylot.core.util.dataprocessing import restitute_data, read_metadata
 from pylot.core.util.utils import getPatternLine
 from pylot.core.io.data import Data
-from pylot.core.analysis.magnitude import WApp
 
 
 def autopickevent(data, param):
@@ -121,8 +119,6 @@ def autopickstation(wfstream, pickparam, verbose=False):
     # parameter to check for spuriously picked S onset
     zfac = pickparam.get('zfac')
     # path to inventory-, dataless- or resp-files
-    invdir = pickparam.get('invdir')
-    invtype, inventory = read_metadata(invdir)
 
     # initialize output
     Pweight = 4  # weight for P onset
@@ -565,24 +561,6 @@ def autopickstation(wfstream, pickparam, verbose=False):
                 # re-create stream object including both horizontal components
                 hdat = edat.copy()
                 hdat += ndat
-                h_copy = hdat.copy()
-                [cordat, restflag] = restitute_data(h_copy, invtype, inventory)
-                # calculate WA-peak-to-peak amplitude
-                # using subclass WApp of superclass Magnitude
-                if restflag:
-                    if Sweight < 4:
-                        wapp = WApp(cordat, mpickS, mpickP + sstop, iplot)
-                    else:
-                        # use larger window for getting peak-to-peak amplitude
-                        # as the S pick is quite unsure
-                        wapp = WApp(cordat, mpickP, mpickP + sstop +
-                                    (0.5 * (mpickP + sstop)), iplot)
-
-                    Ao = wapp.getwapp()
-                else:
-                    print("Go on processing data without source parameter "
-                          "determination!")
-
         else:
             msg = 'Bad initial (AIC) S-pick, skipping this onset!\n' \
                   'AIC-SNR={0}, AIC-Slope={1}counts/s\n' \
@@ -602,16 +580,6 @@ def autopickstation(wfstream, pickparam, verbose=False):
             # re-create stream object including both horizontal components
             hdat = edat.copy()
             hdat += ndat
-            h_copy = hdat.copy()
-            [cordat, restflag] = restitute_data(h_copy, invtype, inventory)
-            if restflag == 1:
-                # calculate WA-peak-to-peak amplitude
-                # using subclass WApp of superclass Magnitude
-                wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP
-                                                                     + sstop)),
-                            iplot)
-                Ao = wapp.getwapp()
-
     else:
         print('autopickstation: No horizontal component data available or ' \
               'bad P onset, skipping S picking!')

pylot/core/pick/utils.py

@@ -103,7 +103,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False):
     # 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 = symmetrize_error(diffti_te, diffti_tl)
 
     if iplot > 1:
         p = plt.figure(iplot)
@@ -325,6 +325,17 @@ def crossings_nonzero_all(data):
     return ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0]
 
 
+def symmetrize_error(dte, dtl):
+    """
+    takes earliest and latest possible pick and returns the symmetrized pick
+    uncertainty value
+    :param dte: relative lower uncertainty
+    :param dtl: relative upper uncertainty
+    :return: symmetrized error
+    """
+    return (dte + 2 * dtl) / 3
+
+
 def getSNR(X, TSNR, t1, tracenum=0):
     '''
     Function to calculate SNR of certain part of seismogram relative to

pylot/core/util/dataprocessing.py

@@ -216,12 +216,8 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False):
     for tr in data:
         seed_id = tr.get_id()
         # check, whether this trace has already been corrected
-        # TODO read actual value of processing key in Trace's stats instead
-        # of just looking for thr key: if processing is setit doesn't
-        # necessarily mean that the trace has been corrected so far but only
-        # processed in some way, e.g. normalized
         if 'processing' in tr.stats.keys() \
-                and np.all(['remove' in p for p in tr.stats.processing]) \
+                and np.any(['remove' in p for p in tr.stats.processing]) \
                 and not force:
             print("Trace {0} has already been corrected!".format(seed_id))
             continue
@@ -254,7 +250,7 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False):
                 finv = invlist[0]
             inventory = read_inventory(finv, format='STATIONXML')
         else:
-            restflag.append(False)
+            data.remove(tr)
             continue
         # apply restitution to data
         try:
@@ -270,7 +266,9 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False):
             if msg0 not in e.message or msg1 not in e.message:
                 raise
             else:
-                restflag.append(False)
+                # restitution done to copies of data thus deleting traces
+                # that failed should not be a problem
+                data.remove(tr)
                 continue
         restflag.append(True)
     # check if ALL traces could be restituted, take care of large datasets

pylot/core/util/errors.py

@@ -24,3 +24,6 @@ class OverwriteError(IOError):
 
 class ParameterError(Exception):
     pass
+
+class ProcessingError(RuntimeError):
+    pass

pylot/core/util/widgets.py

@@ -33,8 +33,8 @@ from pylot.core.pick.compare import Comparison
 from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS, LOCTOOLS, \
     COMPPOSITION_MAP
 from pylot.core.util.utils import prepTimeAxis, full_range, scaleWFData, \
-    demeanTrace, isSorted, findComboBoxIndex, clims, find_horizontals
+    demeanTrace, isSorted, findComboBoxIndex, clims
-
+import icons_rc
 
 def getDataType(parent):
     type = QInputDialog().getItem(parent, "Select phases type", "Type:",
@@ -1377,7 +1377,9 @@ class LocalisationTab(PropTab):
 
     def selectDirectory(self, edit):
         selected_directory = QFileDialog.getExistingDirectory()
-        edit.setText(selected_directory)
+        # check if string is empty
+        if selected_directory:
+            edit.setText(selected_directory)
 
     def getValues(self):
         loctool = self.locToolComboBox.currentText()