diff --git a/QtPyLoT.py b/QtPyLoT.py index 2c36ef61..745cce04 100755 --- a/QtPyLoT.py +++ b/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() diff --git a/autoPyLoT.py b/autoPyLoT.py index 88e09c8b..44b2bd4e 100755 --- a/autoPyLoT.py +++ b/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 ...") diff --git a/pylot/core/io/default_parameters.py b/pylot/core/io/default_parameters.py index a8e648b5..f6176788 100644 --- a/pylot/core/io/default_parameters.py +++ b/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'], } diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index c5ec5d78..5b8c4d96 100644 --- a/pylot/core/pick/autopick.py +++ b/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) diff --git a/pylot/core/pick/charfuns.py b/pylot/core/pick/charfuns.py index 9d91c704..a6fb5841 100644 --- a/pylot/core/pick/charfuns.py +++ b/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) - if len(nn) > 1: - LTA[nn] = 0 + # remove NaN's with first not-NaN-value, + # so autopicker doesnt pick discontinuity at start of the trace + ind = np.where(~np.isnan(LTA))[0] + if ind.size: + first = ind[0] + LTA[:first] = LTA[first] self.cf = LTA self.xcf = x diff --git a/pylot/core/pick/picker.py b/pylot/core/pick/picker.py index 2355d4f3..1715746c 100644 --- a/pylot/core/pick/picker.py +++ b/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 diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index f7cb76cf..006f9837 100644 --- a/pylot/core/pick/utils.py +++ b/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( - (t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick)) + ipick = np.where((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) diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 8112d9f1..46a66eda 100644 --- a/pylot/core/util/utils.py +++ b/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, diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index 4da7968b..6d95c241 100644 --- a/pylot/core/util/widgets.py +++ b/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():