Merge branch 'develop' into update_waveformwidget

QtPyLoT.py

@@ -74,7 +74,7 @@ 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, readFilterInformation, trim_station_components, check4gaps, make_pen, pick_color_plt, \
-    pick_linestyle_plt, remove_underscores, check4doubled, identifyPhaseID, excludeQualityClasses, has_spe
+    pick_linestyle_plt, remove_underscores, check4doubled, identifyPhaseID, excludeQualityClasses, has_spe, check4rotated
 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, \
@@ -1421,6 +1421,8 @@ class MainWindow(QMainWindow):
         # check for gaps and doubled channels
         check4gaps(wfdat)
         check4doubled(wfdat)
+        # check for stations with rotated components
+        wfdat = check4rotated(wfdat, self.metadata)
         # trim station components to same start value
         trim_station_components(wfdat, trim_start=True, trim_end=False)
         self._stime = full_range(self.get_data().getWFData())[0]
@@ -1925,15 +1927,6 @@ class MainWindow(QMainWindow):
         self.apw.insert_log_widget(self.listWidget)
         self.apw.refresh_tooltips()
 
-        # self.logDockWidget = QDockWidget("AutoPickLog", self)
-        # self.logDockWidget.setObjectName("LogDockWidget")
-        # self.logDockWidget.setAllowedAreas(
-        #     Qt.LeftDockWidgetArea | Qt.RightDockWidgetArea)
-        # self.logDockWidget.setWidget(self.listWidget)
-        # self.addDockWidget(Qt.LeftDockWidgetArea, self.logDockWidget)
-        # self.addListItem('Loading default values from PyLoT-input file %s'
-        #                  % self.infile)
-
         self.apw.start.connect(self.start_autopick)
         self.apw.show()
 

autoPyLoT.py

@@ -26,7 +26,8 @@ from pylot.core.util.dataprocessing import restitute_data, read_metadata
 from pylot.core.util.defaults import SEPARATOR
 from pylot.core.util.event import Event
 from pylot.core.util.structure import DATASTRUCTURE
-from pylot.core.util.utils import real_None, remove_underscores, trim_station_components, check4gaps, check4doubled
+from pylot.core.util.utils import real_None, remove_underscores, trim_station_components, check4gaps, check4doubled, \
+    check4rotated
 from pylot.core.util.version import get_git_version as _getVersionString
 
 __version__ = _getVersionString()
@@ -254,6 +255,8 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             wfdat = check4doubled(wfdat)
             wfdat = trim_station_components(wfdat, trim_start=True, trim_end=False)
             metadata = read_metadata(parameter.get('invdir'))
+            # rotate stations to ZNE
+            wfdat = check4rotated(wfdat, metadata)
             corr_dat = None
             if locflag:
                 print("Restitute data ...")

pylot/core/io/default_parameters.py

@@ -348,6 +348,11 @@ defaults = {'rootpath': {'type': str,
                              'value': 1.0,
                              'namestring': 'Wadati tolerance'},
 
+            'jackfactor': {'type': float,
+                             'tooltip': 'pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor',
+                             'value': 5.0,
+                             'namestring': 'Jackknife safety factor'},
+
             'WAscaling': {'type': (float, float, float),
                           'tooltip': 'Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm] \
                           If zeros are set, original Richter magnitude is calculated!',
@@ -386,7 +391,7 @@ defaults = {'rootpath': {'type': str,
                          'namestring': 'Use TauPy'},
 
             'taup_model': {'type': str,
-                           'tooltip': 'define TauPy model for traveltime estimation',
+                           'tooltip': 'define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6',
                            'value': 'iasp91',
                            'namestring': 'TauPy model'}
             }
@@ -481,5 +486,6 @@ settings_special_pick = {
         'minpercent',
         'zfac',
         'mdttolerance',
-        'wdttolerance']
+        'wdttolerance',
+        'jackfactor'],
 }

pylot/core/pick/autopick.py

@@ -41,6 +41,7 @@ def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None,
     # parameter input file (usually autoPyLoT.in).
     wdttolerance = param.get('wdttolerance')
     mdttolerance = param.get('mdttolerance')
+    jackfactor = param.get('jackfactor')
     apverbose = param.get('apverbose')
     for n in range(len(data)):
         station = data[n].stats.station
@@ -80,7 +81,7 @@ def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None,
 
     # quality control
     # median check and jackknife on P-onset times
-    jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, 1, fig_dict_wadatijack)
+    jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, 1, fig_dict_wadatijack)
     #return jk_checked_onsets
     # check S-P times (Wadati)
     wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, 1, fig_dict_wadatijack)
@@ -499,6 +500,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
                                                                      SNRPdB,
                                                                      FM)
                 print(msg)
+                msg = 'autopickstation: Refined P-Pick: {} s | P-Error: {} s'.format(mpickP, Perror)
+                print(msg)
                 Sflag = 1
 
         else:
@@ -759,6 +762,9 @@ def autopickstation(wfstream, pickparam, verbose=False,
                     lpickS = lpick[ipick]
                     Serror = pickerr[ipick]
 
+                    msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(mpickS, Serror)
+                    print(msg)
+
                     # get SNR
                     [SNRS, SNRSdB, Snoiselevel] = getSNR(h_copy, tsnrh, mpickS)
 

pylot/core/pick/charfuns.py

@@ -296,9 +296,12 @@ class HOScf(CharacteristicFunction):
             elif self.getOrder() == 4:
                 LTA[j] = lta / np.power(lta1, 2)
 
-        nn = np.isnan(LTA)
+        # remove NaN's with first not-NaN-value,
-        if len(nn) > 1:
+        # so autopicker doesnt pick discontinuity at start of the trace
-            LTA[nn] = 0
+        ind = np.where(~np.isnan(LTA))[0]
+        if ind.size:
+            first = ind[0]
+            LTA[:first] = LTA[first]
         self.cf = LTA
         self.xcf = x
 

pylot/core/pick/picker.py

@@ -244,7 +244,11 @@ class AICPicker(AutoPicker):
                               & (self.Tcf >= self.Pick)) # TODO: put this in a seperate function like getsignalwin
             # find maximum within slope determination window
             # 'cause slope should be calculated up to first local minimum only!
-            dataslope = self.Data[0].data[islope[0][0]:islope[0][len(islope[0]) - 1]]
+            try:
+                dataslope = self.Data[0].data[islope[0][0:-1]]
+            except IndexError:
+                print("Slope Calculation: empty array islope, check signal window")
+                return
             if len(dataslope) < 1:
                 print('No data in slope window found!')
                 return
@@ -254,7 +258,7 @@ class AICPicker(AutoPicker):
                 # calculate slope from initial onset to maximum of AIC function
                 print("AICPicker: Not enough data samples left for slope calculation!")
                 print("Calculating slope from initial onset to maximum of AIC function ...")
-                imax = np.argmax(aicsmooth[islope[0][0]:islope[0][len(islope[0])-1]])
+                imax = np.argmax(aicsmooth[islope[0][0:-1]])
                 if imax == 0:
                     print("AICPicker: Maximum for slope determination right at the beginning of the window!")
                     print("Choose longer slope determination window!")
@@ -284,10 +288,10 @@ class AICPicker(AutoPicker):
             xslope = np.arange(0, len(dataslope), 1)
             P = np.polyfit(xslope, dataslope, 1)
             datafit = np.polyval(P, xslope)
-            if datafit[0] >= datafit[len(datafit) - 1]:
+            if datafit[0] >= datafit[-1]:
                 print('AICPicker: Negative slope, bad onset skipped!')
                 return
-            self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[len(dataslope) - 1] - datafit[0])
+            self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
 
         else:
             self.SNR = None

pylot/core/pick/utils.py

@@ -135,7 +135,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None):
         ax.axvspan(t[inoise[0]], t[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
         ax.axvspan(t[isignal[0]], t[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window')
         ax.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k', label='Noise Level')
-        ax.plot(t[isignal[zc]], np.zeros(len(zc)), '*g',
+        ax.plot(t[pis[zc]], np.zeros(len(zc)), '*g',
                 markersize=14, label='Zero Crossings')
         ax.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
         ax.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2, label='mpp')
@@ -205,8 +205,10 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None):
         t = np.arange(0, Xraw[0].stats.npts / Xraw[0].stats.sampling_rate,
                       Xraw[0].stats.delta)
         # get pick window
-        ipick = np.where(
+        ipick = np.where((t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick))
-            (t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick))
+        if len(ipick[0]) <= 1:
+            print('fmpicker: Zero crossings window to short!')
+            return
         # remove mean
         xraw[ipick] = xraw[ipick] - np.mean(xraw[ipick])
         xfilt[ipick] = xfilt[ipick] - np.mean(xfilt[ipick])
@@ -801,7 +803,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot=0, fi
     return returnflag
 
 
-def checkPonsets(pickdic, dttolerance, iplot=0, fig_dict=None):
+def checkPonsets(pickdic, dttolerance, jackfactor=5, iplot=0, fig_dict=None):
     '''
     Function to check statistics of P-onset times: Control deviation from
     median (maximum adjusted deviation = dttolerance) and apply pseudo-
@@ -838,9 +840,9 @@ def checkPonsets(pickdic, dttolerance, iplot=0, fig_dict=None):
         return
     # get pseudo variances smaller than average variances
     # (times safety factor), these picks passed jackknife test
-    ij = np.where(PHI_pseudo <= 5 * xjack)
+    ij = np.where(PHI_pseudo <= jackfactor * xjack)
     # these picks did not pass jackknife test
-    badjk = np.where(PHI_pseudo > 5 * xjack)
+    badjk = np.where(PHI_pseudo > jackfactor * xjack)
     badjkstations = np.array(stations)[badjk]
     print("checkPonsets: %d pick(s) did not pass jackknife test!" % len(badjkstations))
     print(badjkstations)

pylot/core/util/utils.py

@@ -9,6 +9,9 @@ import subprocess
 
 import numpy as np
 from obspy import UTCDateTime, read
+from obspy.signal.rotate import rotate2zne
+from obspy.io.xseed.utils import SEEDParserException
+
 from pylot.core.io.inputs import PylotParameter
 
 from scipy.interpolate import splrep, splev
@@ -696,6 +699,94 @@ def get_stations(data):
     return stations
 
 
+def check4rotated(data, metadata=None):
+
+    def rotate_components(wfstream, metadata=None):
+        """rotates components if orientation code is numeric.
+        azimut and dip are fetched from metadata"""
+        try:
+            # indexing fails if metadata is None
+            metadata[0]
+        except:
+            msg = 'Warning: could not rotate traces since no metadata was given\nset Inventory file!'
+            print(msg)
+            return wfstream
+        if metadata[0] is None:
+            # sometimes metadata is (None, (None,))
+            msg = 'Warning: could not rotate traces since no metadata was given\nCheck inventory directory!'
+            print(msg)
+            return wfstream
+        else:
+            parser = metadata[1]
+
+        def get_dip_azimut(parser, trace_id):
+            """gets azimut and dip for a trace out of the metadata parser"""
+            dip = None
+            azimut = None
+            try:
+                blockettes = parser._select(trace_id)
+            except SEEDParserException as e:
+                print(e)
+                raise ValueError
+            for blockette_ in blockettes:
+                if blockette_.id != 52:
+                    continue
+                dip = blockette_.dip
+                azimut = blockette_.azimuth
+                break
+            if dip is None or azimut is None:
+                error_msg = 'Dip and azimuth not available for trace_id {}'.format(trace_id)
+                raise ValueError(error_msg)
+            return dip, azimut
+
+        trace_ids = [trace.id for trace in wfstream]
+        for trace_id in trace_ids:
+            orientation = trace_id[-1]
+            if orientation.isnumeric():
+                # misaligned channels have a number as orientation
+                azimuts = []
+                dips = []
+                for trace_id in trace_ids:
+                    try:
+                        dip, azimut = get_dip_azimut(parser, trace_id)
+                    except ValueError as e:
+                        print(e)
+                        print('Failed to rotate station {}, no azimuth or dip available in metadata'.format(trace_id))
+                        return wfstream
+                    azimuts.append(azimut)
+                    dips.append(dip)
+                # to rotate all traces must have same length
+                wfstream = trim_station_components(wfstream, trim_start=True, trim_end=True)
+                z, n, e = rotate2zne(wfstream[0], azimuts[0], dips[0],
+                                          wfstream[1], azimuts[1], dips[1],
+                                          wfstream[2], azimuts[2], dips[2])
+                print('check4rotated: rotated station {} to ZNE'.format(trace_id))
+                z_index = dips.index(min(dips)) # get z-trace index (dip is measured from 0 to -90
+                wfstream[z_index].data = z
+                wfstream[z_index].stats.channel = wfstream[z_index].stats.channel[0:-1] + 'Z'
+                del trace_ids[z_index]
+                for trace_id in trace_ids:
+                    dip, az = get_dip_azimut(parser, trace_id)
+                    trace = wfstream.select(id=trace_id)[0]
+                    if az > 315 and az <= 45 or az > 135 and az <= 225:
+                        trace.data = n
+                        trace.stats.channel = trace.stats.channel[0:-1] + 'N'
+                    elif az > 45 and az <= 135 or az > 225 and az <= 315:
+                        trace.data = e
+                        trace.stats.channel = trace.stats.channel[0:-1] + 'E'
+                break
+            else:
+                continue
+        return wfstream
+
+    stations = get_stations(data)
+
+    for station in stations:
+        wf_station = data.select(station=station)
+        wf_station = rotate_components(wf_station, metadata)
+    return data
+
+
 def scaleWFData(data, factor=None, components='all'):
     """
     produce scaled waveforms from given waveform data and a scaling factor,

pylot/core/util/widgets.py

@@ -23,7 +23,7 @@ except:
 
 from matplotlib.figure import Figure
 from pylot.core.util.utils import find_horizontals, identifyPhase, loopIdentifyPhase, trim_station_components, \
-    identifyPhaseID
+    identifyPhaseID, check4rotated
 
 try:
     from matplotlib.backends.backend_qt4agg import FigureCanvas
@@ -440,6 +440,7 @@ class WaveformWidgetPG(QtGui.QWidget):
         self.plotWidget.showGrid(x=False, y=True, alpha=0.2)
         self.plotWidget.hideAxis('bottom')
         self.plotWidget.hideAxis('left')
+        self.wfstart, self.wfend = 0, 0
         self.reinitMoveProxy()
         self._proxy = pg.SignalProxy(self.plotWidget.scene().sigMouseMoved, rateLimit=60, slot=self.mouseMoved)
 
@@ -457,8 +458,9 @@ class WaveformWidgetPG(QtGui.QWidget):
             # if x > 0:# and index < len(data1):
             wfID = self._parent.getWFID(y)
             station = self._parent.getStationName(wfID)
+            abstime = self.wfstart + x
             if self._parent.get_current_event():
-                self.label.setText("station = {}, t = {} [s]".format(station, x))
+                self.label.setText("station = {}, T = {}, t = {} [s]".format(station, abstime, x))
             self.vLine.setPos(mousePoint.x())
             self.hLine.setPos(mousePoint.y())
 
@@ -482,7 +484,7 @@ class WaveformWidgetPG(QtGui.QWidget):
                    component='*', nth_sample=1, iniPick=None, verbosity=0):
         self.title = title
         self.clearPlotDict()
-        wfstart, wfend = full_range(wfdata)
+        self.wfstart, self.wfend = full_range(wfdata)
         nmax = 0
 
         settings = QSettings()
@@ -510,7 +512,7 @@ class WaveformWidgetPG(QtGui.QWidget):
 
         try:
             self.plotWidget.getPlotItem().vb.setLimits(xMin=float(0),
-                                                       xMax=float(wfend - wfstart),
+                                                       xMax=float(self.wfend - self.wfstart),
                                                        yMin=-0.5,
                                                        yMax=len(nsc) + 0.5)
         except:
@@ -528,7 +530,7 @@ class WaveformWidgetPG(QtGui.QWidget):
             if verbosity:
                 msg = 'plotting %s channel of station %s' % (channel, station)
                 print(msg)
-            stime = trace.stats.starttime - wfstart
+            stime = trace.stats.starttime - self.wfstart
             time_ax = prepTimeAxis(stime, trace)
             if time_ax is not None:
                 if not scaleddata:
@@ -538,9 +540,9 @@ class WaveformWidgetPG(QtGui.QWidget):
                 data = [datum + n for index, datum in enumerate(trace.data) if not index % nth_sample]
                 plots.append((times, data))
                 self.setPlotDict(n, (station, channel, network))
-        self.xlabel = 'seconds since {0}'.format(wfstart)
+        self.xlabel = 'seconds since {0}'.format(self.wfstart)
         self.ylabel = ''
-        self.setXLims([0, wfend - wfstart])
+        self.setXLims([0, self.wfend - self.wfstart])
         self.setYLims([-0.5, nmax + 0.5])
         return plots
 
@@ -2291,6 +2293,8 @@ class TuneAutopicker(QWidget):
         wfdat = self.data.getWFData()  # all available streams
         # trim station components to same start value
         trim_station_components(wfdat, trim_start=True, trim_end=False)
+        # rotate misaligned stations to ZNE
+        wfdat = check4rotated(wfdat, self.parent.metadata)
         self.stationBox.clear()
         stations = []
         for trace in self.data.getWFData():