diff --git a/QtPyLoT.py b/QtPyLoT.py index 0d596133..ca47bf07 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] 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/util/utils.py b/pylot/core/util/utils.py index 8112d9f1..646a81da 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('checkrotated: 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 47d97057..98ea3af0 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 @@ -2242,6 +2242,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.metadata) self.stationBox.clear() stations = [] for trace in self.data.getWFData():