diff --git a/PyLoT.py b/PyLoT.py index 6e881d77..8ce799bb 100755 --- a/PyLoT.py +++ b/PyLoT.py @@ -47,7 +47,6 @@ from PySide2.QtWidgets import QMainWindow, QInputDialog, QFileDialog, \ QTreeView, QComboBox, QTabWidget, QPushButton, QGridLayout, QTableWidgetItem, QTableWidget import numpy as np from obspy import UTCDateTime, Stream -from obspy.core.event import Magnitude, Origin from obspy.core.util import AttribDict from pylot.core.util.obspyDMT_interface import check_obspydmt_structure @@ -67,21 +66,19 @@ from autoPyLoT import autoPyLoT from pylot.core.pick.compare import Comparison from pylot.core.pick.utils import getQualityFromUncertainty from pylot.core.io.phases import picksdict_from_picks +from pylot.core.io.project import Project import pylot.core.loc.nll as nll from pylot.core.util.errors import DatastructureError, \ OverwriteError from pylot.core.util.connection import checkurl from pylot.core.util.dataprocessing import Metadata, restitute_data -from pylot.core.util.utils import fnConstructor, getLogin, \ +from pylot.core.util.utils import fnConstructor, get_login, \ full_range, readFilterInformation, pick_color_plt, \ pick_linestyle_plt, identifyPhaseID, excludeQualityClasses, \ transform_colors_mpl, transform_colors_mpl_str, getAutoFilteroptions, check_all_obspy, \ check_all_pylot, get_bool, get_None, get_pylot_eventfile_with_extension from pylot.core.util.gui import make_pen -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, \ - PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \ +from pylot.core.util.widgets import FilterOptionsDialog, PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \ ComparisonWidget, TuneAutopicker, PylotParaBox, AutoPickDlg, CanvasWidget, AutoPickWidget, \ CompareEventsWidget, ProgressBarWidget, AddMetadataWidget, SingleTextLineDialog, LogWidget, PickQualitiesFromXml, \ SourceSpecWindow, ChooseWaveFormWindow, SpectrogramTab, SearchFileByExtensionDialog @@ -107,17 +104,14 @@ locateTool = dict(nll=nll) class MainWindow(QMainWindow): + _metadata: Metadata __version__ = _getVersionString() closing = Signal() def __init__(self, parent=None, infile=None, reset_qsettings=False): super(MainWindow, self).__init__(parent) - # check for default pylot.in-file - if not infile: - infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in') - print('Using default input file {}'.format(infile)) - if os.path.isfile(infile) is False: + if infile and os.path.isfile(infile) is False: infile = QFileDialog().getOpenFileName(caption='Choose PyLoT-input file')[0] if not os.path.exists(infile): @@ -181,6 +175,7 @@ class MainWindow(QMainWindow): self.autodata = Data(self) self.fnames = None + self.fnames_comp = None self._stime = None # track deleted picks for logging @@ -196,7 +191,7 @@ class MainWindow(QMainWindow): if settings.value("user/FullName", None) is None: fulluser = QInputDialog.getText(self, "Enter Name:", "Full name") settings.setValue("user/FullName", fulluser) - settings.setValue("user/Login", getLogin()) + settings.setValue("user/Login", get_login()) if settings.value("agency_id", None) is None: agency = QInputDialog.getText(self, "Enter authority/institution name:", @@ -243,20 +238,6 @@ class MainWindow(QMainWindow): self.loc = False - def init_config_files(self, infile): - pylot_config_dir = os.path.join(os.path.expanduser('~'), '.pylot') - if not os.path.exists(pylot_config_dir): - os.mkdir(pylot_config_dir) - - self._inputs = PylotParameter(infile) - if not infile: - self._inputs.reset_defaults() - # check for default pylot.in-file - infile = os.path.join(pylot_config_dir, '.pylot.in') - print('Using default input file {}'.format(infile)) - self._inputs.export2File(infile) - self.infile = infile - def setupUi(self, use_logwidget=False): try: self.startTime = min( @@ -913,7 +894,7 @@ class MainWindow(QMainWindow): return self._metadata @metadata.setter - def metadata(self, value): + def metadata(self, value: Metadata): self._metadata = value def updateFilteroptions(self): @@ -1010,13 +991,13 @@ class MainWindow(QMainWindow): events=events) if not sld.exec_(): return - fext = sld.comboBox.currentText() - # fext = '.xml' + + filenames = sld.getChecked() for event in events: - filename = get_pylot_eventfile_with_extension(event, fext) - if filename: - self.load_data(filename, draw=False, event=event, overwrite=True) - refresh = True + for filename in filenames: + if os.path.isfile(filename) and event.pylot_id in filename: + self.load_data(filename, draw=False, event=event, ask_user=True, merge_strategy=sld.merge_strategy) + refresh = True if not refresh: return if self.get_current_event().pylot_picks: @@ -1024,8 +1005,8 @@ class MainWindow(QMainWindow): self.fill_eventbox() self.setDirty(True) - def load_data(self, fname=None, loc=False, draw=True, event=None, overwrite=False): - if not overwrite: + def load_data(self, fname=None, loc=False, draw=True, event=None, ask_user=False, merge_strategy='Overwrite'): + if not ask_user: if not self.okToContinue(): return if fname is None: @@ -1039,9 +1020,12 @@ class MainWindow(QMainWindow): data = Data(self, event) try: data_new = Data(self, evtdata=str(fname)) - # MP MP commented because adding several picks might cause inconsistencies - data = data_new - # data += data_new + if merge_strategy == 'Overwrite': + data = data_new + elif merge_strategy == 'Merge': + data += data_new + else: + raise NotImplementedError(f'Unknown merge strategy: {merge_strategy}') except ValueError: qmb = QMessageBox(self, icon=QMessageBox.Question, text='Warning: Missmatch in event identifiers {} and {}. Continue?'.format( @@ -1068,9 +1052,6 @@ class MainWindow(QMainWindow): self.refreshEvents() self.setDirty(True) - def add_recentfile(self, event): - self.recentfiles.insert(0, event) - def set_button_border_color(self, button, color=None): ''' Set background color of a button. @@ -1129,16 +1110,19 @@ class MainWindow(QMainWindow): else: return - def getWFFnames_from_eventbox(self, eventbox=None): + def getWFFnames_from_eventbox(self, eventbox: str = None, subpath: str = None) -> list: ''' Return waveform filenames from event in eventbox. ''' # TODO: add dataStructure class for obspyDMT here, this is just a workaround! eventpath = self.get_current_event_path(eventbox) - basepath = eventpath.split(os.path.basename(eventpath))[0] + if subpath: + eventpath = os.path.join(eventpath, subpath) + if not os.path.isdir(eventpath): + return [] if self.dataStructure: if not eventpath: - return + return [] fnames = [os.path.join(eventpath, f) for f in os.listdir(eventpath)] else: raise DatastructureError('not specified') @@ -1176,104 +1160,117 @@ class MainWindow(QMainWindow): def getLastEvent(self): return self.recentfiles[0] - def add_events(self): + def add_events(self) -> None: ''' Creates and adds events by user selection of event folders to GUI. ''' if not self.project: self.createNewProject() - ed = GetExistingDirectories(self, 'Select event directories...') - if ed.exec_(): - eventlist = [event for event in ed.selectedFiles() if not event.endswith('EVENTS-INFO')] - basepath = eventlist[0].split(os.path.basename(eventlist[0]))[0] - # small hack: if a file "eventlist.txt" is found in the basepath use only those events specified inside - eventlist_file = os.path.join(basepath, 'eventlist.txt') - if os.path.isfile(eventlist_file): - with open(eventlist_file, 'r') as infile: - eventlist_subset = [os.path.join(basepath, filename.split('\n')[0]) for filename in - infile.readlines()] - msg = 'Found file "eventlist.txt" in database path. WILL ONLY USE SELECTED EVENTS out of {} events ' \ - 'contained in this subset' - print(msg.format(len(eventlist_subset))) - eventlist = [eventname for eventname in eventlist if eventname in eventlist_subset] - if check_obspydmt_structure(basepath): - print('Recognized obspyDMT structure in selected files. Setting Datastructure to ObspyDMT') - self.dataStructure = DATASTRUCTURE['obspyDMT']() - eventlist = check_all_obspy(eventlist) - else: - print('Setting Datastructure to PILOT') - self.dataStructure = DATASTRUCTURE['PILOT']() - eventlist = check_all_pylot(eventlist) - if not eventlist: - print('No events found! Expected structure for event folders: [eEVID.DOY.YR],\n' - ' e.g. eventID=1, doy=2, yr=2016: e0001.002.16 or obspyDMT database') - return - else: - return - if not self.project: - print('No project found.') + + eventlist, basepath = self.get_event_directories() + if not eventlist: return - # get path from first event in list and split them - path = eventlist[0] + eventlist = self.set_data_structure(basepath, eventlist) + if not eventlist: + self.show_no_events_message() + return + + dirs = self.initialize_directory_structure(eventlist[0]) + if not dirs: + return + + self.update_project_settings(dirs) + + self.project.add_eventlist(eventlist) + self.init_events() + self.setDirty(True) + + def get_event_directories(self): + ed = GetExistingDirectories(self, 'Select event directories...') + if not ed.exec_(): + return None, None + + eventlist = [event for event in ed.selectedFiles() if not event.endswith('EVENTS-INFO')] + basepath = eventlist[0].split(os.path.basename(eventlist[0]))[0] + eventlist = self.filter_eventlist_with_file(eventlist, basepath) + return eventlist, basepath + + def filter_eventlist_with_file(self, eventlist, basepath): + eventlist_file = os.path.join(basepath, 'eventlist.txt') + if os.path.isfile(eventlist_file): + with open(eventlist_file, 'r') as infile: + eventlist_subset = [os.path.join(basepath, filename.strip()) for filename in infile] + print( + f'Found file "eventlist.txt" in database path. WILL ONLY USE SELECTED EVENTS out of {len(eventlist_subset)} events contained in this subset') + eventlist = [eventname for eventname in eventlist if eventname in eventlist_subset] + return eventlist + + def set_data_structure(self, basepath, eventlist): + if check_obspydmt_structure(basepath): + print('Recognized obspyDMT structure in selected files. Setting Datastructure to ObspyDMT') + self.dataStructure = DATASTRUCTURE['obspyDMT']() + eventlist = check_all_obspy(eventlist) + else: + print('Setting Datastructure to PILOT') + self.dataStructure = DATASTRUCTURE['PILOT']() + eventlist = check_all_pylot(eventlist) + return eventlist + + def show_no_events_message(self): + print('No events found! Expected structure for event folders: [eEVID.DOY.YR],\n' + ' e.g. eventID=1, doy=2, yr=2016: e0001.002.16 or obspyDMT database') + + def initialize_directory_structure(self, path): try: system_name = platform.system() if system_name in ["Linux", "Darwin"]: dirs = { 'database': path.split('/')[-2], - 'datapath': os.path.split(path)[0], # path.split('/')[-3], + 'datapath': os.path.split(path)[0], 'rootpath': '/' + os.path.join(*path.split('/')[:-3]) } elif system_name == "Windows": rootpath = path.split('/')[:-3] rootpath[0] += '/' dirs = { - # TODO: Arrange path to meet Win standards 'database': path.split('/')[-2], 'datapath': path.split('/')[-3], 'rootpath': os.path.join(*rootpath) } except Exception as e: - dirs = { - 'database': '', - 'datapath': '', - 'rootpath': '' - } - print('Warning: Could not automatically init folder structure. ({})'.format(e)) + print(f'Warning: Could not automatically init folder structure. ({e})') + return None + self.save_directory_settings(dirs['datapath']) + return dirs + + def save_directory_settings(self, datapath): settings = QSettings() - settings.setValue("data/dataRoot", dirs['datapath']) # d irs['rootpath']) + settings.setValue("data/dataRoot", datapath) settings.sync() + def update_project_settings(self, dirs): if not self.project.eventlist: - # init parameter object - self.setParameter(show=False) - # hide all parameter (show all needed parameter later) + self.set_parameter(show=False) self.paraBox.hide_parameter() - for directory in dirs.keys(): - # set parameter - box = self.paraBox.boxes[directory] - self.paraBox.setValue(box, dirs[directory]) - # show needed parameter in box - self.paraBox.show_parameter(directory) - dirs_box = self.paraBox.get_groupbox_dialog('Directories') - if not dirs_box.exec_(): + self.update_paraBox(dirs) + if not self.paraBox.get_groupbox_dialog('Directories').exec_(): return self.project.rootpath = dirs['rootpath'] self.project.datapath = dirs['datapath'] else: - if hasattr(self.project, 'datapath'): - if not self.project.datapath == dirs['datapath']: - QMessageBox.warning(self, "PyLoT Warning", - 'Datapath missmatch to current project!') - return - else: - self.project.rootpath = dirs['rootpath'] - self.project.datapath = dirs['datapath'] + if hasattr(self.project, 'datapath') and self.project.datapath != dirs['datapath']: + QMessageBox.warning(self, "PyLoT Warning", 'Datapath mismatch to current project!') + return + self.project.rootpath = dirs['rootpath'] + self.project.datapath = dirs['datapath'] - self.project.add_eventlist(eventlist) - self.init_events() - self.setDirty(True) + def update_paraBox(self, dirs): + for directory, path in dirs.items(): + box = self.paraBox.boxes[directory] + self.paraBox.setValue(box, path) + self.paraBox.show_parameter(directory) def remove_event(self, event): qmb = QMessageBox(self, icon=QMessageBox.Question, @@ -1293,11 +1290,6 @@ class MainWindow(QMainWindow): ''' qcb = QComboBox() palette = qcb.palette() - # change highlight color: - # palette.setBrush(QtGui.QPalette.Inactive, QtGui.QPalette.Highlight, - # QtGui.QBrush(QtGui.QColor(0,0,127,127))) - # palette.setBrush(QtGui.QPalette.Active, QtGui.QPalette.Highlight, - # QtGui.QBrush(QtGui.QColor(0,0,127,127))) qcb.setPalette(palette) return qcb @@ -1372,15 +1364,40 @@ class MainWindow(QMainWindow): :type: str ''' - # if pick widget is open, refresh tooltips as well + self.refresh_tooltips() + + if not eventBox: + eventBox = self.eventBox + self.setup_eventBox(eventBox) + + current_event = self.get_current_event() + eventBox.clear() + model = eventBox.model() + + plmax = self.get_max_path_length() + index = eventBox.currentIndex() + + for id, event in enumerate(self.project.eventlist): + event_data = self.collect_event_data(event, plmax) + itemlist = self.create_event_items(event_data) + self.setItemColor(itemlist, id, event, current_event) + + if event.isTestEvent() and select_events == 'ref' or self.isEmpty(event.path): + self.disable_items(itemlist) + + model.appendRow(itemlist) + self.validate_event_path(event.path, id) + + eventBox.setCurrentIndex(index) + self.refreshRefTestButtons() + + def refresh_tooltips(self): if self.apw: self.apw.refresh_tooltips() if hasattr(self, 'cmpw'): self.cmpw.refresh_tooltips() - if not eventBox: - eventBox = self.eventBox - index = eventBox.currentIndex() + def setup_eventBox(self, eventBox): tv = QtWidgets.QTableView() header = tv.horizontalHeader() header.setSectionResizeMode(QtWidgets.QHeaderView.ResizeToContents) @@ -1388,136 +1405,110 @@ class MainWindow(QMainWindow): header.hide() tv.verticalHeader().hide() tv.setSelectionBehavior(QtWidgets.QAbstractItemView.SelectRows) - - current_event = self.get_current_event() - eventBox.setView(tv) - eventBox.clear() - model = eventBox.model() - plmax = 0 - # set maximum length of path string - for event in self.project.eventlist: - pl = len(event.path) - if pl > plmax: - plmax = pl - for id, event in enumerate(self.project.eventlist): - event_path = event.path - #phaseErrors = {'P': self._inputs['timeerrorsP'], - # 'S': self._inputs['timeerrorsS']} + def get_max_path_length(self): + if self.project.eventlist: + return max(len(event.path) for event in self.project.eventlist) + else: + return None - man_au_picks = {'manual': event.pylot_picks, - 'auto': event.pylot_autopicks} - npicks = {'manual': {'P': 0, 'S': 0}, - 'auto': {'P': 0, 'S': 0}} - npicks_total = {'manual': {'P': 0, 'S': 0}, - 'auto': {'P': 0, 'S': 0}} + def collect_event_data(self, event, plmax): + data = { + 'path': event.path, + 'picks': self.get_picks(event), + 'is_ref': event.isRefEvent(), + 'is_test': event.isTestEvent(), + 'time': None, + 'lat': None, + 'lon': None, + 'depth': None, + 'localmag': None, + 'momentmag': None, + 'notes': event.notes, + 'dirty': event.dirty, + 'plmax': plmax, + } - for ma in man_au_picks.keys(): - if man_au_picks[ma]: - for picks in man_au_picks[ma].values(): - for phasename, pick in picks.items(): - if not type(pick) in [dict, AttribDict]: - continue + if event.origins: + origin = event.origins[0] + data['time'] = origin.time + 0 # add 0 to avoid exception for time = 0s + data['lat'] = origin.latitude + data['lon'] = origin.longitude + data['depth'] = origin.depth + + if event.magnitudes: + data['momentmag'] = '{:.1f}'.format(event.magnitudes[0].mag) + data['localmag'] = '{:.1f}'.format(event.magnitudes[1].mag) if len(event.magnitudes) > 1 else ' ' + + return data + + def get_picks(self, event): + man_au_picks = {'manual': event.pylot_picks, 'auto': event.pylot_autopicks} + npicks = {'manual': {'P': 0, 'S': 0}, 'auto': {'P': 0, 'S': 0}} + npicks_total = {'manual': {'P': 0, 'S': 0}, 'auto': {'P': 0, 'S': 0}} + + for ma, picks in man_au_picks.items(): + if picks: + for pick_values in picks.values(): + for phasename, pick in pick_values.items(): + if isinstance(pick, (dict, AttribDict)): phase_ID = identifyPhaseID(phasename) - if not phase_ID in npicks[ma].keys(): - continue - if pick.get('spe'): - npicks[ma][phase_ID] += 1 - npicks_total[ma][phase_ID] += 1 + if phase_ID in npicks[ma]: + if pick.get('spe'): + npicks[ma][phase_ID] += 1 + npicks_total[ma][phase_ID] += 1 - event_ref = event.isRefEvent() - event_test = event.isTestEvent() + return npicks, npicks_total - time = lat = lon = depth = localmag = None - if len(event.origins) == 1: - origin = event.origins[0] - time = origin.time + 0 # add 0 because there was an exception for time = 0s - lat = origin.latitude - lon = origin.longitude - depth = origin.depth - if len(event.magnitudes): # if magnitude information exists, i.e., event.magnitudes has at least 1 entry - moment_magnitude = event.magnitudes[0] - momentmag = '%4.1f' % moment_magnitude.mag - if len(event.magnitudes) > 1: - local_magnitude = event.magnitudes[1] - localmag = '%4.1f' % local_magnitude.mag - else: - localmag = ' ' + def create_event_items(self, data): + items = [ + self.create_standard_item(data['path'], data['plmax'], + icon=self.style().standardIcon(QStyle.SP_DirOpenIcon)), + self.create_standard_item(data['time'].strftime("%Y-%m-%d %H:%M:%S") if data['time'] else ''), + self.create_standard_item(data['lat']), + self.create_standard_item(data['lon']), + self.create_standard_item(data['depth']), + self.create_standard_item(data['localmag']), + self.create_standard_item(data['momentmag']), + self.create_pick_item(data['picks'][0]['manual']), + self.create_pick_item(data['picks'][0]['auto']), + self.create_ref_test_item(data['is_ref'], 'ref'), + self.create_ref_test_item(data['is_test'], 'test'), + self.create_standard_item(data['notes']) + ] - else: - momentmag = ' ' - localmag = ' ' + for item in items: + item.setTextAlignment(Qt.AlignCenter) + if data['dirty'] and item == items[0]: + item.setText(f"{item.text()}*") - # text = '{path:{plen}} | manual: [{p:3d}] | auto: [{a:3d}]' - # text = text.format(path=event_path, - # plen=plmax, - # p=event_npicks, - # a=event_nautopicks) + return items - event_str = '{path:{plen}}'.format(path=event_path, plen=plmax) - if event.dirty: - event_str += '*' - item_path = QStandardItem(event_str) - item_time = QStandardItem('{}'.format(time.strftime("%Y-%m-%d %H:%M:%S") if time else '')) - item_lat = QStandardItem('{}'.format(lat)) - item_lon = QStandardItem('{}'.format(lon)) - item_depth = QStandardItem('{}'.format(depth)) - item_localmag = QStandardItem('{}'.format(localmag)) - item_momentmag = QStandardItem('{}'.format(momentmag)) + def create_standard_item(self, text, max_length=None, icon=None): + item = QStandardItem(f"{text:{max_length}}" if max_length else str(text)) + if icon: + item.setIcon(icon) + return item - item_nmp = QStandardItem() - item_nap = QStandardItem() - item_nmp.setIcon(self.manupicksicon_small) - item_nap.setIcon(self.autopicksicon_small) + def create_pick_item(self, picks): + item = QStandardItem() + item.setText(f"{picks['P']}|{picks['S']}") + return item - for picktype, item_np in [('manual', item_nmp), ('auto', item_nap)]: - npicks_str = f"{npicks[picktype]['P']}|{npicks[picktype]['S']}" - #npicks_str += f"({npicks_total[picktype]['P']}/{npicks_total[picktype]['S']})" - item_np.setText(npicks_str) + def create_ref_test_item(self, condition, color_key): + item = QStandardItem() + if condition: + item.setBackground(self._ref_test_colors[color_key]) + return item - item_ref = QStandardItem() # str(event_ref)) - item_test = QStandardItem() # str(event_test)) - if event_ref: - item_ref.setBackground(self._ref_test_colors['ref']) - if event_test: - item_test.setBackground(self._ref_test_colors['test']) - item_notes = QStandardItem(event.notes) + def disable_items(self, items): + for item in items: + item.setEnabled(False) - openIcon = self.style().standardIcon(QStyle.SP_DirOpenIcon) - item_path.setIcon(openIcon) - # if ref: set different color e.g. - # if event_ref: - # item.setBackground(self._colors['ref']) - # if event_test: - # itemt.setBackground(self._colors['test']) - # item.setForeground(QtGui.QColor('black')) - # font = item.font() - # font.setPointSize(10) - # item.setFont(font) - # item2.setForeground(QtGui.QColor('black')) - # item2.setFont(font) - itemlist = [item_path, item_time, item_lat, item_lon, item_depth, - item_localmag, item_momentmag, item_nmp, item_nap, item_ref, item_test, item_notes] - for item in itemlist: - item.setTextAlignment(Qt.AlignCenter) - if event_test and select_events == 'ref' or self.isEmpty(event_path): - for item in itemlist: - item.setEnabled(False) - - # item color - self.setItemColor(itemlist, id, event, current_event) - - model.appendRow(itemlist) - if not event.path == self.eventBox.itemText(id).split('*')[0].strip(): - message = ('Path missmatch creating eventbox.\n' - '{} unequal {}.' - .format(event.path, self.eventBox.itemText(id))) - raise ValueError(message) - # not working with obspy events - # eventBox.setItemData(id, event) - eventBox.setCurrentIndex(index) - self.refreshRefTestButtons() + def validate_event_path(self, path, id): + if path != self.eventBox.itemText(id).split('*')[0].strip(): + raise ValueError(f"Path mismatch creating eventbox. {path} unequal {self.eventBox.itemText(id)}") def isEmpty(self, event_path): wf_stat = {True: 'processed', @@ -1549,11 +1540,6 @@ class MainWindow(QMainWindow): fname = str(action.data().toString()) return fname - def getEventFileName(self, type='manual'): - if self.get_fnames(type) is None: - self.set_fname(self.get_data().getEventFileName(), type) - return self.get_fnames(type) - def saveData(self, event=None, directory=None, outformats=['.xml', '.cnv', '.obs', '_focmec.in', '.pha']): ''' Save event data to directory with specified output formats. @@ -1787,7 +1773,7 @@ class MainWindow(QMainWindow): def getStime(self): if self.get_data(): - return full_range(self.get_data().getWFData())[0] + return full_range(self.get_data().get_wf_data())[0] def addActions(self, target, actions): for action in actions: @@ -1960,26 +1946,25 @@ class MainWindow(QMainWindow): def prepareLoadWaveformData(self): self.fnames = self.getWFFnames_from_eventbox() - self.fnames_syn = [] + self.fnames_comp = [] + fnames_comp = self.getWFFnames_from_eventbox(subpath='compare') + self.dataPlot.activateCompareOptions(bool(fnames_comp)) + if fnames_comp: + if self.dataPlot.comp_checkbox.isChecked(): + self.fnames_comp = fnames_comp + eventpath = self.get_current_event_path() basepath = eventpath.split(os.path.basename(eventpath))[0] self.obspy_dmt = check_obspydmt_structure(basepath) self.dataPlot.activateObspyDMToptions(self.obspy_dmt) if self.obspy_dmt: self.prepareObspyDMT_data(eventpath) + self.dataPlot.activateCompareOptions(True) def loadWaveformData(self): ''' Load waveform data corresponding to current selected event. ''' - # if self.fnames and self.okToContinue(): - # self.setDirty(True) - # ans = self.data.setWFData(self.fnames) - # elif self.fnames is None and self.okToContinue(): - # ans = self.data.setWFData(self.getWFFnames()) - # else: - # ans = False - settings = QSettings() # process application events to wait for event items to appear in event box QApplication.processEvents() @@ -1990,16 +1975,16 @@ class MainWindow(QMainWindow): if len(curr_event.origins) > 0: origin_time = curr_event.origins[0].time - tstart = settings.value('tstart') if get_None(settings.value('tstart')) else 0 - tstop = settings.value('tstop') if get_None(settings.value('tstop')) else 0 + tstart = settings.value('tstart') if get_none(settings.value('tstart')) else 0 + tstop = settings.value('tstop') if get_none(settings.value('tstop')) else 0 tstart = origin_time + float(tstart) tstop = origin_time + float(tstop) else: tstart = None tstop = None - self.data.setWFData(self.fnames, - self.fnames_syn, + self.data.set_wf_data(self.fnames, + self.fnames_comp, checkRotated=True, metadata=self.metadata, tstart=tstart, @@ -2007,7 +1992,7 @@ class MainWindow(QMainWindow): def prepareObspyDMT_data(self, eventpath): qcbox_processed = self.dataPlot.qcombo_processed - qcheckb_syn = self.dataPlot.syn_checkbox + qcheckb_syn = self.dataPlot.comp_checkbox qcbox_processed.setEnabled(False) qcheckb_syn.setEnabled(False) for fpath in os.listdir(eventpath): @@ -2015,8 +2000,8 @@ class MainWindow(QMainWindow): if 'syngine' in fpath: eventpath_syn = os.path.join(eventpath, fpath) qcheckb_syn.setEnabled(True) - if self.dataPlot.syn_checkbox.isChecked(): - self.fnames_syn = [os.path.join(eventpath_syn, filename) for filename in os.listdir(eventpath_syn)] + if self.dataPlot.comp_checkbox.isChecked(): + self.fnames_comp = [os.path.join(eventpath_syn, filename) for filename in os.listdir(eventpath_syn)] if 'processed' in fpath: qcbox_processed.setEnabled(True) if qcbox_processed.isEnabled(): @@ -2027,21 +2012,6 @@ class MainWindow(QMainWindow): eventpath_dmt = os.path.join(eventpath, wftype) self.fnames = [os.path.join(eventpath_dmt, filename) for filename in os.listdir(eventpath_dmt)] - # def setWFstatus(self): - # ''' - # show status of current data, can be either 'raw' or 'processed' - # :param status: - # :return: - # ''' - # status = self.data.processed - # wf_stat_color = {True: 'green', - # False: 'black', - # None: None} - # wf_stat = {True: 'processed', - # False: 'raw', - # None: None} - # self.dataPlot.setQCboxItem(wf_stat[status]) - def check_plot_quantity(self): """ Check the amount of samples to be plotted and ask user to reduce the amount if it is too large. @@ -2061,29 +2031,11 @@ class MainWindow(QMainWindow): self.update_status('Dataset is very large. Using fast plotting method (MIN/MAX)', 10000) settings.sync() return True - # nth_sample_new = int(np.ceil(npts/npts_max)) - # message = "You are about to plot a huge dataset with {npts} datapoints. With a current setting of " \ - # "nth_sample = {nth_sample} a total of {npts2plot} points will be plotted which is more " \ - # "than the maximum setting of {npts_max}. " \ - # "PyLoT recommends to raise nth_sample from {nth_sample} to {nth_sample_new}. Do you want "\ - # "to change nth_sample to {nth_sample_new} now?" - # - # ans = QMessageBox.question(self, self.tr("Optimize plot performance..."), - # self.tr(message.format(npts=npts, - # nth_sample=nth_sample, - # npts_max=npts_max, - # nth_sample_new=nth_sample_new, - # npts2plot=npts2plot)), - # QMessageBox.Yes | QMessageBox.No, - # QMessageBox.Yes) - # if ans == QMessageBox.Yes: - # settings.setValue("nth_sample", nth_sample_new) - # settings.sync() def get_npts_to_plot(self): if not hasattr(self.data, 'wfdata'): return 0 - return sum(trace.stats.npts for trace in self.data.getWFData()) + return sum(trace.stats.npts for trace in self.data.get_wf_data()) def connectWFplotEvents(self): ''' @@ -2296,17 +2248,17 @@ class MainWindow(QMainWindow): zne_text = {'Z': 'vertical', 'N': 'north-south', 'E': 'east-west'} comp = self.getComponent() title = 'section: {0} components'.format(zne_text[comp]) - wfst = self.get_data().getWFData() - wfsyn = self.get_data().getSynWFData() + wfst = self.get_data().get_wf_data() + wfsyn = self.get_data().getAltWFdata() if self.filterActionP.isChecked() and filter: self.filterWaveformData(plot=False, phase='P') elif self.filterActionS.isChecked() and filter: self.filterWaveformData(plot=False, phase='S') - # wfst = self.get_data().getWFData().select(component=comp) - # wfst += self.get_data().getWFData().select(component=alter_comp) + # wfst = self.get_data().get_wf_data().select(component=comp) + # wfst += self.get_data().get_wf_data().select(component=alter_comp) plotWidget = self.getPlotWidget() self.adjustPlotHeight() - if get_bool(settings.value('large_dataset')): + if get_bool(settings.value('large_dataset')) == True: self.plot_method = 'fast' else: self.plot_method = 'normal' @@ -2318,7 +2270,7 @@ class MainWindow(QMainWindow): def adjustPlotHeight(self): if self.pg: return - height_need = len(self.data.getWFData()) * self.height_factor + height_need = len(self.data.get_wf_data()) * self.height_factor plotWidget = self.getPlotWidget() if self.tabs.widget(0).frameSize().height() < height_need: plotWidget.setMinimumHeight(height_need) @@ -2328,40 +2280,34 @@ class MainWindow(QMainWindow): def plotZ(self): self.setComponent('Z') self.plotWaveformDataThread() - # self.drawPicks() - # self.draw() def plotN(self): self.setComponent('N') self.plotWaveformDataThread() - # self.drawPicks() - # self.draw() def plotE(self): self.setComponent('E') self.plotWaveformDataThread() - # self.drawPicks() - # self.draw() def pushFilterWF(self, param_args): - self.get_data().filterWFData(param_args) + self.get_data().filter_wf_data(param_args) def filterP(self): self.filterActionS.setChecked(False) if self.filterActionP.isChecked(): self.filterWaveformData(phase='P') else: - self.resetWFData() + self.reset_wf_data() def filterS(self): self.filterActionP.setChecked(False) if self.filterActionS.isChecked(): self.filterWaveformData(phase='S') else: - self.resetWFData() + self.reset_wf_data() - def resetWFData(self): - self.get_data().resetWFData() + def reset_wf_data(self): + self.get_data().reset_wf_data() self.plotWaveformDataThread() def filterWaveformData(self, plot=True, phase=None): @@ -2380,15 +2326,13 @@ class MainWindow(QMainWindow): kwargs = self.getFilterOptions()[phase].parseFilterOptions() self.pushFilterWF(kwargs) else: - self.get_data().resetWFData() + self.get_data().reset_wf_data() elif self.filterActionP.isChecked() or self.filterActionS.isChecked(): self.adjustFilterOptions() else: - self.get_data().resetWFData() + self.get_data().reset_wf_data() if plot: self.plotWaveformDataThread(filter=False) - # self.drawPicks() - # self.draw() def getAutoFilteroptions(self, phase): return getAutoFilteroptions(phase, self._inputs) @@ -2420,20 +2364,8 @@ class MainWindow(QMainWindow): def getFilterOptions(self): return self.filteroptions - # try: - # return self.filteroptions[self.getSeismicPhase()] - # except AttributeError as e: - # print(e) - # return FilterOptions(None, None, None) - def getFilters(self): - return self.filteroptions - - def setFilterOptions(self, filterOptions): # , seismicPhase=None): - # if seismicPhase is None: - # self.getFilterOptions()[self.getSeismicPhase()] = filterOptions - # else: - # self.getFilterOptions()[seismicPhase] = filterOptions + def setFilterOptions(self, filterOptions): self.filterOptions = filterOptions filterP = filterOptions['P'] filterS = filterOptions['S'] @@ -2462,28 +2394,6 @@ class MainWindow(QMainWindow): self.checkFilterOptions() - # def updateFilterOptions(self): - # try: - # settings = QSettings() - # if settings.value("filterdefaults", - # None) is None and not self.getFilters(): - # for key, value in FILTERDEFAULTS.items(): - # self.setFilterOptions(FilterOptions(**value), key) - # elif settings.value("filterdefaults", None) is not None: - # for key, value in settings.value("filterdefaults"): - # self.setFilterOptions(FilterOptions(**value), key) - # except Exception as e: - # self.update_status('Error ...') - # emsg = QErrorMessage(self) - # emsg.showMessage('Error: {0}'.format(e)) - # else: - # self.update_status('Filter loaded ... ' - # '[{0}: {1} Hz]'.format( - # self.getFilterOptions().getFilterType(), - # self.getFilterOptions().getFreq())) - # if self.filterActionP.isChecked() or self.filterActionS.isChecked(): - # self.filterWaveformData() - def getSeismicPhase(self): return self.seismicPhase @@ -2509,11 +2419,6 @@ class MainWindow(QMainWindow): def alterPhase(self): pass - def setSeismicPhase(self, phase): - self.seismicPhase = self.seismicPhaseButtonGroup.getValue() - self.update_status('Seismic phase changed to ' - '{0}'.format(self.getSeismicPhase())) - def scrollPlot(self, gui_event): ''' Function connected to mouse wheel scrolling inside WFdataPlot. @@ -2609,24 +2514,27 @@ class MainWindow(QMainWindow): print("Warning! No network, station, and location info available!") return self.update_status('picking on station {0}'.format(station)) - data = self.get_data().getOriginalWFData().copy() + wfdata = self.get_data().getOriginalWFData().copy() + wfdata_comp = self.get_data().getAltWFdata().copy() event = self.get_current_event() wftype = self.dataPlot.qcombo_processed.currentText() if self.obspy_dmt else None pickDlg = PickDlg(self, parameter=self._inputs, - data=data.select(station=station), + data=wfdata.select(station=station), + data_compare=wfdata_comp.select(station=station), station=station, network=network, location=location, picks=self.getPicksOnStation(station, 'manual'), autopicks=self.getPicksOnStation(station, 'auto'), metadata=self.metadata, event=event, model=self.inputs.get('taup_model'), - filteroptions=self.filteroptions, wftype=wftype) + filteroptions=self.filteroptions, wftype=wftype, + show_comp_data=self.dataPlot.comp_checkbox.isChecked()) if self.filterActionP.isChecked(): pickDlg.currentPhase = "P" - pickDlg.filterWFData() + pickDlg.filter_wf_data() elif self.filterActionS.isChecked(): pickDlg.currentPhase = "S" - pickDlg.filterWFData() + pickDlg.filter_wf_data() pickDlg.nextStation.setChecked(self.nextStation) pickDlg.nextStation.stateChanged.connect(self.toggle_next_station) if pickDlg.exec_(): @@ -2926,11 +2834,6 @@ class MainWindow(QMainWindow): automanu[type](station=station, pick=picks) return rval - def deletePicks(self, station, deleted_pick, type): - deleted_pick['station'] = station - self.addPicks(station, {}, type=type) - self.log_deleted_picks([deleted_pick]) - def log_deleted_picks(self, deleted_picks, event_path=None): ''' Log deleted picks to list self.deleted_picks @@ -2999,10 +2902,16 @@ class MainWindow(QMainWindow): event = self.get_current_event() event.pylot_picks = {} event.pylot_autopicks = {} - picksdict = picksdict_from_picks(evt=self.get_data().get_evt_data()) + picksdict = picksdict_from_picks(evt=self.get_data().get_evt_data(), parameter=self.getParameter()) event.addPicks(picksdict['manual']) event.addAutopicks(picksdict['auto']) + def getParameter(self): + if hasattr(self.project, 'parameter') and isinstance(self.project.parameter, PylotParameter): + return self.project.parameter + else: + return self._inputs + def drawPicks(self, station=None, picktype=None, stime=None): # if picktype not specified, draw both if not stime: @@ -3060,10 +2969,6 @@ class MainWindow(QMainWindow): elif phaseID == 'S': quality = getQualityFromUncertainty(picks['spe'], self._inputs['timeerrorsS']) - # quality = getPickQuality(self.get_data().getWFData(), - # stat_picks, self._inputs, phaseID, - # compclass) - if not picks['mpp']: continue mpp = picks['mpp'] - stime @@ -3573,7 +3478,7 @@ class MainWindow(QMainWindow): if not self.metadata: return None - wf_copy = self.get_data().getWFData().copy() + wf_copy = self.get_data().get_wf_data().copy() wf_select = Stream() # restitute only picked traces @@ -3613,19 +3518,6 @@ class MainWindow(QMainWindow): def check4Loc(self): return self.picksNum() >= 4 - # def check4Comparison(self): - # mpicks = self.getPicks() - # apicks = self.getPicks('auto') - # for station, phases in mpicks.items(): - # try: - # aphases = apicks[station] - # for phase in phases.keys(): - # if phase in aphases.keys(): - # return True - # except KeyError: - # continue - # return False - def picksNum(self, type='manual'): num = 0 for phases in self.getPicks(type).values(): @@ -3666,16 +3558,6 @@ class MainWindow(QMainWindow): def show_event_information(self): pass - def createNewEvent(self): - if self.okToContinue(): - new = NewEventDlg() - if new.exec_() != QDialog.Rejected: - evtpar = new.getValues() - cinfo = create_creation_info(agency_id=self.agency) - event = create_event(evtpar['origintime'], cinfo) - self.data = Data(self, evtdata=event) - self.setDirty(True) - def createNewProject(self): ''' Create new project file. @@ -3819,16 +3701,6 @@ class MainWindow(QMainWindow): self.dirty = value self.fill_eventbox() - def closeEvent(self, event): - if self.okToContinue(): - if hasattr(self, 'logwidget'): - self.logwidget.close() - event.accept() - else: - event.ignore() - # self.closing.emit() - # QMainWindow.closeEvent(self, event) - def setParameter(self, checked=0, show=True): if checked: pass # dummy argument to receive trigger signal (checked) if called by QAction if not self.paraBox: @@ -3878,181 +3750,6 @@ class MainWindow(QMainWindow): form = HelpForm(self, ':/help.html') form.show() - -class Project(object): - ''' - Pickable class containing information of a PyLoT project, like event lists and file locations. - ''' - - # TODO: remove rootpath - def __init__(self): - self.eventlist = [] - self.location = None - self.rootpath = None - self.datapath = None - self.dirty = False - self.parameter = None - self._table = None - - def add_eventlist(self, eventlist): - ''' - Add events from an eventlist containing paths to event directories. - Will skip existing paths. - ''' - if len(eventlist) == 0: - return - for item in eventlist: - event = Event(item) - event.rootpath = self.parameter['rootpath'] - event.database = self.parameter['database'] - event.datapath = self.parameter['datapath'] - if not event.path in self.getPaths(): - self.eventlist.append(event) - self.setDirty() - else: - print('Skipping event with path {}. Already part of project.'.format(event.path)) - self.eventlist.sort(key=lambda x: x.pylot_id) - self.search_eventfile_info() - - def remove_event(self, event): - self.eventlist.remove(event) - - def remove_event_by_id(self, eventID): - for event in self.eventlist: - if eventID in str(event.resource_id): - self.remove_event(event) - break - - def read_eventfile_info(self, filename, separator=','): - ''' - Try to read event information from file (:param:filename) comparing specific event datetimes. - File structure (each row): event, date, time, magnitude, latitude, longitude, depth - separated by :param:separator each. - ''' - with open(filename, 'r') as infile: - for line in infile.readlines(): - eventID, date, time, mag, lat, lon, depth = line.split(separator)[:7] - # skip first line - try: - day, month, year = date.split('/') - except: - continue - year = int(year) - # hardcoded, if year only consists of 2 digits (e.g. 16 instead of 2016) - if year < 100: - year += 2000 - datetime = '{}-{}-{}T{}'.format(year, month, day, time) - try: - datetime = UTCDateTime(datetime) - except Exception as e: - print(e, datetime, filename) - continue - for event in self.eventlist: - if eventID in str(event.resource_id) or eventID in event.origins: - if event.origins: - origin = event.origins[0] # should have only one origin - if origin.time == datetime: - origin.latitude = float(lat) - origin.longitude = float(lon) - origin.depth = float(depth) - else: - continue - elif not event.origins: - origin = Origin(resource_id=event.resource_id, - time=datetime, latitude=float(lat), - longitude=float(lon), depth=float(depth)) - event.origins.append(origin) - event.magnitudes.append(Magnitude(resource_id=event.resource_id, - mag=float(mag), - mag_type='M')) - break - - def search_eventfile_info(self): - ''' - Search all datapaths in rootpath for filenames with given file extension fext - and try to read event info from it - ''' - datapaths = [] - fext = '.csv' - for event in self.eventlist: - if not event.datapath in datapaths: - datapaths.append(event.datapath) - for datapath in datapaths: - # datapath = os.path.join(self.rootpath, datapath) - if os.path.isdir(datapath): - for filename in os.listdir(datapath): - filename = os.path.join(datapath, filename) - if os.path.isfile(filename) and filename.endswith(fext): - try: - self.read_eventfile_info(filename) - except Exception as e: - print('Failed on reading eventfile info from file {}: {}'.format(filename, e)) - else: - print("Directory %s does not exist!" % datapath) - - def getPaths(self): - ''' - Returns paths (eventlist) of all events saved in the project. - ''' - paths = [] - for event in self.eventlist: - paths.append(event.path) - return paths - - def setDirty(self, value=True): - self.dirty = value - - def getEventFromPath(self, path): - ''' - Search for an event in the project by event path. - ''' - for event in self.eventlist: - if event.path == path: - return event - - def save(self, filename=None): - ''' - Save PyLoT Project to a file. - Can be loaded by using project.load(filename). - ''' - try: - import pickle - except ImportError: - import _pickle as pickle - - if filename: - self.location = filename - else: - filename = self.location - - table = self._table # MP: see below - try: - outfile = open(filename, 'wb') - self._table = [] # MP: Workaround as long as table cannot be saved as part of project - pickle.dump(self, outfile, protocol=pickle.HIGHEST_PROTOCOL) - self.setDirty(False) - self._table = table # MP: see above - return True - except Exception as e: - print('Could not pickle PyLoT project. Reason: {}'.format(e)) - self.setDirty() - self._table = table # MP: see above - return False - - @staticmethod - def load(filename): - ''' - Load project from filename. - ''' - import pickle - infile = open(filename, 'rb') - project = pickle.load(infile) - infile.close() - project.location = filename - print('Loaded %s' % filename) - return project - - class GetExistingDirectories(QFileDialog): ''' File dialog with possibility to select multiple folders. @@ -4080,8 +3777,6 @@ def create_window(): app.setOrganizationDomain("rub.de") app.setApplicationName("PyLoT") app.references = set() - # app.references.add(window) - # window.show() return app, app_created @@ -4089,7 +3784,6 @@ def main(project_filename=None, pylot_infile=None, reset_qsettings=False): # create the Qt application pylot_app, app_created = create_window() - # pylot_app = QApplication(sys.argv) pixmap = QPixmap(":/splash/splash.png") splash = QSplashScreen(pixmap) splash.show() diff --git a/autoPyLoT.py b/autoPyLoT.py index 00da0797..cc86a786 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -28,7 +28,7 @@ from pylot.core.util.dataprocessing import restitute_data, 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 get_None, trim_station_components, check4gapsAndRemove, check4doubled, \ +from pylot.core.util.utils import get_none, trim_station_components, check4gapsAndRemove, check4doubled, \ check4rotated from pylot.core.util.version import get_git_version as _getVersionString @@ -91,9 +91,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even sp=sp_info) print(splash) - parameter = get_None(parameter) - inputfile = get_None(inputfile) - eventid = get_None(eventid) + parameter = get_none(parameter) + inputfile = get_none(inputfile) + eventid = get_none(eventid) fig_dict = None fig_dict_wadatijack = None @@ -119,13 +119,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even obspyDMT_wfpath = input_dict['obspyDMT_wfpath'] if not parameter: - if inputfile: - parameter = PylotParameter(inputfile) - # iplot = parameter['iplot'] - else: - infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in') - print('Using default input file {}'.format(infile)) - parameter = PylotParameter(infile) + if not inputfile: + print('Using default input parameter') + parameter = PylotParameter(inputfile) else: if not type(parameter) == PylotParameter: print('Wrong input type for parameter: {}'.format(type(parameter))) @@ -154,7 +150,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even datastructure.setExpandFields(exf) # check if default location routine NLLoc is available and all stations are used - if get_None(parameter['nllocbin']) and station == 'all': + if get_none(parameter['nllocbin']) and station == 'all': locflag = 1 # get NLLoc-root path nllocroot = parameter.get('nllocroot') @@ -247,7 +243,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even pylot_event = Event(eventpath) # event should be path to event directory data.setEvtData(pylot_event) if fnames == 'None': - data.setWFData(glob.glob(os.path.join(datapath, event_datapath, '*'))) + data.set_wf_data(glob.glob(os.path.join(datapath, event_datapath, '*'))) # the following is necessary because within # multiple event processing no event ID is provided # in autopylot.in @@ -262,7 +258,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even now.minute) parameter.setParam(eventID=eventID) else: - data.setWFData(fnames) + data.set_wf_data(fnames) eventpath = events[0] # now = datetime.datetime.now() @@ -272,7 +268,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even # now.hour, # now.minute) parameter.setParam(eventID=eventid) - wfdat = data.getWFData() # all available streams + wfdat = data.get_wf_data() # all available streams if not station == 'all': wfdat = wfdat.select(station=station) if not wfdat: diff --git a/pylot/core/io/data.py b/pylot/core/io/data.py index b085d181..26a2e72a 100644 --- a/pylot/core/io/data.py +++ b/pylot/core/io/data.py @@ -1,658 +1,69 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -import copy import os +from dataclasses import dataclass, field +from typing import List, Union -from PySide2.QtWidgets import QMessageBox -from obspy import read_events -from obspy.core import read, Stream, UTCDateTime -from obspy.core.event import Event as ObsPyEvent -from obspy.io.sac import SacIOError +from obspy import UTCDateTime -import pylot.core.loc.focmec as focmec -import pylot.core.loc.hypodd as hypodd -import pylot.core.loc.velest as velest -from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \ - picksdict_from_pilot, merge_picks, PylotParameter -from pylot.core.util.errors import FormatError, OverwriteError -from pylot.core.util.event import Event -from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT -from pylot.core.util.utils import fnConstructor, full_range, check4rotated, \ - check_for_gaps_and_merge, trim_station_components, check_for_nan +from pylot.core.io.event import EventData +from pylot.core.io.waveformdata import WaveformData +from pylot.core.util.dataprocessing import Metadata - -class Data(object): - """ - Data container with attributes wfdata holding ~obspy.core.stream. - - :type parent: PySide2.QtWidgets.QWidget object, optional - :param parent: A PySide2.QtWidgets.QWidget object utilized when - called by a GUI to display a PySide2.QtWidgets.QMessageBox instead of printing - to standard out. - :type evtdata: ~obspy.core.event.Event object, optional - :param evtdata ~obspy.core.event.Event object containing all derived or - loaded event. Container object holding, e.g. phase arrivals, etc. - """ +@dataclass +class Data: + event_data: EventData = field(default_factory=EventData) + waveform_data: WaveformData = field(default_factory=WaveformData) + metadata: Metadata = field(default_factory=Metadata) + _parent: Union[None, 'QtWidgets.QWidget'] = None def __init__(self, parent=None, evtdata=None): self._parent = parent - if self.getParent(): - self.comp = parent.getComponent() - else: - self.comp = 'Z' - self.wfdata = Stream() - self._new = False - if isinstance(evtdata, ObsPyEvent) or isinstance(evtdata, Event): - pass - elif isinstance(evtdata, dict): - evt = readPILOTEvent(**evtdata) - evtdata = evt - elif isinstance(evtdata, str): - try: - cat = read_events(evtdata) - if len(cat) != 1: - raise ValueError('ambiguous event information for file: ' - '{file}'.format(file=evtdata)) - evtdata = cat[0] - except TypeError as e: - if 'Unknown format for file' in e.message: - if 'PHASES' in evtdata: - picks = picksdict_from_pilot(evtdata) - evtdata = ObsPyEvent() - evtdata.picks = picks_from_picksdict(picks) - elif 'LOC' in evtdata: - raise NotImplementedError('PILOT location information ' - 'read support not yet ' - 'implemented.') - elif 'event.pkl' in evtdata: - evtdata = qml_from_obspyDMT(evtdata) - else: - raise e - else: - raise e - else: # create an empty Event object - self.setNew() - evtdata = ObsPyEvent() - evtdata.picks = [] - self.evtdata = evtdata - self.wforiginal = None - self.cuttimes = None - self.dirty = False - self.processed = None + self.event_data = EventData(evtdata) + self.waveform_data = WaveformData() def __str__(self): - return str(self.wfdata) + return str(self.waveform_data.wfdata) def __add__(self, other): - assert isinstance(other, Data), "operands must be of same type 'Data'" - rs_id = self.get_evt_data().get('resource_id') - rs_id_other = other.get_evt_data().get('resource_id') - if other.isNew() and not self.isNew(): - picks_to_add = other.get_evt_data().picks - old_picks = self.get_evt_data().picks - wf_ids_old = [pick.waveform_id for pick in old_picks] - for new_pick in picks_to_add: - wf_id = new_pick.waveform_id - if wf_id in wf_ids_old: - for old_pick in old_picks: - comparison = [old_pick.waveform_id == new_pick.waveform_id, - old_pick.phase_hint == new_pick.phase_hint, - old_pick.method_id == new_pick.method_id] - if all(comparison): - del (old_pick) - old_picks.append(new_pick) - elif not other.isNew() and self.isNew(): - new = other + self - self.evtdata = new.get_evt_data() - elif self.isNew() and other.isNew(): + if not isinstance(other, Data): + raise TypeError("Operands must be of type 'Data'") + if self.event_data.is_new() and other.event_data.is_new(): pass - elif rs_id == rs_id_other: - other.setNew() + elif other.event_data.is_new(): + new_picks = other.event_data.evtdata.picks + old_picks = self.event_data.evtdata.picks + old_picks.extend([pick for pick in new_picks if pick not in old_picks]) + elif self.event_data.is_new(): + return other + self + elif self.event_data.get_id() == other.event_data.get_id(): + other.event_data.set_new() return self + other else: - raise ValueError("both Data objects have differing " - "unique Event identifiers") + raise ValueError("Both Data objects have differing unique Event identifiers") return self - def getPicksStr(self): - """ - Return picks in event data - :return: picks seperated by newlines - :rtype: str - """ - picks_str = '' - for pick in self.get_evt_data().picks: - picks_str += str(pick) + '\n' - return picks_str - - def getParent(self): - """ - Get PySide.QtGui.QWidget parent object - """ + def get_parent(self): return self._parent - def isNew(self): - return self._new + def filter_wf_data(self, **kwargs): + self.waveform_data.wfdata.detrend('linear') + self.waveform_data.wfdata.taper(0.02, type='cosine') + self.waveform_data.wfdata.filter(**kwargs) + self.waveform_data.dirty = True - def setNew(self): - self._new = True + def set_wf_data(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0): + return self.waveform_data.load_waveforms(fnames, fnames_alt, check_rotated, metadata, tstart, tstop) - def getCutTimes(self): - """ - Returns earliest start and latest end of all waveform data - :return: minimum start time and maximum end time as a tuple - :rtype: (UTCDateTime, UTCDateTime) - """ - if self.cuttimes is None: - self.updateCutTimes() - return self.cuttimes + def reset_wf_data(self): + self.waveform_data.reset() - def updateCutTimes(self): - """ - Update cuttimes to contain earliest start and latest end time - of all waveform data - :rtype: None - """ - self.cuttimes = full_range(self.getWFData()) + def get_wf_data(self): + return self.waveform_data.wfdata - def getEventFileName(self): - ID = self.getID() - # handle forbidden filenames especially on windows systems - return fnConstructor(str(ID)) - - def checkEvent(self, event, fcheck, forceOverwrite=False): - """ - Check information in supplied event and own event and replace with own - information if no other information are given or forced by forceOverwrite - :param event: Event that supplies information for comparison - :type event: pylot.core.util.event.Event - :param fcheck: check and delete existing information - can be a str or a list of strings of ['manual', 'auto', 'origin', 'magnitude'] - :type fcheck: str, [str] - :param forceOverwrite: Set to true to force overwrite own information. If false, - supplied information from event is only used if there is no own information in that - category (given in fcheck: manual, auto, origin, magnitude) - :type forceOverwrite: bool - :return: - :rtype: None - """ - if 'origin' in fcheck: - self.replaceOrigin(event, forceOverwrite) - if 'magnitude' in fcheck: - self.replaceMagnitude(event, forceOverwrite) - if 'auto' in fcheck: - self.replacePicks(event, 'auto') - if 'manual' in fcheck: - self.replacePicks(event, 'manual') - - def replaceOrigin(self, event, forceOverwrite=False): - """ - Replace own origin with the one supplied in event if own origin is not - existing or forced by forceOverwrite = True - :param event: Event that supplies information for comparison - :type event: pylot.core.util.event.Event - :param forceOverwrite: always replace own information with supplied one if true - :type forceOverwrite: bool - :return: - :rtype: None - """ - if self.get_evt_data().origins or forceOverwrite: - if event.origins: - print("Found origin, replace it by new origin.") - event.origins = self.get_evt_data().origins - - def replaceMagnitude(self, event, forceOverwrite=False): - """ - Replace own magnitude with the one supplied in event if own magnitude is not - existing or forced by forceOverwrite = True - :param event: Event that supplies information for comparison - :type event: pylot.core.util.event.Event - :param forceOverwrite: always replace own information with supplied one if true - :type forceOverwrite: bool - :return: - :rtype: None - """ - if self.get_evt_data().magnitudes or forceOverwrite: - if event.magnitudes: - print("Found magnitude, replace it by new magnitude") - event.magnitudes = self.get_evt_data().magnitudes - - def replacePicks(self, event, picktype): - """ - Replace picks in event with own picks - :param event: Event that supplies information for comparison - :type event: pylot.core.util.event.Event - :param picktype: 'auto' or 'manual' picks - :type picktype: str - :return: - :rtype: None - """ - checkflag = 1 - picks = event.picks - # remove existing picks - for j, pick in reversed(list(enumerate(picks))): - try: - if picktype in str(pick.method_id.id): - picks.pop(j) - checkflag = 2 - except AttributeError as e: - msg = '{}'.format(e) - print(e) - checkflag = 0 - if checkflag > 0: - if checkflag == 1: - print("Write new %s picks to catalog." % picktype) - if checkflag == 2: - print("Found %s pick(s), remove them and append new picks to catalog." % picktype) - - # append new picks - for pick in self.get_evt_data().picks: - if picktype in str(pick.method_id.id): - picks.append(pick) - - def exportEvent(self, fnout, fnext='.xml', fcheck='auto', upperErrors=None): - """ - Export event to file - :param fnout: basename of file - :param fnext: file extensions xml, cnv, obs, focmec, or/and pha - :param fcheck: check and delete existing information - can be a str or a list of strings of ['manual', 'auto', 'origin', 'magnitude'] - """ - from pylot.core.util.defaults import OUTPUTFORMATS - if not type(fcheck) == list: - fcheck = [fcheck] - - try: - evtformat = OUTPUTFORMATS[fnext] - except KeyError as e: - errmsg = '{0}; selected file extension {1} not ' \ - 'supported'.format(e, fnext) - raise FormatError(errmsg) - - if hasattr(self.get_evt_data(), 'notes'): - try: - with open(os.path.join(os.path.dirname(fnout), 'notes.txt'), 'w') as notes_file: - notes_file.write(self.get_evt_data().notes) - except Exception as e: - print('Warning: Could not save notes.txt: ', str(e)) - - # check for already existing xml-file - if fnext == '.xml': - if os.path.isfile(fnout + fnext): - print("xml-file already exists! Check content ...") - cat = read_events(fnout + fnext) - if len(cat) > 1: - raise IOError('Ambigious event information in file {}'.format(fnout + fnext)) - if len(cat) < 1: - raise IOError('No event information in file {}'.format(fnout + fnext)) - event = cat[0] - if not event.resource_id == self.get_evt_data().resource_id: - QMessageBox.warning(self, 'Warning', 'Different resource IDs!') - return - self.checkEvent(event, fcheck) - self.setEvtData(event) - - self.get_evt_data().write(fnout + fnext, format=evtformat) - - # try exporting event - else: - evtdata_org = self.get_evt_data() - picks = evtdata_org.picks - eventpath = evtdata_org.path - picks_copy = copy.deepcopy(picks) - evtdata_copy = Event(eventpath) - evtdata_copy.picks = picks_copy - - # check for stations picked automatically as well as manually - # Prefer manual picks! - for i in range(len(picks)): - if picks[i].method_id == 'manual': - mstation = picks[i].waveform_id.station_code - mstation_ext = mstation + '_' - for k in range(len(picks_copy)): - if ((picks_copy[k].waveform_id.station_code == mstation) or - (picks_copy[k].waveform_id.station_code == mstation_ext)) and \ - (picks_copy[k].method_id == 'auto'): - del picks_copy[k] - break - lendiff = len(picks) - len(picks_copy) - if lendiff != 0: - print("Manual as well as automatic picks available. Prefered the {} manual ones!".format(lendiff)) - - - no_uncertainties_p = [] - no_uncertainties_s = [] - if upperErrors: - # check for pick uncertainties exceeding adjusted upper errors - # Picks with larger uncertainties will not be saved in output file! - for j in range(len(picks)): - for i in range(len(picks_copy)): - if picks_copy[i].phase_hint[0] == 'P': - # Skipping pick if no upper_uncertainty is found and warning user - if picks_copy[i].time_errors['upper_uncertainty'] is None: - #print("{1} P-Pick of station {0} does not have upper_uncertainty and cant be checked".format( - # picks_copy[i].waveform_id.station_code, - # picks_copy[i].method_id)) - if not picks_copy[i].waveform_id.station_code in no_uncertainties_p: - no_uncertainties_p.append(picks_copy[i].waveform_id.station_code) - continue - - #print ("checking for upper_uncertainty") - if (picks_copy[i].time_errors['uncertainty'] is None) or \ - (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[0]): - print("Uncertainty exceeds or equal adjusted upper time error!") - print("Adjusted uncertainty: {}".format(upperErrors[0])) - print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty'])) - print("{1} P-Pick of station {0} will not be saved in outputfile".format( - picks_copy[i].waveform_id.station_code, - picks_copy[i].method_id)) - del picks_copy[i] - break - if picks_copy[i].phase_hint[0] == 'S': - - # Skipping pick if no upper_uncertainty is found and warning user - if picks_copy[i].time_errors['upper_uncertainty'] is None: - #print("{1} S-Pick of station {0} does not have upper_uncertainty and cant be checked".format( - #picks_copy[i].waveform_id.station_code, - #picks_copy[i].method_id)) - if not picks_copy[i].waveform_id.station_code in no_uncertainties_s: - no_uncertainties_s.append(picks_copy[i].waveform_id.station_code) - continue - - - if (picks_copy[i].time_errors['uncertainty'] is None) or \ - (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[1]): - print("Uncertainty exceeds or equal adjusted upper time error!") - print("Adjusted uncertainty: {}".format(upperErrors[1])) - print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty'])) - print("{1} S-Pick of station {0} will not be saved in outputfile".format( - picks_copy[i].waveform_id.station_code, - picks_copy[i].method_id)) - del picks_copy[i] - break - for s in no_uncertainties_p: - print("P-Pick of station {0} does not have upper_uncertainty and cant be checked".format(s)) - for s in no_uncertainties_s: - print("S-Pick of station {0} does not have upper_uncertainty and cant be checked".format(s)) - - if fnext == '.obs': - try: - evtdata_copy.write(fnout + fnext, format=evtformat) - # write header afterwards - evid = str(evtdata_org.resource_id).split('/')[1] - header = '# EQEVENT: Label: EQ%s Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n' % evid - nllocfile = open(fnout + fnext) - l = nllocfile.readlines() - # Adding A0/Generic Amplitude to .obs file - # l2 = [] - # for li in l: - # for amp in evtdata_org.amplitudes: - # if amp.waveform_id.station_code == li[0:5].strip(): - # li = li[0:64] + '{:0.2e}'.format(amp.generic_amplitude) + li[73:-1] + '\n' - # l2.append(li) - # l = l2 - nllocfile.close() - l.insert(0, header) - nllocfile = open(fnout + fnext, 'w') - nllocfile.write("".join(l)) - nllocfile.close() - except KeyError as e: - raise KeyError('''{0} export format - not implemented: {1}'''.format(evtformat, e)) - if fnext == '.cnv': - try: - velest.export(picks_copy, fnout + fnext, eventinfo=self.get_evt_data()) - except KeyError as e: - raise KeyError('''{0} export format - not implemented: {1}'''.format(evtformat, e)) - if fnext == '_focmec.in': - try: - infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in') - print('Using default input file {}'.format(infile)) - parameter = PylotParameter(infile) - focmec.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data()) - except KeyError as e: - raise KeyError('''{0} export format - not implemented: {1}'''.format(evtformat, e)) - if fnext == '.pha': - try: - infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in') - print('Using default input file {}'.format(infile)) - parameter = PylotParameter(infile) - hypodd.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data()) - except KeyError as e: - raise KeyError('''{0} export format - not implemented: {1}'''.format(evtformat, e)) - - def getComp(self): - """ - Get component (ZNE) - """ - return self.comp - - def getID(self): - """ - Get unique resource id - """ - try: - return self.evtdata.get('resource_id').id - except: - return None - - def filterWFData(self, kwargs): - """ - Filter waveform data - :param kwargs: arguments to pass through to filter function - """ - data = self.getWFData() - data.detrend('linear') - data.taper(0.02, type='cosine') - data.filter(**kwargs) - self.dirty = True - - def setWFData(self, fnames, fnames_syn=None, checkRotated=False, metadata=None, tstart=0, tstop=0): - """ - Clear current waveform data and set given waveform data - :param fnames: waveform data names to append - :type fnames: list - """ - def check_fname_exists(filenames: list) -> list: - if filenames: - filenames = [fn for fn in filenames if os.path.isfile(fn)] - return filenames - - self.wfdata = Stream() - self.wforiginal = None - self.wfsyn = Stream() - if tstart == tstop: - tstart = tstop = None - self.tstart = tstart - self.tstop = tstop - - fnames = check_fname_exists(fnames) - fnames_syn = check_fname_exists(fnames_syn) - # if obspy_dmt: - # wfdir = 'raw' - # self.processed = False - # for fname in fnames: - # if fname.endswith('processed'): - # wfdir = 'processed' - # self.processed = True - # break - # for fpath in fnames: - # if fpath.endswith(wfdir): - # wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)] - # if 'syngine' in fpath.split('/')[-1]: - # wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)] - # else: - # wffnames = fnames - if fnames is not None: - self.appendWFData(fnames) - if fnames_syn is not None: - self.appendWFData(fnames_syn, synthetic=True) - else: - return False - - # various pre-processing steps: - # remove possible underscores in station names - # self.wfdata = remove_underscores(self.wfdata) - # check for gaps and merge - self.wfdata, _ = check_for_gaps_and_merge(self.wfdata) - # check for nans - check_for_nan(self.wfdata) - # check for stations with rotated components - if checkRotated and metadata is not None: - self.wfdata = check4rotated(self.wfdata, metadata, verbosity=0) - # trim station components to same start value - trim_station_components(self.wfdata, trim_start=True, trim_end=False) - - # make a copy of original data - self.wforiginal = self.getWFData().copy() - self.dirty = False - return True - - def appendWFData(self, fnames, synthetic=False): - """ - Read waveform data from fnames and append it to current wf data - :param fnames: waveform data to append - :type fnames: list - """ - assert isinstance(fnames, list), "input parameter 'fnames' is " \ - "supposed to be of type 'list' " \ - "but is actually" \ - " {0}".format(type(fnames)) - if self.dirty: - self.resetWFData() - - real_or_syn_data = {True: self.wfsyn, - False: self.wfdata} - - warnmsg = '' - for fname in set(fnames): - try: - real_or_syn_data[synthetic] += read(fname, starttime=self.tstart, endtime=self.tstop) - except TypeError: - try: - real_or_syn_data[synthetic] += read(fname, format='GSE2', starttime=self.tstart, endtime=self.tstop) - except Exception as e: - try: - real_or_syn_data[synthetic] += read(fname, format='SEGY', starttime=self.tstart, - endtime=self.tstop) - except Exception as e: - warnmsg += '{0}\n{1}\n'.format(fname, e) - except SacIOError as se: - warnmsg += '{0}\n{1}\n'.format(fname, se) - if warnmsg: - warnmsg = 'WARNING in appendWFData: unable to read waveform data\n' + warnmsg - print(warnmsg) - - def getWFData(self): - return self.wfdata - - def getOriginalWFData(self): - return self.wforiginal - - def getSynWFData(self): - return self.wfsyn - - def resetWFData(self): - """ - Set waveform data to original waveform data - """ - if self.getOriginalWFData(): - self.wfdata = self.getOriginalWFData().copy() - else: - self.wfdata = Stream() - self.dirty = False - - def resetPicks(self): - """ - Clear all picks from event - """ - self.get_evt_data().picks = [] - - def get_evt_data(self): - return self.evtdata - - def setEvtData(self, event): - self.evtdata = event - - def applyEVTData(self, data, typ='pick'): - """ - Either takes an `obspy.core.event.Event` object and applies all new - information on the event to the actual data if typ is 'event or - creates ObsPy pick objects and append it to the picks list from the - PyLoT dictionary contain all picks if type is pick - :param data: data to apply, either picks or complete event - :type data: - :param typ: which event data to apply, 'pick' or 'event' - :type typ: str - :param authority_id: (currently unused) - :type: str - :raise OverwriteError: - """ - - def applyPicks(picks): - """ - Creates ObsPy pick objects and append it to the picks list from the - PyLoT dictionary contain all picks. - :param picks: - :raise OverwriteError: raises an OverwriteError if the picks list is - not empty. The GUI will then ask for a decision. - """ - # firstonset = find_firstonset(picks) - # check for automatic picks - print("Writing phases to ObsPy-quakeml file") - for key in picks: - if not picks[key].get('P'): - continue - if picks[key]['P']['picker'] == 'auto': - print("Existing auto-picks will be overwritten in pick-dictionary!") - picks = picks_from_picksdict(picks) - break - else: - if self.get_evt_data().picks: - raise OverwriteError('Existing picks would be overwritten!') - else: - picks = picks_from_picksdict(picks) - break - self.get_evt_data().picks = picks - # if 'smi:local' in self.getID() and firstonset: - # fonset_str = firstonset.strftime('%Y_%m_%d_%H_%M_%S') - # ID = ResourceIdentifier('event/' + fonset_str) - # ID.convertIDToQuakeMLURI(authority_id=authority_id) - # self.get_evt_data().resource_id = ID - - def applyEvent(event): - """ - takes an `obspy.core.event.Event` object and applies all new - information on the event to the actual data - :param event: - """ - if event is None: - print("applyEvent: Received None") - return - if self.isNew(): - self.setEvtData(event) - else: - # prevent overwriting original pick information - event_old = self.get_evt_data() - if not event_old.resource_id == event.resource_id: - print("WARNING: Missmatch in event resource id's: {} and {}".format( - event_old.resource_id, - event.resource_id)) - else: - picks = copy.deepcopy(event_old.picks) - event = merge_picks(event, picks) - # apply event information from location - event_old.update(event) - - applydata = {'pick': applyPicks, - 'event': applyEvent} - - applydata[typ](data) - self._new = False + def rotate_wf_data(self): + self.waveform_data.rotate_zne(self.metadata) class GenericDataStructure(object): @@ -837,22 +248,6 @@ class PilotDataStructure(GenericDataStructure): self.setExpandFields(['root', 'database']) -class ObspyDMTdataStructure(GenericDataStructure): - """ - Object containing the data access information for the old PILOT data - structure. - """ - - def __init__(self, **fields): - if not fields: - fields = {'database': '', - 'root': ''} - - GenericDataStructure.__init__(self, **fields) - - self.setExpandFields(['root', 'database']) - - class SeiscompDataStructure(GenericDataStructure): """ Dictionary containing the data access information for an SDS data archive: diff --git a/pylot/core/io/event.py b/pylot/core/io/event.py new file mode 100644 index 00000000..2f7a669b --- /dev/null +++ b/pylot/core/io/event.py @@ -0,0 +1,106 @@ +import copy +from dataclasses import dataclass +from typing import Union + +from obspy import read_events +from obspy.core.event import Event as ObsPyEvent + + +@dataclass +class EventData: + evtdata: Union[ObsPyEvent, None] = None + _new: bool = False + + def __init__(self, evtdata=None): + self.set_event_data(evtdata) + + def set_event_data(self, evtdata): + if isinstance(evtdata, ObsPyEvent): + self.evtdata = evtdata + elif isinstance(evtdata, dict): + self.evtdata = self.read_pilot_event(**evtdata) + elif isinstance(evtdata, str): + self.evtdata = self.read_event_file(evtdata) + else: + self.set_new() + self.evtdata = ObsPyEvent(picks=[]) + + def read_event_file(self, evtdata: str) -> ObsPyEvent: + try: + cat = read_events(evtdata) + if len(cat) != 1: + raise ValueError(f'ambiguous event information for file: {evtdata}') + return cat[0] + except TypeError as e: + self.handle_event_file_error(e, evtdata) + + def handle_event_file_error(self, e: TypeError, evtdata: str): + if 'Unknown format for file' in str(e): + if 'PHASES' in evtdata: + picks = self.picksdict_from_pilot(evtdata) + evtdata = ObsPyEvent(picks=self.picks_from_picksdict(picks)) + elif 'LOC' in evtdata: + raise NotImplementedError('PILOT location information read support not yet implemented.') + elif 'event.pkl' in evtdata: + evtdata = self.qml_from_obspy_dmt(evtdata) + else: + raise e + else: + raise e + + def set_new(self): + self._new = True + + def is_new(self) -> bool: + return self._new + + def get_picks_str(self) -> str: + return '\n'.join(str(pick) for pick in self.evtdata.picks) + + def replace_origin(self, event: ObsPyEvent, force_overwrite: bool = False): + if self.evtdata.origins or force_overwrite: + event.origins = self.evtdata.origins + + def replace_magnitude(self, event: ObsPyEvent, force_overwrite: bool = False): + if self.evtdata.magnitudes or force_overwrite: + event.magnitudes = self.evtdata.magnitudes + + def replace_picks(self, event: ObsPyEvent, picktype: str): + checkflag = 1 + picks = event.picks + for j, pick in reversed(list(enumerate(picks))): + if picktype in str(pick.method_id.id): + picks.pop(j) + checkflag = 2 + + if checkflag > 0: + for pick in self.evtdata.picks: + if picktype in str(pick.method_id.id): + picks.append(pick) + + def get_id(self) -> Union[str, None]: + try: + return self.evtdata.resource_id.id + except: + return None + + def apply_event_data(self, data, typ='pick'): + if typ == 'pick': + self.apply_picks(data) + elif typ == 'event': + self.apply_event(data) + + def apply_picks(self, picks): + self.evtdata.picks = picks + + def apply_event(self, event: ObsPyEvent): + if self.is_new(): + self.evtdata = event + else: + old_event = self.evtdata + if old_event.resource_id == event.resource_id: + picks = copy.deepcopy(old_event.picks) + event = self.merge_picks(event, picks) + old_event.update(event) + else: + print(f"WARNING: Mismatch in event resource id's: {old_event.resource_id} and {event.resource_id}") diff --git a/pylot/core/io/location.py b/pylot/core/io/location.py index 6fa2bbac..3ef79519 100644 --- a/pylot/core/io/location.py +++ b/pylot/core/io/location.py @@ -1,7 +1,7 @@ from obspy import UTCDateTime from obspy.core import event as ope -from pylot.core.util.utils import getLogin, getHash +from pylot.core.util.utils import get_login, get_hash def create_amplitude(pickID, amp, unit, category, cinfo): @@ -61,7 +61,7 @@ def create_creation_info(agency_id=None, creation_time=None, author=None): :return: ''' if author is None: - author = getLogin() + author = get_login() if creation_time is None: creation_time = UTCDateTime() return ope.CreationInfo(agency_id=agency_id, author=author, @@ -210,7 +210,7 @@ def create_resourceID(timetohash, restype, authority_id=None, hrstr=None): ''' assert isinstance(timetohash, UTCDateTime), "'timetohash' is not an ObsPy" \ "UTCDateTime object" - hid = getHash(timetohash) + hid = get_hash(timetohash) if hrstr is None: resID = ope.ResourceIdentifier(restype + '/' + hid[0:6]) else: diff --git a/pylot/core/io/phases.py b/pylot/core/io/phases.py index 97e3e4c1..1aa4ee89 100644 --- a/pylot/core/io/phases.py +++ b/pylot/core/io/phases.py @@ -1,6 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- import glob +import logging import os import warnings @@ -16,29 +17,10 @@ from pylot.core.io.inputs import PylotParameter from pylot.core.io.location import create_event, \ create_magnitude from pylot.core.pick.utils import select_for_phase, get_quality_class -from pylot.core.util.utils import getOwner, full_range, four_digits, transformFilterString4Export, \ +from pylot.core.util.utils import get_owner, full_range, four_digits, transformFilterString4Export, \ backtransformFilterString, loopIdentifyPhase, identifyPhase -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 @@ -58,7 +40,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs): if phasfn is not None and os.path.isfile(phasfn): phases = sio.loadmat(phasfn) phasctime = UTCDateTime(os.path.getmtime(phasfn)) - phasauthor = getOwner(phasfn) + phasauthor = get_owner(phasfn) else: phases = None phasctime = None @@ -66,7 +48,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs): if locfn is not None and os.path.isfile(locfn): loc = sio.loadmat(locfn) locctime = UTCDateTime(os.path.getmtime(locfn)) - locauthor = getOwner(locfn) + locauthor = get_owner(locfn) else: loc = None locctime = None @@ -192,32 +174,7 @@ def convert_pilot_times(time_array): return UTCDateTime(*times) -def picksdict_from_obs(fn): - """ - create pick dictionary from obs file - :param fn: filename - :type fn: - :return: - :rtype: - """ - picks = dict() - station_name = str() - for line in open(fn, 'r'): - if line.startswith('#'): - continue - else: - phase_line = line.split() - if not station_name == phase_line[0]: - phase = dict() - station_name = phase_line[0] - phase_name = phase_line[4].upper() - pick = UTCDateTime(phase_line[6] + phase_line[7] + phase_line[8]) - phase[phase_name] = dict(mpp=pick, fm=phase_line[5]) - picks[station_name] = phase - return picks - - -def picksdict_from_picks(evt): +def picksdict_from_picks(evt, parameter=None): """ Takes an Event object and return the pick dictionary commonly used within PyLoT @@ -230,6 +187,7 @@ def picksdict_from_picks(evt): 'auto': {} } for pick in evt.picks: + errors = None phase = {} station = pick.waveform_id.station_code if pick.waveform_id.channel_code is None: @@ -273,32 +231,28 @@ def picksdict_from_picks(evt): phase['epp'] = epp phase['lpp'] = lpp phase['spe'] = spe - try: - phase['weight'] = weight - except: - # get onset weight from uncertainty - infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in') - print('Using default input file {}'.format(infile)) - parameter = PylotParameter(infile) + weight = phase.get('weight') + if not weight: + if not parameter: + logging.warning('Using ') + logging.warning('Using default input parameter') + parameter = PylotParameter() pick.phase_hint = identifyPhase(pick.phase_hint) if pick.phase_hint == 'P': errors = parameter['timeerrorsP'] elif pick.phase_hint == 'S': errors = parameter['timeerrorsS'] - weight = get_quality_class(spe, errors) - phase['weight'] = weight + if errors: + weight = get_quality_class(spe, errors) + phase['weight'] = weight phase['channel'] = channel phase['network'] = network phase['picker'] = pick_method - try: - if pick.polarity == 'positive': - phase['fm'] = 'U' - elif pick.polarity == 'negative': - phase['fm'] = 'D' - else: - phase['fm'] = 'N' - except: - print("No FM info available!") + if pick.polarity == 'positive': + phase['fm'] = 'U' + elif pick.polarity == 'negative': + phase['fm'] = 'D' + else: phase['fm'] = 'N' phase['filter_id'] = filter_id if filter_id is not None else '' @@ -375,636 +329,228 @@ def picks_from_picksdict(picks, creation_info=None): return picks_list -def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0): - # TODO: change root to datapath - db_root = os.path.join(root_dir, db_dir) - evt_list = glob.glob1(db_root, 'e????.???.??') - - for evt in evt_list: - if verbosity > 0: - print('Reassessing event {0}'.format(evt)) - reassess_pilot_event(root_dir, db_dir, evt, out_dir, fn_param, verbosity) - - -def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None, verbosity=0): - from obspy import read - - from pylot.core.io.inputs import PylotParameter - from pylot.core.pick.utils import earllatepicker - # TODO: change root to datapath - - default = PylotParameter(fn_param, verbosity) - - search_base = os.path.join(root_dir, db_dir, event_id) - phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat')) - if not phases_file: - return - if verbosity > 1: - print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0])) - picks_dict = picksdict_from_pilot(phases_file[0]) - if verbosity > 0: - print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict)) - datacheck = list() - info = None - for station in picks_dict.keys(): - fn_pattern = os.path.join(search_base, '{0}*'.format(station)) - try: - st = read(fn_pattern) - except TypeError as e: - if 'Unknown format for file' in e.message: - try: - st = read(fn_pattern, format='GSE2') - except ValueError as e: - if e.message == 'second must be in 0..59': - info = 'A known Error was raised. Please find the list of corrupted files and double-check these files.' - datacheck.append(fn_pattern + ' (time info)\n') - continue - else: - raise ValueError(e.message) - except Exception as e: - if 'No file matching file pattern:' in e.message: - if verbosity > 0: - warnings.warn('no waveform data found for station {station}'.format(station=station), - RuntimeWarning) - datacheck.append(fn_pattern + ' (no data)\n') - continue - else: - raise e - else: - raise e - for phase in picks_dict[station].keys(): - try: - mpp = picks_dict[station][phase]['mpp'] - except KeyError as e: - print(e.message, station) - continue - sel_st = select_for_phase(st, phase) - if not sel_st: - msg = 'no waveform data found for station {station}'.format(station=station) - warnings.warn(msg, RuntimeWarning) - continue - stime, etime = full_range(sel_st) - rel_pick = mpp - stime - epp, lpp, spe = earllatepicker(sel_st, - default.get('nfac{0}'.format(phase)), - default.get('tsnrz' if phase == 'P' else 'tsnrh'), - Pick1=rel_pick, - iplot=0, - verbosity=0) - if epp is None or lpp is None: - continue - epp = stime + epp - lpp = stime + lpp - min_diff = 3 * st[0].stats.delta - if lpp - mpp < min_diff: - lpp = mpp + min_diff - if mpp - epp < min_diff: - epp = mpp - min_diff - picks_dict[station][phase] = dict(epp=epp, mpp=mpp, lpp=lpp, spe=spe) - if datacheck: - if info: - if verbosity > 0: - print(info + ': {0}'.format(search_base)) - fncheck = open(os.path.join(search_base, 'datacheck_list'), 'w') - fncheck.writelines(datacheck) - fncheck.close() - del datacheck - # create Event object for export - evt = ope.Event(resource_id=event_id) - evt.picks = picks_from_picksdict(picks_dict) - # write phase information to file - if not out_dir: - fnout_prefix = os.path.join(root_dir, db_dir, event_id, 'PyLoT_{0}.'.format(event_id)) - else: - out_dir = os.path.join(out_dir, db_dir) - if not os.path.isdir(out_dir): - os.makedirs(out_dir) - fnout_prefix = os.path.join(out_dir, 'PyLoT_{0}.'.format(event_id)) - evt.write(fnout_prefix + 'xml', format='QUAKEML') - - -def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None): +def write_phases(arrivals, fformat, filename, parameter=None, eventinfo=None): """ - Function of methods to write phases to the following standard file - formats used for locating earthquakes: + Writes earthquake phase data to different file formats. - HYPO71, NLLoc, VELEST, HYPOSAT, FOCMEC, and hypoDD - - :param arrivals:dictionary containing all phase information including - station ID, phase, first motion, weight (uncertainty), ... + :param arrivals: Dictionary containing phase information (station ID, phase, first motion, weight, etc.) :type arrivals: dict - - :param fformat: chosen file format (location routine), - choose between NLLoc, HYPO71, HYPOSAT, VELEST, - HYPOINVERSE, FOCMEC, and hypoDD + :param fformat: File format to write to (e.g., 'NLLoc', 'HYPO71', 'HYPOSAT', 'VELEST', 'HYPODD', 'FOCMEC') :type fformat: str - - :param filename: full path and name of phase file - :type filename: string - - :param parameter: all input information - :type parameter: object - - :param eventinfo: optional, needed for VELEST-cnv file - and FOCMEC- and HASH-input files - :type eventinfo: `obspy.core.event.Event` object + :param filename: Path and name of the output phase file + :type filename: str + :param parameter: Additional parameters for writing the phase data + :type parameter: object + :param eventinfo: Event information needed for specific formats like VELEST, FOCMEC, and HASH + :type eventinfo: obspy.core.event.Event """ - if fformat == 'NLLoc': - print("Writing phases to %s for NLLoc" % filename) - fid = open("%s" % filename, 'w') - # write header - fid.write('# EQEVENT: %s Label: EQ%s Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n' % - (parameter.get('database'), parameter.get('eventID'))) - arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere - for key in arrivals: - # P onsets - if 'P' in arrivals[key]: - try: - fm = arrivals[key]['P']['fm'] - except KeyError as e: - print(e) - fm = None - if fm is None: - fm = '?' - onset = arrivals[key]['P']['mpp'] - year = onset.year - month = onset.month - day = onset.day - hh = onset.hour - mm = onset.minute - ss = onset.second - ms = onset.microsecond - ss_ms = ss + ms / 1000000.0 - pweight = 1 # use pick - try: - if arrivals[key]['P']['weight'] >= 4: - pweight = 0 # do not use pick - print("Station {}: Uncertain pick, do not use it!".format(key)) - except KeyError as e: - print(e.message + '; no weight set during processing') - fid.write('%s ? ? ? P %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 %d \n' % (key, - fm, - year, - month, - day, - hh, - mm, - ss_ms, - pweight)) - # S onsets - if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None: - fm = '?' - onset = arrivals[key]['S']['mpp'] - year = onset.year - month = onset.month - day = onset.day - hh = onset.hour - mm = onset.minute - ss = onset.second - ms = onset.microsecond - ss_ms = ss + ms / 1000000.0 - sweight = 1 # use pick - try: - if arrivals[key]['S']['weight'] >= 4: - sweight = 0 # do not use pick - except KeyError as e: - print(str(e) + '; no weight set during processing') - Ao = arrivals[key]['S']['Ao'] # peak-to-peak amplitude - if Ao == None: - Ao = 0.0 - # fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 %d \n' % (key, - fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 %9.2f 0 0 %d \n' % (key, - fm, - year, - month, - day, - hh, - mm, - ss_ms, - Ao, - sweight)) - fid.close() - elif fformat == 'HYPO71': - print("Writing phases to %s for HYPO71" % filename) - fid = open("%s" % filename, 'w') - # write header - fid.write(' %s\n' % - parameter.get('eventID')) - arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere - for key in arrivals: - if arrivals[key]['P']['weight'] < 4: - stat = key - if len(stat) > 4: # HYPO71 handles only 4-string station IDs - stat = stat[1:5] - Ponset = arrivals[key]['P']['mpp'] - Sonset = arrivals[key]['S']['mpp'] - pweight = arrivals[key]['P']['weight'] - sweight = arrivals[key]['S']['weight'] - fm = arrivals[key]['P']['fm'] - if fm is None: - fm = '-' - Ao = arrivals[key]['S']['Ao'] - if Ao is None: - Ao = '' - else: - Ao = str('%7.2f' % Ao) - year = Ponset.year - if year >= 2000: - year = year - 2000 - else: - year = year - 1900 - month = Ponset.month - day = Ponset.day - hh = Ponset.hour - mm = Ponset.minute - ss = Ponset.second - ms = Ponset.microsecond - ss_ms = ss + ms / 1000000.0 - if pweight < 2: - pstr = 'I' - elif pweight >= 2: - pstr = 'E' - if arrivals[key]['S']['weight'] < 4: - Sss = Sonset.second - Sms = Sonset.microsecond - Sss_ms = Sss + Sms / 1000000.0 - Sss_ms = str('%5.02f' % Sss_ms) - if sweight < 2: - sstr = 'I' - elif sweight >= 2: - sstr = 'E' - fid.write('%-4s%sP%s%d %02d%02d%02d%02d%02d%5.2f %s%sS %d %s\n' % (stat, - pstr, - fm, - pweight, - year, - month, - day, - hh, - mm, - ss_ms, - Sss_ms, - sstr, - sweight, - Ao)) - else: - fid.write('%-4s%sP%s%d %02d%02d%02d%02d%02d%5.2f %s\n' % (stat, - pstr, - fm, - pweight, - year, - month, - day, - hh, - mm, - ss_ms, - Ao)) + def write_nlloc(): + with open(filename, 'w') as fid: + fid.write('# EQEVENT: {} Label: EQ{} Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n'.format( + parameter.get('database'), parameter.get('eventID'))) + for key, value in arrivals.items(): + for phase in ['P', 'S']: + if phase in value: + fm = value[phase].get('fm', '?') + onset = value[phase]['mpp'] + ss_ms = onset.second + onset.microsecond / 1000000.0 + weight = 1 if value[phase].get('weight', 0) < 4 else 0 + amp = value[phase].get('Ao', 0.0) if phase == 'S' else '' + fid.write('{} ? ? ? {} {}{}{} {}{} {:7.4f} GAU 0 {} 0 0 {}\n'.format( + key, phase, fm, onset.year, onset.month, onset.day, onset.hour, onset.minute, ss_ms, amp, + weight)) - fid.close() + def write_hypo71(): + with open(filename, 'w') as fid: + fid.write( + ' {}\n'.format(parameter.get('eventID'))) + for key, value in arrivals.items(): + if value['P'].get('weight', 0) < 4: + stat = key[:4] + Ponset = value['P']['mpp'] + Sonset = value.get('S', {}).get('mpp') + pweight = value['P'].get('weight', 0) + sweight = value.get('S', {}).get('weight', 0) + fm = value['P'].get('fm', '-') + Ao = value.get('S', {}).get('Ao', '') + year = Ponset.year - 2000 if Ponset.year >= 2000 else Ponset.year - 1900 + ss_ms = Ponset.second + Ponset.microsecond / 1000000.0 + if Sonset: + Sss_ms = Sonset.second + Sonset.microsecond / 1000000.0 + fid.write('{}P{}{}{} {}{}{}{}{} {:5.2f} {}{}S {} {}\n'.format( + stat, 'I' if pweight < 2 else 'E', fm, pweight, year, Ponset.month, Ponset.day, + Ponset.hour, Ponset.minute, ss_ms, Sss_ms, 'I' if sweight < 2 else 'E', sweight, Ao)) + else: + fid.write('{}P{}{}{} {}{}{}{}{} {:5.2f} {}\n'.format( + stat, 'I' if pweight < 2 else 'E', fm, pweight, year, Ponset.month, Ponset.day, + Ponset.hour, Ponset.minute, ss_ms, Ao)) - elif fformat == 'HYPOSAT': - print("Writing phases to %s for HYPOSAT" % filename) - fid = open("%s" % filename, 'w') - # write header - fid.write('%s, event %s \n' % (parameter.get('database'), parameter.get('eventID'))) - arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere - for key in arrivals: - # P onsets - if 'P' in arrivals[key] and arrivals[key]['P']['mpp'] is not None: - if arrivals[key]['P']['weight'] < 4: - Ponset = arrivals[key]['P']['mpp'] - pyear = Ponset.year - pmonth = Ponset.month - pday = Ponset.day - phh = Ponset.hour - pmm = Ponset.minute - pss = Ponset.second - pms = Ponset.microsecond - Pss = pss + pms / 1000000.0 - # use symmetrized picking error as std - # (read the HYPOSAT manual) - pstd = arrivals[key]['P']['spe'] - if pstd is None: - errorsP = parameter.get('timeerrorsP') - if arrivals[key]['P']['weight'] == 0: - pstd = errorsP[0] - elif arrivals[key]['P']['weight'] == 1: - pstd = errorsP[1] - elif arrivals[key]['P']['weight'] == 2: - pstd = errorsP[2] - elif arrivals[key]['P']['weight'] == 3: - psrd = errorsP[3] - else: - pstd = errorsP[4] - fid.write('%-5s P1 %4.0f %02d %02d %02d %02d %05.02f %5.3f -999. 0.00 -999. 0.00\n' - % (key, pyear, pmonth, pday, phh, pmm, Pss, pstd)) - # S onsets - if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None: - if arrivals[key]['S']['weight'] < 4: - Sonset = arrivals[key]['S']['mpp'] - syear = Sonset.year - smonth = Sonset.month - sday = Sonset.day - shh = Sonset.hour - smm = Sonset.minute - sss = Sonset.second - sms = Sonset.microsecond - Sss = sss + sms / 1000000.0 - sstd = arrivals[key]['S']['spe'] - if pstd is None: - errorsS = parameter.get('timeerrorsS') - if arrivals[key]['S']['weight'] == 0: - pstd = errorsS[0] - elif arrivals[key]['S']['weight'] == 1: - pstd = errorsS[1] - elif arrivals[key]['S']['weight'] == 2: - pstd = errorsS[2] - elif arrivals[key]['S']['weight'] == 3: - psrd = errorsS[3] - else: - pstd = errorsP[4] - fid.write('%-5s S1 %4.0f %02d %02d %02d %02d %05.02f %5.3f -999. 0.00 -999. 0.00\n' - % (key, syear, smonth, sday, shh, smm, Sss, sstd)) - fid.close() + def write_hyposat(): + with open(filename, 'w') as fid: + fid.write('{}, event {} \n'.format(parameter.get('database'), parameter.get('eventID'))) + for key, value in arrivals.items(): + for phase in ['P', 'S']: + if phase in value and value[phase].get('weight', 0) < 4: + onset = value[phase]['mpp'] + ss_ms = onset.second + onset.microsecond / 1000000.0 + std = value[phase].get('spe', parameter.get('timeerrorsP')[value[phase].get('weight', 0)]) + fid.write( + '{:<5} {}1 {:4} {:02} {:02} {:02} {:02} {:05.02f} {:5.3f} -999. 0.00 -999. 0.00\n'.format( + key, phase, onset.year, onset.month, onset.day, onset.hour, onset.minute, ss_ms, std)) - elif fformat == 'VELEST': - print("Writing phases to %s for VELEST" % filename) - fid = open("%s" % filename, 'w') - # get informations needed in cnv-file - # check, whether latitude is N or S and longitude is E or W - try: - eventsource = eventinfo.origins[0] - except: + def write_velest(): + if not eventinfo: print("No source origin calculated yet, thus no cnv-file creation possible!") return - if eventsource['latitude'] < 0: - cns = 'S' - else: - cns = 'N' - if eventsource['longitude'] < 0: - cew = 'W' - else: - cew = 'E' - # get last two integers of origin year - stime = eventsource['time'] - if stime.year - 2000 >= 0: - syear = stime.year - 2000 - else: - syear = stime.year - 1900 - ifx = 0 # default value, see VELEST manual, pp. 22-23 - # write header - fid.write('%s%02d%02d %02d%02d %05.2f %7.4f%c %8.4f%c %7.2f %6.2f %02.0f 0.0 0.03 1.0 1.0\n' % ( - syear, stime.month, stime.day, stime.hour, stime.minute, stime.second, eventsource['latitude'], - cns, eventsource['longitude'], cew, eventsource['depth'], eventinfo.magnitudes[0]['mag'], ifx)) - n = 0 - # check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT) - if isinstance(arrivals, dict) == False: - # convert pick object (PyLoT) into dictionary - evt = ope.Event(resource_id=eventinfo['resource_id']) - evt.picks = arrivals - arrivals = picksdict_from_picks(evt) - # check for automatic and manual picks - # prefer manual picks - usedarrivals = chooseArrivals(arrivals) - for key in usedarrivals: - # P onsets - if 'P' in usedarrivals[key]: - if usedarrivals[key]['P']['weight'] < 4: - n += 1 - stat = key - if len(stat) > 4: # VELEST handles only 4-string station IDs - stat = stat[1:5] - Ponset = usedarrivals[key]['P']['mpp'] - Pweight = usedarrivals[key]['P']['weight'] - Prt = Ponset - stime # onset time relative to source time - if n % 6 != 0: - fid.write('%-4sP%d%6.2f' % (stat, Pweight, Prt)) - else: - fid.write('%-4sP%d%6.2f\n' % (stat, Pweight, Prt)) - # S onsets - if 'S' in usedarrivals[key]: - if usedarrivals[key]['S']['weight'] < 4: - n += 1 - stat = key - if len(stat) > 4: # VELEST handles only 4-string station IDs - stat = stat[1:5] - Sonset = usedarrivals[key]['S']['mpp'] - Sweight = usedarrivals[key]['S']['weight'] - Srt = Ponset - stime # onset time relative to source time - if n % 6 != 0: - fid.write('%-4sS%d%6.2f' % (stat, Sweight, Srt)) - else: - fid.write('%-4sS%d%6.2f\n' % (stat, Sweight, Srt)) - fid.close() + with open(filename, 'w') as fid: + origin = eventinfo.origins[0] + lat_dir = 'S' if origin.latitude < 0 else 'N' + lon_dir = 'W' if origin.longitude < 0 else 'E' + year = origin.time.year - 2000 if origin.time.year >= 2000 else origin.time.year - 1900 + fid.write( + '{}{}{} {}{} {} {:05.2f} {:7.4f}{} {:8.4f}{} {:7.2f} {:6.2f} {:02.0f} 0.0 0.03 1.0 1.0\n'.format( + year, origin.time.month, origin.time.day, origin.time.hour, origin.time.minute, origin.time.second, + origin.latitude, lat_dir, origin.longitude, lon_dir, origin.depth, eventinfo.magnitudes[0].mag, 0)) + for key, value in arrivals.items(): + for phase in ['P', 'S']: + if phase in value and value[phase].get('weight', 0) < 4: + onset = value[phase]['mpp'] + rt = (onset - origin.time).total_seconds() + fid.write('{:<4}{}{}{:6.2f}\n'.format(key[:4], phase, value[phase].get('weight', 0), rt)) - elif fformat == 'HYPODD': - print("Writing phases to %s for hypoDD" % filename) - fid = open("%s" % filename, 'w') - # get event information needed for hypoDD-phase file - try: - eventsource = eventinfo.origins[0] - except: + def write_hypodd(): + if not eventinfo: print("No source origin calculated yet, thus no hypoDD-infile creation possible!") return - stime = eventsource['time'] - try: - event = eventinfo['pylot_id'] - hddID = event.split('.')[0][1:5] - except: - print("Error 1111111!") - hddID = "00000" - # write header - fid.write('# %d %d %d %d %d %5.2f %7.4f +%6.4f %7.4f %4.2f 0.1 0.5 %4.2f %s\n' % ( - stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second, - eventsource['latitude'], eventsource['longitude'], eventsource['depth'] / 1000, - eventinfo.magnitudes[0]['mag'], eventsource['quality']['standard_error'], hddID)) - # check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT) - if isinstance(arrivals, dict) == False: - # convert pick object (PyLoT) into dictionary - evt = ope.Event(resource_id=eventinfo['resource_id']) - evt.picks = arrivals - arrivals = picksdict_from_picks(evt) - # check for automatic and manual picks - # prefer manual picks - usedarrivals = chooseArrivals(arrivals) - for key in usedarrivals: - if 'P' in usedarrivals[key]: - # P onsets - if usedarrivals[key]['P']['weight'] < 4: - Ponset = usedarrivals[key]['P']['mpp'] - Prt = Ponset - stime # onset time relative to source time - fid.write('%s %6.3f 1 P\n' % (key, Prt)) - if 'S' in usedarrivals[key]: - # S onsets - if usedarrivals[key]['S']['weight'] < 4: - Sonset = usedarrivals[key]['S']['mpp'] - Srt = Sonset - stime # onset time relative to source time - fid.write('%-5s %6.3f 1 S\n' % (key, Srt)) + with open(filename, 'w') as fid: + origin = eventinfo.origins[0] + stime = origin.time + fid.write('# {} {} {} {} {} {} {:7.4f} +{:6.4f} {:7.4f} {:4.2f} 0.1 0.5 {:4.2f} {}\n'.format( + stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second, + origin.latitude, origin.longitude, origin.depth / 1000, eventinfo.magnitudes[0].mag, + origin.quality.standard_error, "00000")) + for key, value in arrivals.items(): + for phase in ['P', 'S']: + if phase in value and value[phase].get('weight', 0) < 4: + onset = value[phase]['mpp'] + rt = (onset - stime).total_seconds() + fid.write('{} {:6.3f} 1 {}\n'.format(key, rt, phase)) - fid.close() - - elif fformat == 'FOCMEC': - print("Writing phases to %s for FOCMEC" % filename) - fid = open("%s" % filename, 'w') - # get event information needed for FOCMEC-input file - try: - eventsource = eventinfo.origins[0] - except: + def write_focmec(): + if not eventinfo: print("No source origin calculated yet, thus no FOCMEC-infile creation possible!") return - stime = eventsource['time'] - - # avoid printing '*' in focmec-input file - if parameter.get('eventid') == '*' or parameter.get('eventid') is None: - evID = 'e0000' - else: - evID = parameter.get('eventid') - - # write header line including event information - fid.write('%s %d%02d%02d%02d%02d%02.0f %7.4f %6.4f %3.1f %3.1f\n' % (evID, - stime.year, stime.month, stime.day, - stime.hour, stime.minute, stime.second, - eventsource['latitude'], - eventsource['longitude'], - eventsource['depth'] / 1000, - eventinfo.magnitudes[0]['mag'])) - picks = eventinfo.picks - # check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT) - if isinstance(arrivals, dict) == False: - # convert pick object (PyLoT) into dictionary - evt = ope.Event(resource_id=eventinfo['resource_id']) - evt.picks = arrivals - arrivals = picksdict_from_picks(evt) - # check for automatic and manual picks - # prefer manual picks - usedarrivals = chooseArrivals(arrivals) - for key in usedarrivals: - if 'P' in usedarrivals[key]: - if usedarrivals[key]['P']['weight'] < 4 and usedarrivals[key]['P']['fm'] is not None: - stat = key - for i in range(len(picks)): - station = picks[i].waveform_id.station_code - if station == stat: - # get resource ID - resid_picks = picks[i].get('resource_id') - # find same ID in eventinfo - # there it is the pick_id!! - for j in range(len(eventinfo.origins[0].arrivals)): - resid_eventinfo = eventinfo.origins[0].arrivals[j].get('pick_id') - if resid_eventinfo == resid_picks and eventinfo.origins[0].arrivals[j].phase == 'P': - if len(stat) > 4: # FOCMEC handles only 4-string station IDs - stat = stat[1:5] - az = eventinfo.origins[0].arrivals[j].get('azimuth') - inz = eventinfo.origins[0].arrivals[j].get('takeoff_angle') - fid.write('%-4s %6.2f %6.2f%s \n' % (stat, - az, - inz, - usedarrivals[key]['P']['fm'])) + with open(filename, 'w') as fid: + origin = eventinfo.origins[0] + stime = origin.time + fid.write('{} {}{:02d}{:02d}{:02d}{:02d}{:02.0f} {:7.4f} {:6.4f} {:3.1f} {:3.1f}\n'.format( + parameter.get('eventid', 'e0000'), stime.year, stime.month, stime.day, stime.hour, stime.minute, + stime.second, origin.latitude, origin.longitude, origin.depth / 1000, eventinfo.magnitudes[0].mag)) + for key, value in arrivals.items(): + if 'P' in value and value['P'].get('weight', 0) < 4 and value['P'].get('fm'): + for pick in eventinfo.picks: + if pick.waveform_id.station_code == key: + for arrival in origin.arrivals: + if arrival.pick_id == pick.resource_id and arrival.phase == 'P': + stat = key[:4] + az = arrival.azimuth + inz = arrival.takeoff_angle + fid.write('{:<4} {:6.2f} {:6.2f}{}\n'.format(stat, az, inz, value['P']['fm'])) break - fid.close() + def write_hash(): + # Define filenames for HASH driver 1 and 2 + filename1 = f"{filename}drv1.phase" + filename2 = f"{filename}drv2.phase" + print(f"Writing phases to {filename1} for HASH-driver 1") + print(f"Writing phases to {filename2} for HASH-driver 2") + + # Open files for writing + with open(filename1, 'w') as fid1, open(filename2, 'w') as fid2: + # Get event information needed for HASH-input file + try: + eventsource = eventinfo.origins[0] + except IndexError: + print("No source origin calculated yet, thus no cnv-file creation possible!") + return + + event = parameter.get('eventID') + hashID = event.split('.')[0][1:5] + latdeg = eventsource['latitude'] + latmin = (eventsource['latitude'] * 60) / 10000 + londeg = eventsource['longitude'] + lonmin = (eventsource['longitude'] * 60) / 10000 + + erh = (eventsource.origin_uncertainty['min_horizontal_uncertainty'] + + eventsource.origin_uncertainty['max_horizontal_uncertainty']) / 2000 + erz = eventsource.depth_errors['uncertainty'] + + stime = eventsource['time'] + syear = stime.year % 100 # Calculate two-digit year + + picks = eventinfo.picks + + # Write header line including event information for HASH-driver 1 + fid1.write(f"{syear:02d}{stime.month:02d}{stime.day:02d}{stime.hour:02d}{stime.minute:02d}" + f"{stime.second:05.2f}{latdeg:2d}N{latmin:05.2f}{londeg:3d}E{lonmin:05.2f}" + f"{eventsource['depth']:6.2f}{eventinfo.magnitudes[0]['mag']:4.2f}{erh:5.2f}{erz:5.2f}{hashID}\n") + + # Write header line including event information for HASH-driver 2 + fid2.write(f"{syear:02d}{stime.month:02d}{stime.day:02d}{stime.hour:02d}{stime.minute:02d}" + f"{stime.second:05.2f}{latdeg}N{latmin:05.2f}{londeg}E{lonmin:6.2f}{eventsource['depth']:5.2f}" + f"{eventsource['quality']['used_phase_count']:3d}{erh:5.2f}{erz:5.2f}" + f"{eventinfo.magnitudes[0]['mag']:4.2f}{hashID}\n") + + # Write phase lines + for key, arrival in arrivals.items(): + if 'P' in arrival and arrival['P']['weight'] < 4 and arrival['P']['fm'] is not None: + stat = key + ccode = arrival['P']['channel'] + ncode = arrival['P']['network'] + Pqual = 'I' if arrival['P']['weight'] < 2 else 'E' + + for pick in picks: + if pick.waveform_id.station_code == stat: + resid_picks = pick.get('resource_id') + for origin_arrival in eventinfo.origins[0].arrivals: + if (origin_arrival.get('pick_id') == resid_picks and + origin_arrival.phase == 'P'): + if len(stat) > 4: # HASH handles only 4-character station IDs + stat = stat[1:5] + + az = origin_arrival.get('azimuth') + inz = origin_arrival.get('takeoff_angle') + dist = origin_arrival.get('distance') + + # Write phase line for HASH-driver 1 + fid1.write(f"{stat:<4}{Pqual}P{arrival['P']['fm']}{arrival['P']['weight']:d}" + f"{dist:3.1f}{inz:03d}{az:03d}{ccode}\n") + + # Write phase line for HASH-driver 2 + fid2.write(f"{stat:<4} {ncode} {ccode} {Pqual} {arrival['P']['fm']}\n") + break + + fid1.write(f"{'':<36}{hashID}") + + # Prefer Manual Picks over automatic ones if possible + arrivals = chooseArrivals(arrivals) # Function not defined, assumed to exist + + if fformat == 'NLLoc': + write_nlloc() + elif fformat == 'HYPO71': + write_hypo71() + elif fformat == 'HYPOSAT': + write_hyposat() + elif fformat == 'VELEST': + write_velest() + elif fformat == 'HYPODD': + write_hypodd() + elif fformat == 'FOCMEC': + write_focmec() elif fformat == 'HASH': - # two different input files for - # HASH-driver 1 and 2 (see HASH manual!) - filename1 = filename + 'drv1' + '.phase' - filename2 = filename + 'drv2' + '.phase' - print("Writing phases to %s for HASH for HASH-driver 1" % filename1) - fid1 = open("%s" % filename1, 'w') - print("Writing phases to %s for HASH for HASH-driver 2" % filename2) - fid2 = open("%s" % filename2, 'w') - # get event information needed for HASH-input file - try: - eventsource = eventinfo.origins[0] - except: - print("No source origin calculated yet, thus no cnv-file creation possible!") - return - eventsource = eventinfo.origins[0] - event = parameter.get('eventID') - hashID = event.split('.')[0][1:5] - latdeg = eventsource['latitude'] - latmin = eventsource['latitude'] * 60 / 10000 - londeg = eventsource['longitude'] - lonmin = eventsource['longitude'] * 60 / 10000 - erh = 1 / 2 * (eventsource.origin_uncertainty['min_horizontal_uncertainty'] + - eventsource.origin_uncertainty['max_horizontal_uncertainty']) / 1000 - erz = eventsource.depth_errors['uncertainty'] - stime = eventsource['time'] - if stime.year - 2000 >= 0: - syear = stime.year - 2000 - else: - syear = stime.year - 1900 - picks = eventinfo.picks - # write header line including event information - # for HASH-driver 1 - fid1.write('%s%02d%02d%02d%02d%5.2f%2dN%5.2f%3dE%5.2f%6.3f%4.2f%5.2f%5.2f%s\n' % (syear, - stime.month, stime.day, - stime.hour, stime.minute, - stime.second, - latdeg, latmin, londeg, - lonmin, eventsource['depth'], - eventinfo.magnitudes[0][ - 'mag'], erh, erz, - hashID)) - # write header line including event information - # for HASH-driver 2 - fid2.write( - '%d%02d%02d%02d%02d%5.2f%dN%5.2f%3dE%6.2f%5.2f %d %5.2f %5.2f %4.2f %s \n' % ( - syear, stime.month, stime.day, - stime.hour, stime.minute, stime.second, - latdeg, latmin, londeg, lonmin, - eventsource['depth'], - eventsource['quality']['used_phase_count'], - erh, erz, eventinfo.magnitudes[0]['mag'], - hashID)) - # Prefer Manual Picks over automatic ones if possible - arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere - # write phase lines - for key in arrivals: - if 'P' in arrivals[key]: - if arrivals[key]['P']['weight'] < 4 and arrivals[key]['P']['fm'] is not None: - stat = key - ccode = arrivals[key]['P']['channel'] - ncode = arrivals[key]['P']['network'] - - if arrivals[key]['P']['weight'] < 2: - Pqual = 'I' - else: - Pqual = 'E' - - for i in range(len(picks)): - station = picks[i].waveform_id.station_code - if station == stat: - # get resource ID - resid_picks = picks[i].get('resource_id') - # find same ID in eventinfo - # there it is the pick_id!! - for j in range(len(eventinfo.origins[0].arrivals)): - resid_eventinfo = eventinfo.origins[0].arrivals[j].get('pick_id') - if resid_eventinfo == resid_picks and eventinfo.origins[0].arrivals[j].phase == 'P': - if len(stat) > 4: # HASH handles only 4-string station IDs - stat = stat[1:5] - az = eventinfo.origins[0].arrivals[j].get('azimuth') - inz = eventinfo.origins[0].arrivals[j].get('takeoff_angle') - dist = eventinfo.origins[0].arrivals[j].get('distance') - # write phase line for HASH-driver 1 - fid1.write( - '%-4s%sP%s%d 0 %3.1f %03d %03d 2 1 %s\n' % ( - stat, Pqual, arrivals[key]['P']['fm'], arrivals[key]['P']['weight'], - dist, inz, az, ccode)) - # write phase line for HASH-driver 2 - fid2.write('%-4s %s %s %s %s \n' % ( - stat, - ncode, - ccode, - Pqual, - arrivals[key]['P']['fm'])) - break - - fid1.write(' %s' % hashID) - fid1.close() - fid2.close() + write_hash() def chooseArrivals(arrivals): @@ -1124,7 +670,7 @@ def getQualitiesfromxml(path, errorsP, errorsS, plotflag=1, figure=None, verbosi mstation = pick.waveform_id.station_code mstation_ext = mstation + '_' for mpick in arrivals_copy: - phase = identifyPhase(loopIdentifyPhase(pick.phase_hint)) # MP MP catch if this fails? + phase = identifyPhase(loopIdentifyPhase(pick.phase_hint)) # MP MP catch if this fails? if ((mpick.waveform_id.station_code == mstation) or (mpick.waveform_id.station_code == mstation_ext)) and \ (mpick.method_id.id.split('/')[1] == 'auto') and \ diff --git a/pylot/core/io/project.py b/pylot/core/io/project.py new file mode 100644 index 00000000..764b004a --- /dev/null +++ b/pylot/core/io/project.py @@ -0,0 +1,177 @@ +import os + +from pylot.core.util.event import Event + + +class Project(object): + ''' + Pickable class containing information of a PyLoT project, like event lists and file locations. + ''' + + # TODO: remove rootpath + def __init__(self): + self.eventlist = [] + self.location = None + self.rootpath = None + self.datapath = None + self.dirty = False + self.parameter = None + self._table = None + + def add_eventlist(self, eventlist): + ''' + Add events from an eventlist containing paths to event directories. + Will skip existing paths. + ''' + if len(eventlist) == 0: + return + for item in eventlist: + event = Event(item) + event.rootpath = self.parameter['rootpath'] + event.database = self.parameter['database'] + event.datapath = self.parameter['datapath'] + if not event.path in self.getPaths(): + self.eventlist.append(event) + self.setDirty() + else: + print('Skipping event with path {}. Already part of project.'.format(event.path)) + self.eventlist.sort(key=lambda x: x.pylot_id) + self.search_eventfile_info() + + def remove_event(self, event): + self.eventlist.remove(event) + + def remove_event_by_id(self, eventID): + for event in self.eventlist: + if eventID in str(event.resource_id): + self.remove_event(event) + break + + def read_eventfile_info(self, filename, separator=','): + ''' + Try to read event information from file (:param:filename) comparing specific event datetimes. + File structure (each row): event, date, time, magnitude, latitude, longitude, depth + separated by :param:separator each. + ''' + with open(filename, 'r') as infile: + for line in infile.readlines(): + eventID, date, time, mag, lat, lon, depth = line.split(separator)[:7] + # skip first line + try: + day, month, year = date.split('/') + except: + continue + year = int(year) + # hardcoded, if year only consists of 2 digits (e.g. 16 instead of 2016) + if year < 100: + year += 2000 + datetime = '{}-{}-{}T{}'.format(year, month, day, time) + try: + datetime = UTCDateTime(datetime) + except Exception as e: + print(e, datetime, filename) + continue + for event in self.eventlist: + if eventID in str(event.resource_id) or eventID in event.origins: + if event.origins: + origin = event.origins[0] # should have only one origin + if origin.time == datetime: + origin.latitude = float(lat) + origin.longitude = float(lon) + origin.depth = float(depth) + else: + continue + elif not event.origins: + origin = Origin(resource_id=event.resource_id, + time=datetime, latitude=float(lat), + longitude=float(lon), depth=float(depth)) + event.origins.append(origin) + event.magnitudes.append(Magnitude(resource_id=event.resource_id, + mag=float(mag), + mag_type='M')) + break + + def search_eventfile_info(self): + ''' + Search all datapaths in rootpath for filenames with given file extension fext + and try to read event info from it + ''' + datapaths = [] + fext = '.csv' + for event in self.eventlist: + if not event.datapath in datapaths: + datapaths.append(event.datapath) + for datapath in datapaths: + # datapath = os.path.join(self.rootpath, datapath) + if os.path.isdir(datapath): + for filename in os.listdir(datapath): + filename = os.path.join(datapath, filename) + if os.path.isfile(filename) and filename.endswith(fext): + try: + self.read_eventfile_info(filename) + except Exception as e: + print('Failed on reading eventfile info from file {}: {}'.format(filename, e)) + else: + print("Directory %s does not exist!" % datapath) + + def getPaths(self): + ''' + Returns paths (eventlist) of all events saved in the project. + ''' + paths = [] + for event in self.eventlist: + paths.append(event.path) + return paths + + def setDirty(self, value=True): + self.dirty = value + + def getEventFromPath(self, path): + ''' + Search for an event in the project by event path. + ''' + for event in self.eventlist: + if event.path == path: + return event + + def save(self, filename=None): + ''' + Save PyLoT Project to a file. + Can be loaded by using project.load(filename). + ''' + try: + import pickle + except ImportError: + import _pickle as pickle + + if filename: + self.location = filename + else: + filename = self.location + + table = self._table # MP: see below + try: + outfile = open(filename, 'wb') + self._table = [] # MP: Workaround as long as table cannot be saved as part of project + pickle.dump(self, outfile, protocol=pickle.HIGHEST_PROTOCOL) + self.setDirty(False) + self._table = table # MP: see above + return True + except Exception as e: + print('Could not pickle PyLoT project. Reason: {}'.format(e)) + self.setDirty() + self._table = table # MP: see above + return False + + @staticmethod + def load(filename): + ''' + Load project from filename. + ''' + import pickle + infile = open(filename, 'rb') + project = pickle.load(infile) + infile.close() + project.location = filename + print('Loaded %s' % filename) + return project diff --git a/pylot/core/io/utils.py b/pylot/core/io/utils.py new file mode 100644 index 00000000..109ee8a2 --- /dev/null +++ b/pylot/core/io/utils.py @@ -0,0 +1,13 @@ +import os +from typing import List + + +def validate_filenames(filenames: List[str]) -> List[str]: + """ + validate a list of filenames for file abundance + :param filenames: list of possible filenames + :type filenames: List[str] + :return: list of valid filenames + :rtype: List[str] + """ + return [fn for fn in filenames if os.path.isfile(fn)] \ No newline at end of file diff --git a/pylot/core/io/waveformdata.py b/pylot/core/io/waveformdata.py new file mode 100644 index 00000000..c2b0078c --- /dev/null +++ b/pylot/core/io/waveformdata.py @@ -0,0 +1,123 @@ +import logging +from dataclasses import dataclass, field +from typing import Union, List + +from obspy import Stream, read +from obspy.io.sac import SacIOError + +from pylot.core.io.utils import validate_filenames +from pylot.core.util.dataprocessing import Metadata +from pylot.core.util.utils import get_stations, check_for_nan, check4rotated + + +@dataclass +class WaveformData: + wfdata: Stream = field(default_factory=Stream) + wforiginal: Union[Stream, None] = None + wf_alt: Stream = field(default_factory=Stream) + dirty: bool = False + + def load_waveforms(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0): + fn_list = validate_filenames(fnames) + if not fn_list: + logging.warning('No valid filenames given for loading waveforms') + else: + self.clear() + self.add_waveforms(fn_list) + + if fnames_alt is None: + pass + else: + alt_fn_list = validate_filenames(fnames_alt) + if not alt_fn_list: + logging.warning('No valid alternative filenames given for loading waveforms') + else: + self.add_waveforms(alt_fn_list, alternative=True) + + if not fn_list and not alt_fn_list: + logging.error('No filenames or alternative filenames given for loading waveforms') + return False + + self.merge() + self.replace_nan() + if not check_rotated or not metadata: + pass + else: + self.rotate_zne() + self.trim_station_traces() + self.wforiginal = self.wfdata.copy() + self.dirty = False + return True + + def add_waveforms(self, fnames: List[str], alternative: bool = False): + data_stream = self.wf_alt if alternative else self.wfdata + warnmsg = '' + for fname in set(fnames): + try: + data_stream += read(fname) + except TypeError: + try: + data_stream += read(fname, format='GSE2') + except Exception as e: + try: + data_stream += read(fname, format='SEGY') + except Exception as e: + warnmsg += f'{fname}\n{e}\n' + except SacIOError as se: + warnmsg += f'{fname}\n{se}\n' + + if warnmsg: + print(f'WARNING in add_waveforms: unable to read waveform data\n{warnmsg}') + + def clear(self): + self.wfdata = Stream() + self.wforiginal = None + self.wf_alt = Stream() + + def reset(self): + """ + Resets the waveform data to its original state. + """ + if self.wforiginal: + self.wfdata = self.wforiginal.copy() + else: + self.wfdata = Stream() + self.dirty = False + + def merge(self): + """ + check for gaps in Stream and merge if gaps are found + """ + gaps = self.wfdata.get_gaps() + if gaps: + merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps] + self.wfdata.merge(method=1) + logging.info('Merged the following stations because of gaps:') + for station in merged: + logging.info(station) + + def replace_nan(self): + """ + Replace all NaNs in data with 0. (in place) + """ + self.wfdata = check_for_nan(self.wfdata) + + def rotate_zne(self, metadata: Metadata = None): + """ + Check all traces in stream for rotation. If a trace is not in ZNE rotation (last symbol of channel code is numeric) and the trace + is in the metadata with azimuth and dip, rotate it to classical ZNE orientation. + Rotating the traces requires them to be of the same length, so, all traces will be trimmed to a common length as a + side effect. + """ + + self.wfdata = check4rotated(self.wfdata, metadata) + + def trim_station_traces(self): + """ + trim data stream to common time window + """ + + for station in get_stations(self.wfdata): + station_traces = self.wfdata.select(station=station) + station_traces.trim(starttime=max([trace.stats.starttime for trace in station_traces]), + endtime=min([trace.stats.endtime for trace in station_traces])) diff --git a/pylot/core/loc/focmec.py b/pylot/core/loc/focmec.py index 2df845cf..29cb7496 100644 --- a/pylot/core/loc/focmec.py +++ b/pylot/core/loc/focmec.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -25,4 +25,4 @@ def export(picks, fnout, parameter, eventinfo): :type eventinfo: list object ''' # write phases to FOCMEC-phase file - writephases(picks, 'FOCMEC', fnout, parameter, eventinfo) + write_phases(picks, 'FOCMEC', fnout, parameter, eventinfo) diff --git a/pylot/core/loc/hash.py b/pylot/core/loc/hash.py index 5ba18faa..d39e5076 100644 --- a/pylot/core/loc/hash.py +++ b/pylot/core/loc/hash.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -25,4 +25,4 @@ def export(picks, fnout, parameter, eventinfo): :type eventinfo: list object ''' # write phases to HASH-phase file - writephases(picks, 'HASH', fnout, parameter, eventinfo) + write_phases(picks, 'HASH', fnout, parameter, eventinfo) diff --git a/pylot/core/loc/hypo71.py b/pylot/core/loc/hypo71.py index 430e2c14..e96c7bc3 100644 --- a/pylot/core/loc/hypo71.py +++ b/pylot/core/loc/hypo71.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -22,4 +22,4 @@ def export(picks, fnout, parameter): :type parameter: object ''' # write phases to HYPO71-phase file - writephases(picks, 'HYPO71', fnout, parameter) + write_phases(picks, 'HYPO71', fnout, parameter) diff --git a/pylot/core/loc/hypodd.py b/pylot/core/loc/hypodd.py index c194ed90..7fe2a744 100644 --- a/pylot/core/loc/hypodd.py +++ b/pylot/core/loc/hypodd.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -25,4 +25,4 @@ def export(picks, fnout, parameter, eventinfo): :type eventinfo: list object ''' # write phases to hypoDD-phase file - writephases(picks, 'HYPODD', fnout, parameter, eventinfo) + write_phases(picks, 'HYPODD', fnout, parameter, eventinfo) diff --git a/pylot/core/loc/hyposat.py b/pylot/core/loc/hyposat.py index e9c05a87..49cbe098 100644 --- a/pylot/core/loc/hyposat.py +++ b/pylot/core/loc/hyposat.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -22,4 +22,4 @@ def export(picks, fnout, parameter): :type parameter: object ''' # write phases to HYPOSAT-phase file - writephases(picks, 'HYPOSAT', fnout, parameter) + write_phases(picks, 'HYPOSAT', fnout, parameter) diff --git a/pylot/core/loc/nll.py b/pylot/core/loc/nll.py index 658bb368..ad3eb571 100644 --- a/pylot/core/loc/nll.py +++ b/pylot/core/loc/nll.py @@ -7,7 +7,7 @@ import subprocess from obspy import read_events -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.gui import which from pylot.core.util.utils import getPatternLine, runProgram from pylot.core.util.version import get_git_version as _getVersionString @@ -34,7 +34,7 @@ def export(picks, fnout, parameter): :type parameter: object ''' # write phases to NLLoc-phase file - writephases(picks, 'NLLoc', fnout, parameter) + write_phases(picks, 'NLLoc', fnout, parameter) def modify_inputs(ctrfn, root, nllocoutn, phasefn, tttn): diff --git a/pylot/core/loc/velest.py b/pylot/core/loc/velest.py index 8dd5d88f..040bacab 100644 --- a/pylot/core/loc/velest.py +++ b/pylot/core/loc/velest.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from pylot.core.io.phases import writephases +from pylot.core.io.phases import write_phases from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -25,4 +25,4 @@ def export(picks, fnout, eventinfo, parameter=None): :type parameter: object ''' # write phases to VELEST-phase file - writephases(picks, 'VELEST', fnout, parameter, eventinfo) + write_phases(picks, 'VELEST', fnout, parameter, eventinfo) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 6cca5691..70f71658 100644 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -22,7 +22,7 @@ from pylot.core.pick.picker import AICPicker, PragPicker from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \ getSNR, fmpicker, checkPonsets, wadaticheck, get_quality_class, PickingFailedException, MissingTraceException from pylot.core.util.utils import getPatternLine, gen_Pool, \ - get_bool, identifyPhaseID, get_None, correct_iplot + get_bool, identifyPhaseID, get_none, correct_iplot def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None, ncores=0, metadata=None, origin=None): @@ -258,7 +258,7 @@ class AutopickStation(object): self.pickparams = copy.deepcopy(pickparam) self.verbose = verbose self.iplot = correct_iplot(iplot) - self.fig_dict = get_None(fig_dict) + self.fig_dict = get_none(fig_dict) self.metadata = metadata self.origin = origin diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index e2bb8e7a..196ac2ee 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -15,7 +15,7 @@ import numpy as np from obspy.core import Stream, UTCDateTime from scipy.signal import argrelmax -from pylot.core.util.utils import get_bool, get_None, SetChannelComponents +from pylot.core.util.utils import get_bool, get_none, SetChannelComponents def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecolor='k'): @@ -136,7 +136,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol PickError = symmetrize_error(diffti_te, diffti_tl) if iplot > 1: - if get_None(fig) is None: + if get_none(fig) is None: fig = plt.figure() # iplot) plt_flag = 1 fig._tight = True @@ -275,7 +275,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'): try: P1 = np.polyfit(xslope1, xraw[islope1], 1) datafit1 = np.polyval(P1, xslope1) - except Exception as e: + except ValueError as e: print("fmpicker: Problems with data fit! {}".format(e)) print("Skip first motion determination!") return FM @@ -321,7 +321,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'): try: P2 = np.polyfit(xslope2, xfilt[islope2], 1) datafit2 = np.polyval(P2, xslope2) - except Exception as e: + except ValueError as e: emsg = 'fmpicker: polyfit failed: {}'.format(e) print(emsg) return FM @@ -344,7 +344,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'): print("fmpicker: Found polarity %s" % FM) if iplot > 1: - if get_None(fig) is None: + if get_none(fig) is None: fig = plt.figure() # iplot) plt_flag = 1 fig._tight = True @@ -868,7 +868,7 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line returnflag = 0 if iplot > 1: - if get_None(fig) is None: + if get_none(fig) is None: fig = plt.figure() # iplot) plt_flag = 1 fig._tight = True @@ -1213,14 +1213,14 @@ def checkZ4S(X, pick, pickparams, iplot, fig=None, linecolor='k'): t = np.linspace(diff_dict[key], trace.stats.endtime - trace.stats.starttime + diff_dict[key], trace.stats.npts) if i == 0: - if get_None(fig) is None: + if get_none(fig) is None: fig = plt.figure() # self.iplot) ### WHY? MP MP plt_flag = 1 ax1 = fig.add_subplot(3, 1, i + 1) ax = ax1 ax.set_title('CheckZ4S, Station %s' % zdat[0].stats.station) else: - if get_None(fig) is None: + if get_none(fig) is None: fig = plt.figure() # self.iplot) ### WHY? MP MP plt_flag = 1 ax = fig.add_subplot(3, 1, i + 1, sharex=ax1) @@ -1494,7 +1494,7 @@ def getQualityFromUncertainty(uncertainty, Errors): # set initial quality to 4 (worst) and change only if one condition is hit quality = 4 - if get_None(uncertainty) is None: + if get_none(uncertainty) is None: return quality if uncertainty <= Errors[0]: diff --git a/pylot/core/util/array_map.py b/pylot/core/util/array_map.py index 3407ba8a..ea8d0c00 100644 --- a/pylot/core/util/array_map.py +++ b/pylot/core/util/array_map.py @@ -124,8 +124,8 @@ class Array_map(QtWidgets.QWidget): self.cmaps_box = QtWidgets.QComboBox() self.cmaps_box.setMaxVisibleItems(20) [self.cmaps_box.addItem(map_name) for map_name in sorted(plt.colormaps())] - # try to set to hsv as default - self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('hsv')) + # try to set to viridis as default + self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('viridis')) self.top_row.addWidget(QtWidgets.QLabel('Select a phase: ')) self.top_row.addWidget(self.comboBox_phase) @@ -474,17 +474,19 @@ class Array_map(QtWidgets.QWidget): transform=ccrs.PlateCarree(), label='deleted')) def openPickDlg(self, ind): - data = self._parent.get_data().getWFData() + wfdata = self._parent.get_data().get_wf_data() + wfdata_comp = self._parent.get_data().get_wf_dataComp() for index in ind: network, station = self._station_onpick_ids[index].split('.')[:2] pyl_mw = self._parent try: - data = data.select(station=station) - if not data: + wfdata = wfdata.select(station=station) + wfdata_comp = wfdata_comp.select(station=station) + if not wfdata: self._warn('No data for station {}'.format(station)) return pickDlg = PickDlg(self._parent, parameter=self.parameter, - data=data, network=network, station=station, + data=wfdata.copy(), data_compare=wfdata_comp.copy(), network=network, station=station, picks=self._parent.get_current_event().getPick(station), autopicks=self._parent.get_current_event().getAutopick(station), filteroptions=self._parent.filteroptions, metadata=self.metadata, diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index 09d2ba13..0c1cf1d1 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -338,19 +338,6 @@ class Metadata(object): return inv, exc -def time_from_header(header): - """ - Function takes in the second line from a .gse file and takes out the date and time from that line. - :param header: second line from .gse file - :type header: string - :return: a list of integers of form [year, month, day, hour, minute, second, microsecond] - """ - timeline = header.split(' ') - time = timeline[1].split('/') + timeline[2].split(':') - time = time[:-1] + time[-1].split('.') - return [int(t) for t in time] - - def check_time(datetime): """ Function takes in date and time as list and validates it's values by trying to make an UTCDateTime object from it @@ -388,167 +375,6 @@ def check_time(datetime): return False -# TODO: change root to datapath -def get_file_list(root_dir): - """ - Function uses a directorie to get all the *.gse files from it. - :param root_dir: a directorie leading to the .gse files - :type root_dir: string - :return: returns a list of filenames (without path to them) - """ - file_list = glob.glob1(root_dir, '*.gse') - return file_list - - -def checks_station_second(datetime, file): - """ - Function uses the given list to check if the parameter 'second' is set to 60 by mistake - and sets the time correctly if so. Can only correct time if no date change would be necessary. - :param datetime: [year, month, day, hour, minute, second, microsecond] - :return: returns the input with the correct value for second - """ - if datetime[5] == 60: - if datetime[4] == 59: - if datetime[3] == 23: - err_msg = 'Date should be next day. ' \ - 'File not changed: {0}'.format(file) - raise ValueError(err_msg) - else: - datetime[3] += 1 - datetime[4] = 0 - datetime[5] = 0 - else: - datetime[4] += 1 - datetime[5] = 0 - return datetime - - -def make_time_line(line, datetime): - """ - Function takes in the original line from a .gse file and a list of date and - time values to make a new line with corrected date and time. - :param line: second line from .gse file. - :type line: string - :param datetime: list of integers [year, month, day, hour, minute, second, microsecond] - :type datetime: list - :return: returns a string to write it into a file. - """ - ins_form = '{0:02d}:{1:02d}:{2:02d}.{3:03d}' - insertion = ins_form.format(int(datetime[3]), - int(datetime[4]), - int(datetime[5]), - int(datetime[6] * 1e-3)) - newline = line[:16] + insertion + line[28:] - return newline - - -def evt_head_check(root_dir, out_dir=None): - """ - A function to make sure that an arbitrary number of .gse files have correct values in their header. - :param root_dir: a directory leading to the .gse files. - :type root_dir: string - :param out_dir: a directory to store the new files somwhere els. - :return: returns nothing - """ - if not out_dir: - print('WARNING files are going to be overwritten!') - inp = str(input('Continue? [y/N]')) - if not inp == 'y': - sys.exit() - filelist = get_file_list(root_dir) - nfiles = 0 - for file in filelist: - infile = open(os.path.join(root_dir, file), 'r') - lines = infile.readlines() - infile.close() - datetime = time_from_header(lines[1]) - if check_time(datetime): - continue - else: - nfiles += 1 - datetime = checks_station_second(datetime, file) - print('writing ' + file) - # write File - lines[1] = make_time_line(lines[1], datetime) - if not out_dir: - out = open(os.path.join(root_dir, file), 'w') - out.writelines(lines) - out.close() - else: - out = open(os.path.join(out_dir, file), 'w') - out.writelines(lines) - out.close() - print(nfiles) - - -def read_metadata(path_to_inventory): - """ - take path_to_inventory and return either the corresponding list of files - found or the Parser object for a network dataless seed volume to prevent - read overhead for large dataless seed volumes - :param path_to_inventory: - :return: tuple containing a either list of files or `obspy.io.xseed.Parser` - object and the inventory type found - :rtype: tuple - """ - dlfile = list() - invfile = list() - respfile = list() - # possible file extensions specified here: - inv = dict(dless=dlfile, xml=invfile, resp=respfile, dseed=dlfile[:]) - if os.path.isfile(path_to_inventory): - ext = os.path.splitext(path_to_inventory)[1].split('.')[1] - inv[ext] += [path_to_inventory] - else: - for ext in inv.keys(): - inv[ext] += glob.glob1(path_to_inventory, '*.{0}'.format(ext)) - - invtype = key_for_set_value(inv) - - if invtype is None: - print("Neither dataless-SEED file, inventory-xml file nor " - "RESP-file found!") - print("!!WRONG CALCULATION OF SOURCE PARAMETERS!!") - robj = None, - elif invtype == 'dless': # prevent multiple read of large dlsv - print("Reading metadata information from dataless-SEED file ...") - if len(inv[invtype]) == 1: - fullpath_inv = os.path.join(path_to_inventory, inv[invtype][0]) - robj = Parser(fullpath_inv) - else: - robj = inv[invtype] - else: - print("Reading metadata information from inventory-xml file ...") - robj = read_inventory(inv[invtype]) - return invtype, robj - - -# idea to optimize read_metadata -# def read_metadata_new(path_to_inventory): -# metadata_objects = [] -# # read multiple files from directory -# if os.path.isdir(path_to_inventory): -# fnames = os.listdir(path_to_inventory) -# # read single file -# elif os.path.isfile(path_to_inventory): -# fnames = [path_to_inventory] -# else: -# print("Neither dataless-SEED file, inventory-xml file nor " -# "RESP-file found!") -# print("!!WRONG CALCULATION OF SOURCE PARAMETERS!!") -# fnames = [] -# -# for fname in fnames: -# path_to_inventory_filename = os.path.join(path_to_inventory, fname) -# try: -# ftype, robj = read_metadata_file(path_to_inventory_filename) -# metadata_objects.append((ftype, robj)) -# except Exception as e: -# print('Could not read metadata file {} ' -# 'because of the following Exception: {}'.format(path_to_inventory_filename, e)) -# return metadata_objects - - def restitute_trace(input_tuple): def no_metadata(tr, seed_id): print('no metadata file found ' diff --git a/pylot/core/util/generate_array_maps.py b/pylot/core/util/generate_array_maps.py index 0e7a9d64..69d3bb59 100755 --- a/pylot/core/util/generate_array_maps.py +++ b/pylot/core/util/generate_array_maps.py @@ -20,7 +20,7 @@ import matplotlib matplotlib.use('Qt5Agg') sys.path.append(os.path.join('/'.join(sys.argv[0].split('/')[:-1]), '../../..')) -from PyLoT import Project +from pylot.core.io.project import Project from pylot.core.util.dataprocessing import Metadata from pylot.core.util.array_map import Array_map diff --git a/pylot/core/util/structure.py b/pylot/core/util/structure.py index dd6cffd9..571f5dab 100644 --- a/pylot/core/util/structure.py +++ b/pylot/core/util/structure.py @@ -6,7 +6,7 @@ Created on Wed Jan 26 17:47:25 2015 @author: sebastianw """ -from pylot.core.io.data import SeiscompDataStructure, PilotDataStructure, ObspyDMTdataStructure +from pylot.core.io.data import SeiscompDataStructure, PilotDataStructure DATASTRUCTURE = {'PILOT': PilotDataStructure, 'SeisComP': SeiscompDataStructure, - 'obspyDMT': ObspyDMTdataStructure, None: PilotDataStructure} + 'obspyDMT': PilotDataStructure, None: PilotDataStructure} diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 0557db83..0c74ddfe 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -20,6 +20,10 @@ from pylot.core.io.inputs import PylotParameter, FilterOptions from pylot.core.util.obspyDMT_interface import check_obspydmt_eventfolder from pylot.styles import style_settings +Rgba: Type[tuple] = Tuple[int, int, int, int] +Mplrgba: Type[tuple] = Tuple[float, float, float, float] +Mplrgbastr: Type[tuple] = Tuple[str, str, str, str] + def _pickle_method(m): if m.im_self is None: @@ -83,25 +87,6 @@ def fit_curve(x, y): return splev, splrep(x, y) -def getindexbounds(f, eta): - """ - Get indices of values closest below and above maximum value in an array - :param f: array - :type f: `~numpy.ndarray` - :param eta: look for value in array that is closes to max_value * eta - :type eta: float - :return: tuple containing index of max value, index of value closest below max value, - index of value closest above max value - :rtype: (int, int, int) - """ - mi = f.argmax() # get indices of max values - m = max(f) # get maximum value - b = m * eta # - l = find_nearest(f[:mi], b) # find closest value below max value - u = find_nearest(f[mi:], b) + mi # find closest value above max value - return mi, l, u - - def gen_Pool(ncores=0): """ Generate mulitprocessing pool object utilizing ncores amount of cores @@ -121,7 +106,7 @@ def gen_Pool(ncores=0): print('gen_Pool: Generated multiprocessing Pool with {} cores\n'.format(ncores)) - pool = multiprocessing.Pool(ncores, maxtasksperchild=100) + pool = multiprocessing.Pool(ncores) return pool @@ -167,11 +152,11 @@ def clims(lim1, lim2): """ takes two pairs of limits and returns one pair of common limts :param lim1: limit 1 - :type lim1: int + :type lim1: List[int] :param lim2: limit 2 - :type lim2: int + :type lim2: List[int] :return: new upper and lower limit common to both given limits - :rtype: [int, int] + :rtype: List[int] >>> clims([0, 4], [1, 3]) [0, 4] @@ -303,7 +288,7 @@ def fnConstructor(s): if type(s) is str: s = s.split(':')[-1] else: - s = getHash(UTCDateTime()) + s = get_hash(UTCDateTime()) badchars = re.compile(r'[^A-Za-z0-9_. ]+|^\.|\.$|^ | $|^$') badsuffix = re.compile(r'(aux|com[1-9]|con|lpt[1-9]|prn)(\.|$)') @@ -315,15 +300,32 @@ def fnConstructor(s): return fn -def get_None(value): +def get_none(value): """ Convert "None" to None :param value: - :type value: str, bool + :type value: str, NoneType :return: - :rtype: bool + :rtype: type(value) or NoneType + + >>> st = read() + >>> print(get_none(st)) + 3 Trace(s) in Stream: + BW.RJOB..EHZ | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples + BW.RJOB..EHN | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples + BW.RJOB..EHE | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples + >>> get_none('Stream') + 'Stream' + >>> get_none(0) + 0 + >>> get_none(0.) + 0.0 + >>> print(get_none('None')) + None + >>> print(get_none(None)) + None """ - if value == 'None': + if value is None or (type(value) is str and value == 'None'): return None else: return value @@ -331,20 +333,47 @@ def get_None(value): def get_bool(value): """ - Convert string representations of bools to their true boolean value + Convert string representations of bools to their true boolean value. Return value if it cannot be identified as bool. :param value: - :type value: str, bool + :type value: str, bool, int, float :return: true boolean value :rtype: bool + + >>> get_bool(True) + True + >>> get_bool(False) + False + >>> get_bool(0) + False + >>> get_bool(0.) + False + >>> get_bool(0.1) + True + >>> get_bool(2) + True + >>> get_bool(-1) + False + >>> get_bool(-0.3) + False + >>> get_bool(None) + None + >>> get_bool('Stream') + 'Stream' """ - if type(value) is bool: + if type(value) == bool: return value elif value in ['True', 'true']: return True elif value in ['False', 'false']: return False + elif isinstance(value, float) or isinstance(value, int): + if value > 0. or value > 0: + return True + else: + return False else: - return bool(value) + return value + def four_digits(year): """ @@ -355,8 +384,8 @@ def four_digits(year): :return: four digit year correspondent :rtype: int - >>> four_digits(20) - 1920 + >>> four_digits(75) + 1975 >>> four_digits(16) 2016 >>> four_digits(00) @@ -438,36 +467,53 @@ def backtransformFilterString(st): return st -def getHash(time): +def get_hash(time): """ takes a time object and returns the corresponding SHA1 hash of the formatted date string :param time: time object for which a hash should be calculated :type time: `~obspy.core.utcdatetime.UTCDateTime` :return: SHA1 hash :rtype: str + + >>> time = UTCDateTime(0) + >>> get_hash(time) + '7627cce3b1b58dd21b005dac008b34d18317dd15' + >>> get_hash(0) + Traceback (most recent call last): + ... + AssertionError: 'time' is not an ObsPy UTCDateTime object """ + assert isinstance(time, UTCDateTime), '\'time\' is not an ObsPy UTCDateTime object' hg = hashlib.sha1() - hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f')) + hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f').encode('utf-8')) return hg.hexdigest() -def getLogin(): +def get_login(): """ - returns the actual user's login ID - :return: login ID + returns the actual user's name + :return: login name :rtype: str """ import getpass return getpass.getuser() -def getOwner(fn): +def get_owner(fn): """ takes a filename and return the login ID of the actual owner of the file :param fn: filename of the file tested :type fn: str :return: login ID of the file's owner :rtype: str + + >>> import tempfile + >>> with tempfile.NamedTemporaryFile() as tmpfile: + ... tmpfile.write(b'') and True + ... tmpfile.flush() + ... get_owner(tmpfile.name) == os.path.expanduser('~').split('/')[-1] + 0 + True """ system_name = platform.system() if system_name in ["Linux", "Darwin"]: @@ -513,6 +559,11 @@ def is_executable(fn): :param fn: path to the file to be tested :return: True or False :rtype: bool + + >>> is_executable('/bin/ls') + True + >>> is_executable('/var/log/system.log') + False """ return os.path.isfile(fn) and os.access(fn, os.X_OK) @@ -539,24 +590,36 @@ def isSorted(iterable): >>> isSorted([2,3,1,4]) False """ - assert isIterable(iterable), 'object is not iterable; object: {' \ - '}'.format(iterable) + assert is_iterable(iterable), "object is not iterable; object: {}".format(iterable) if type(iterable) is str: iterable = [s for s in iterable] return sorted(iterable) == iterable -def isIterable(obj): +def is_iterable(obj): """ takes a python object and returns True is the object is iterable and False otherwise :param obj: a python object - :type obj: object + :type obj: obj :return: True of False :rtype: bool + + >>> is_iterable(1) + False + >>> is_iterable(True) + False + >>> is_iterable(0.) + False + >>> is_iterable((0,1,3,4)) + True + >>> is_iterable([1]) + True + >>> is_iterable('a') + True """ try: - iterator = iter(obj) + iter(obj) except TypeError as te: return False return True @@ -565,13 +628,19 @@ def isIterable(obj): def key_for_set_value(d): """ takes a dictionary and returns the first key for which's value the - boolean is True + boolean representation is True :param d: dictionary containing values :type d: dict :return: key to the first non-False value found; None if no value's boolean equals True - :rtype: + :rtype: bool or NoneType + + >>> key_for_set_value({'one': 0, 'two': 1}) + 'two' + >>> print(key_for_set_value({1: 0, 2: False})) + None """ + assert type(d) is dict, "Function only defined for inputs of type 'dict'." r = None for k, v in d.items(): if v: @@ -579,32 +648,53 @@ def key_for_set_value(d): return r -def prepTimeAxis(stime, trace, verbosity=0): +def prep_time_axis(offset, trace, verbosity=0): """ - takes a starttime and a trace object and returns a valid time axis for + takes an offset and a trace object and returns a valid time axis for plotting - :param stime: start time of the actual seismogram as UTCDateTime - :type stime: `~obspy.core.utcdatetime.UTCDateTime` + :param offset: offset of the actual seismogram on plotting axis + :type offset: float or int :param trace: seismic trace object :type trace: `~obspy.core.trace.Trace` :param verbosity: if != 0, debug output will be written to console :type verbosity: int :return: valid numpy array with time stamps for plotting :rtype: `~numpy.ndarray` + + >>> tr = read()[0] + >>> prep_time_axis(0., tr) + array([0.00000000e+00, 1.00033344e-02, 2.00066689e-02, ..., + 2.99799933e+01, 2.99899967e+01, 3.00000000e+01]) + >>> prep_time_axis(22.5, tr) + array([22.5 , 22.51000333, 22.52000667, ..., 52.47999333, + 52.48999667, 52.5 ]) + >>> prep_time_axis(tr.stats.starttime, tr) + Traceback (most recent call last): + ... + AssertionError: 'offset' is not of type 'float' or 'int'; type: + >>> tr.stats.npts -= 1 + >>> prep_time_axis(0, tr) + array([0.00000000e+00, 1.00033356e-02, 2.00066711e-02, ..., + 2.99699933e+01, 2.99799967e+01, 2.99900000e+01]) + >>> tr.stats.npts += 2 + >>> prep_time_axis(0, tr) + array([0.00000000e+00, 1.00033333e-02, 2.00066667e-02, ..., + 2.99899933e+01, 2.99999967e+01, 3.00100000e+01]) """ + assert isinstance(offset, (float, int)), "'offset' is not of type 'float' or 'int'; type: {}".format(type(offset)) nsamp = trace.stats.npts srate = trace.stats.sampling_rate tincr = trace.stats.delta - etime = stime + nsamp / srate - time_ax = np.linspace(stime, etime, nsamp) + etime = offset + nsamp / srate + time_ax = np.linspace(offset, etime, nsamp) if len(time_ax) < nsamp: if verbosity: print('elongate time axes by one datum') - time_ax = np.arange(stime, etime + tincr, tincr) + time_ax = np.arange(offset, etime + tincr, tincr) elif len(time_ax) > nsamp: if verbosity: print('shorten time axes by one datum') - time_ax = np.arange(stime, etime - tincr, tincr) + time_ax = np.arange(offset, etime - tincr, tincr) if len(time_ax) != nsamp: print('Station {0}, {1} samples of data \n ' '{2} length of time vector \n' @@ -620,13 +710,13 @@ def find_horizontals(data): :param data: waveform data :type data: `obspy.core.stream.Stream` :return: components list - :rtype: list + :rtype: List(str) ..example:: >>> st = read() >>> find_horizontals(st) - [u'N', u'E'] + ['N', 'E'] """ rval = [] for tr in data: @@ -637,7 +727,7 @@ def find_horizontals(data): return rval -def pick_color(picktype, phase, quality=0): +def pick_color(picktype: Literal['manual', 'automatic'], phase: Literal['P', 'S'], quality: int = 0) -> Rgba: """ Create pick color by modifying the base color by the quality. @@ -650,7 +740,7 @@ def pick_color(picktype, phase, quality=0): :param quality: quality of pick. Decides the new intensity of the modifier color :type quality: int :return: tuple containing modified rgba color values - :rtype: (int, int, int, int) + :rtype: Rgba """ min_quality = 3 bpc = base_phase_colors(picktype, phase) # returns dict like {'modifier': 'g', 'rgba': (0, 0, 255, 255)} @@ -706,17 +796,17 @@ def pick_linestyle_plt(picktype, key): return linestyles[picktype][key] -def modify_rgba(rgba, modifier, intensity): +def modify_rgba(rgba: Rgba, modifier: Literal['r', 'g', 'b'], intensity: float) -> Rgba: """ Modify rgba color by adding the given intensity to the modifier color :param rgba: tuple containing rgba values - :type rgba: (int, int, int, int) - :param modifier: which color should be modified, eg. 'r', 'g', 'b' - :type modifier: str + :type rgba: Rgba + :param modifier: which color should be modified; options: 'r', 'g', 'b' + :type modifier: Literal['r', 'g', 'b'] :param intensity: intensity to be added to selected color :type intensity: float :return: tuple containing rgba values - :rtype: (int, int, int, int) + :rtype: Rgba """ rgba = list(rgba) index = {'r': 0, @@ -750,18 +840,20 @@ def transform_colors_mpl_str(colors, no_alpha=False): Transforms rgba color values to a matplotlib string of color values with a range of [0, 1] :param colors: tuple of rgba color values ranging from [0, 255] :type colors: (float, float, float, float) - :param no_alpha: Wether to return a alpha value in the matplotlib color string + :param no_alpha: Whether to return an alpha value in the matplotlib color string :type no_alpha: bool :return: String containing r, g, b values and alpha value if no_alpha is False (default) :rtype: str + + >>> transform_colors_mpl_str((255., 255., 255., 255.), True) + '(1.0, 1.0, 1.0)' + >>> transform_colors_mpl_str((255., 255., 255., 255.)) + '(1.0, 1.0, 1.0, 1.0)' """ - colors = list(colors) - colors_mpl = tuple([color / 255. for color in colors]) if no_alpha: - colors_mpl = '({}, {}, {})'.format(*colors_mpl) + return '({}, {}, {})'.format(*transform_colors_mpl(colors)) else: - colors_mpl = '({}, {}, {}, {})'.format(*colors_mpl) - return colors_mpl + return '({}, {}, {}, {})'.format(*transform_colors_mpl(colors)) def transform_colors_mpl(colors): @@ -771,27 +863,16 @@ def transform_colors_mpl(colors): :type colors: (float, float, float, float) :return: tuple of rgba color values ranging from [0, 1] :rtype: (float, float, float, float) + + >>> transform_colors_mpl((127.5, 0., 63.75, 255.)) + (0.5, 0.0, 0.25, 1.0) + >>> transform_colors_mpl(()) """ colors = list(colors) colors_mpl = tuple([color / 255. for color in colors]) return colors_mpl -def remove_underscores(data): - """ - takes a `obspy.core.stream.Stream` object and removes all underscores - from station names - :param data: stream of seismic data - :type data: `~obspy.core.stream.Stream` - :return: data stream - :rtype: `~obspy.core.stream.Stream` - """ - # for tr in data: - # # remove underscores - # tr.stats.station = tr.stats.station.strip('_') - return data - - def trim_station_components(data, trim_start=True, trim_end=True): """ cut a stream so only the part common to all three traces is kept to avoid dealing with offsets @@ -820,6 +901,19 @@ def trim_station_components(data, trim_start=True, trim_end=True): return data +def merge_stream(stream): + gaps = stream.get_gaps() + if gaps: + # list of merged stations (seed_ids) + merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps] + stream.merge(method=1) + print('Merged the following stations because of gaps:') + for merged_station in merged: + print(merged_station) + + return stream, gaps + + def check4gapsAndRemove(data): """ check for gaps in Stream and remove them @@ -840,12 +934,12 @@ def check4gapsAndRemove(data): return data -def check_for_gaps_and_merge(data): +def check4gapsAndMerge(data): """ check for gaps in Stream and merge if gaps are found :param data: stream of seismic data :type data: `~obspy.core.stream.Stream` - :return: data stream, gaps returned from obspy get_gaps + :return: data stream :rtype: `~obspy.core.stream.Stream` """ gaps = data.get_gaps() @@ -856,7 +950,7 @@ def check_for_gaps_and_merge(data): for merged_station in merged: print(merged_station) - return data, gaps + return data def check4doubled(data): @@ -928,11 +1022,11 @@ def get_possible_pylot_eventfile_extensions(event, fext): def get_stations(data): """ - Get list of all station names in data stream + Get list of all station names in data-stream :param data: stream containing seismic traces :type data: `~obspy.core.stream.Stream` :return: list of all station names in data, no duplicates - :rtype: list of str + :rtype: List(str) """ stations = [] for tr in data: @@ -959,66 +1053,87 @@ def check4rotated(data, metadata=None, verbosity=1): :rtype: `~obspy.core.stream.Stream` """ - def rotate_components(wfstream, metadata=None): + def rotation_required(trace_ids): + """ + Derive if any rotation is required from the orientation code of the input. + + :param trace_ids: string identifier of waveform data trace + :type trace_ids: List(str) + :return: boolean representing if rotation is necessary for any of the traces + :rtype: bool + """ + orientations = [trace_id[-1] for trace_id in trace_ids] + return any([orientation.isnumeric() for orientation in orientations]) + + def rotate_components(wfs_in, metadata=None): """ Rotate components if orientation code is numeric (= non traditional orientation). Azimut and dip are fetched from metadata. To be rotated, traces of a station have to be cut to the same length. Returns unrotated traces of no metadata is provided - :param wfstream: stream containing seismic traces of a station - :type wfstream: `~obspy.core.stream.Stream` + :param wfs_in: stream containing seismic traces of a station + :type wfs_in: `~obspy.core.stream.Stream` :param metadata: tuple containing metadata type string and metadata parser object :type metadata: (str, `~obspy.io.xseed.parser.Parser`) :return: stream object with traditionally oriented traces (ZNE) :rtype: `~obspy.core.stream.Stream` """ - # check if any traces in this station need to be rotated - trace_ids = [trace.id for trace in wfstream] - orientations = [trace_id[-1] for trace_id in trace_ids] - rotation_required = [orientation.isnumeric() for orientation in orientations] - if any(rotation_required): - t_start = full_range(wfstream) - try: - azimuts = [] - dips = [] - for tr_id in trace_ids: - azimuts.append(metadata.get_coordinates(tr_id, t_start)['azimuth']) - dips.append(metadata.get_coordinates(tr_id, t_start)['dip']) - except (KeyError, TypeError) as e: - print('Failed to rotate trace {}, no azimuth or dip available in metadata'.format(tr_id)) - return wfstream - if len(wfstream) < 3: - print('Failed to rotate Stream {}, not enough components available.'.format(wfstream)) - return wfstream - # to rotate all traces must have same length, so trim them - wfstream = trim_station_components(wfstream, trim_start=True, trim_end=True) - try: - 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 trace {} to ZNE'.format(trace_ids)) - # replace old data with rotated data, change the channel code to ZNE - z_index = dips.index(min( - dips)) # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90 being vertical) - 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: - coordinates = metadata.get_coordinates(trace_id, t_start) - dip, az = coordinates['dip'], coordinates['azimuth'] - trace = wfstream.select(id=trace_id)[0] - if az > 315 or 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' - except (ValueError) as e: - print(e) - return wfstream + if len(wfs_in) < 3: + print(f"Stream {wfs_in=}, has not enough components to rotate.") + return wfs_in - return wfstream + # check if any traces in this station need to be rotated + trace_ids = [trace.id for trace in wfs_in] + if not rotation_required(trace_ids): + logging.debug(f"Stream does not need any rotation: Traces are {trace_ids=}") + return wfs_in + + # check metadata quality + t_start = full_range(wfs_in) + try: + azimuths = [] + dips = [] + for tr_id in trace_ids: + azimuths.append(metadata.get_coordinates(tr_id, t_start)['azimuth']) + dips.append(metadata.get_coordinates(tr_id, t_start)['dip']) + except (KeyError, TypeError) as err: + logging.error(f"{type(err)=} occurred: {err=} Rotating not possible, not all azimuth and dip information " + f"available in metadata. Stream remains unchanged.") + return wfs_in + except Exception as err: + print(f"Unexpected {err=}, {type(err)=}") + raise + + # to rotate all traces must have same length, so trim them + wfs_out = trim_station_components(wfs_in, trim_start=True, trim_end=True) + try: + z, n, e = rotate2zne(wfs_out[0], azimuths[0], dips[0], + wfs_out[1], azimuths[1], dips[1], + wfs_out[2], azimuths[2], dips[2]) + print('check4rotated: rotated trace {} to ZNE'.format(trace_ids)) + # replace old data with rotated data, change the channel code to ZNE + z_index = dips.index(min( + dips)) # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90 + # being vertical) + wfs_out[z_index].data = z + wfs_out[z_index].stats.channel = wfs_out[z_index].stats.channel[0:-1] + 'Z' + del trace_ids[z_index] + for trace_id in trace_ids: + coordinates = metadata.get_coordinates(trace_id, t_start) + dip, az = coordinates['dip'], coordinates['azimuth'] + trace = wfs_out.select(id=trace_id)[0] + if az > 315 or az <= 45 or 135 < az <= 225: + trace.data = n + trace.stats.channel = trace.stats.channel[0:-1] + 'N' + elif 45 < az <= 135 or 225 < az <= 315: + trace.data = e + trace.stats.channel = trace.stats.channel[0:-1] + 'E' + except ValueError as err: + print(f"{err=} Rotation failed. Stream remains unchanged.") + return wfs_in + + return wfs_out if metadata is None: if verbosity: @@ -1032,38 +1147,6 @@ def check4rotated(data, metadata=None, verbosity=1): return data -def scaleWFData(data, factor=None, components='all'): - """ - produce scaled waveforms from given waveform data and a scaling factor, - waveform may be selected by their components name - :param data: waveform data to be scaled - :type data: `~obspy.core.stream.Stream` object - :param factor: scaling factor - :type factor: float - :param components: components labels for the traces in data to be scaled by - the scaling factor (optional, default: 'all') - :type components: tuple - :return: scaled waveform data - :rtype: `~obspy.core.stream.Stream` object - """ - if components != 'all': - for comp in components: - if factor is None: - max_val = np.max(np.abs(data.select(component=comp)[0].data)) - data.select(component=comp)[0].data /= 2 * max_val - else: - data.select(component=comp)[0].data /= 2 * factor - else: - for tr in data: - if factor is None: - max_val = float(np.max(np.abs(tr.data))) - tr.data /= 2 * max_val - else: - tr.data /= 2 * factor - - return data - - def runProgram(cmd, parameter=None): """ run an external program specified by cmd with parameters input returning the @@ -1135,7 +1218,6 @@ def identifyPhase(phase): return False -@lru_cache def identifyPhaseID(phase): """ Returns phase id (capital P or S) diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index a9a3e6f4..7ef01394 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -7,6 +7,7 @@ Created on Wed Mar 19 11:27:35 2014 import copy import datetime import getpass +import glob import multiprocessing import os import subprocess @@ -16,6 +17,7 @@ import traceback import matplotlib import numpy as np +from pylot.core.io.phases import getQualitiesfromxml matplotlib.use('QT5Agg') @@ -49,11 +51,11 @@ from pylot.core.pick.utils import getSNR, earllatepicker, getnoisewin, \ from pylot.core.pick.compare import Comparison from pylot.core.pick.autopick import fmpicker from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS -from pylot.core.util.utils import prepTimeAxis, full_range, demeanTrace, isSorted, findComboBoxIndex, clims, \ +from pylot.core.util.utils import prep_time_axis, full_range, demeanTrace, isSorted, findComboBoxIndex, clims, \ pick_linestyle_plt, pick_color_plt, \ check4rotated, check4doubled, check_for_gaps_and_merge, check_for_nan, identifyPhase, \ loopIdentifyPhase, trim_station_components, transformFilteroptions2String, \ - identifyPhaseID, get_bool, get_None, pick_color, getAutoFilteroptions, SetChannelComponents, \ + identifyPhaseID, get_bool, get_none, pick_color, getAutoFilteroptions, SetChannelComponents, \ station_id_remove_channel, get_pylot_eventfile_with_extension, get_possible_pylot_eventfile_extensions from autoPyLoT import autoPyLoT from pylot.core.util.thread import Thread @@ -138,18 +140,6 @@ class LogWidget(QtWidgets.QWidget): self.stderr.append(60 * '#' + '\n\n') -def getDataType(parent): - type = QInputDialog().getItem(parent, "Select phases type", "Type:", - ["manual", "automatic"]) - - if type[0].startswith('auto'): - type = 'auto' - else: - type = type[0] - - return type - - def plot_pdf(_axes, x, y, annotation, bbox_props, xlabel=None, ylabel=None, title=None): # try method or data @@ -793,7 +783,7 @@ class WaveformWidgetPG(QtWidgets.QWidget): def connect_signals(self): self.qcombo_processed.activated.connect(self.parent().newWF) - self.syn_checkbox.clicked.connect(self.parent().newWF) + self.comp_checkbox.clicked.connect(self.parent().newWF) def init_labels(self): self.label_layout.addWidget(self.status_label) @@ -804,13 +794,13 @@ class WaveformWidgetPG(QtWidgets.QWidget): # use widgets as placeholder, so that child widgets keep position when others are hidden mid_layout = QHBoxLayout() right_layout = QHBoxLayout() - mid_layout.addWidget(self.syn_checkbox) + mid_layout.addWidget(self.comp_checkbox) right_layout.addWidget(self.qcombo_processed) mid_widget.setLayout(mid_layout) right_widget.setLayout(right_layout) self.label_layout.addWidget(mid_widget) self.label_layout.addWidget(right_widget) - self.syn_checkbox.setLayoutDirection(Qt.RightToLeft) + self.comp_checkbox.setLayoutDirection(Qt.RightToLeft) self.label_layout.setStretch(0, 4) self.label_layout.setStretch(1, 0) self.label_layout.setStretch(2, 0) @@ -825,7 +815,7 @@ class WaveformWidgetPG(QtWidgets.QWidget): label = QtWidgets.QLabel() self.perm_labels.append(label) self.qcombo_processed = QtWidgets.QComboBox() - self.syn_checkbox = QtWidgets.QCheckBox('synthetics') + self.comp_checkbox = QtWidgets.QCheckBox('Load comparison data') self.addQCboxItem('processed', 'green') self.addQCboxItem('raw', 'black') # self.perm_qcbox_right.setAlignment(2) @@ -834,9 +824,11 @@ class WaveformWidgetPG(QtWidgets.QWidget): def getPlotDict(self): return self.plotdict - def activateObspyDMToptions(self, activate): - self.syn_checkbox.setVisible(activate) - self.qcombo_processed.setVisible(activate) + def activateObspyDMToptions(self, activate: bool) -> None: + self.qcombo_processed.setEnabled(activate) + + def activateCompareOptions(self, activate: bool) -> None: + self.comp_checkbox.setEnabled(activate) def setPermText(self, number, text=None, color='black'): if not 0 <= number < len(self.perm_labels): @@ -936,10 +928,10 @@ class WaveformWidgetPG(QtWidgets.QWidget): msg = 'plotting %s channel of station %s' % (channel, station) print(msg) stime = trace.stats.starttime - self.wfstart - time_ax = prepTimeAxis(stime, trace) + time_ax = prep_time_axis(stime, trace) if st_syn: stime_syn = trace_syn.stats.starttime - self.wfstart - time_ax_syn = prepTimeAxis(stime_syn, trace_syn) + time_ax_syn = prep_time_axis(stime_syn, trace_syn) if method == 'fast': trace.data, time_ax = self.minMax(trace, time_ax) @@ -959,7 +951,7 @@ class WaveformWidgetPG(QtWidgets.QWidget): [time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else []) trace.data = np.array( [datum * gain + n for index, datum in enumerate(trace.data) if not index % nth_sample]) - trace_syn.data = np.array([datum + n for index, datum in enumerate(trace_syn.data) + trace_syn.data = np.array([datum + n + shift_syn for index, datum in enumerate(trace_syn.data) if not index % nth_sample] if st_syn else []) plots.append((times, trace.data, times_syn, trace_syn.data)) @@ -1005,15 +997,6 @@ class WaveformWidgetPG(QtWidgets.QWidget): time_ax = np.linspace(time_ax[0], time_ax[-1], num=len(data)) return data, time_ax - # def getAxes(self): - # return self.axes - - # def getXLims(self): - # return self.getAxes().get_xlim() - - # def getYLims(self): - # return self.getAxes().get_ylim() - def setXLims(self, lims): vb = self.plotWidget.getPlotItem().getViewBox() vb.setXRange(float(lims[0]), float(lims[1]), padding=0) @@ -1148,12 +1131,12 @@ class PylotCanvas(FigureCanvas): ax.set_xlim(self.cur_xlim) ax.set_ylim(self.cur_ylim) self.refreshPickDlgText() - ax.figure.canvas.draw() + ax.figure.canvas.draw_idle() def panRelease(self, gui_event): self.press = None self.press_rel = None - self.figure.canvas.draw() + self.figure.canvas.draw_idle() def panZoom(self, gui_event, threshold=2., factor=1.1): if not gui_event.x and not gui_event.y: @@ -1167,8 +1150,6 @@ class PylotCanvas(FigureCanvas): break if not ax_check: return - # self.updateCurrentLimits() #maybe put this down to else: - # calculate delta (relative values in axis) old_x, old_y = self.press_rel xdiff = gui_event.x - old_x @@ -1371,110 +1352,145 @@ class PylotCanvas(FigureCanvas): plot_positions[channel] = plot_pos return plot_positions - def plotWFData(self, wfdata, title=None, zoomx=None, zoomy=None, + def plotWFData(self, wfdata, wfdata_compare=None, title=None, zoomx=None, zoomy=None, noiselevel=None, scaleddata=False, mapping=True, component='*', nth_sample=1, iniPick=None, verbosity=0, plot_additional=False, additional_channel=None, scaleToChannel=None, snr=None): + ax = self.prepare_plot() + self.clearPlotDict() + + wfstart, wfend = self.get_wf_range(wfdata) + compclass = self.get_comp_class() + plot_streams = self.get_plot_streams(wfdata, wfdata_compare, component, compclass) + + st_main = plot_streams['wfdata']['data'] + if mapping: + plot_positions = self.calcPlotPositions(st_main, compclass) + + nslc = self.get_sorted_nslc(st_main) + nmax = self.plot_traces(ax, plot_streams, nslc, wfstart, mapping, plot_positions, + scaleToChannel, noiselevel, scaleddata, nth_sample, verbosity) + + if plot_additional and additional_channel: + self.plot_additional_trace(ax, wfdata, additional_channel, scaleToChannel, + scaleddata, nth_sample, wfstart) + + self.finalize_plot(ax, wfstart, wfend, nmax, zoomx, zoomy, iniPick, title, snr) + + def prepare_plot(self): ax = self.axes[0] ax.cla() + return ax - self.clearPlotDict() - wfstart, wfend = full_range(wfdata) - nmax = 0 + def get_wf_range(self, wfdata): + return full_range(wfdata) + def get_comp_class(self): settings = QSettings() - compclass = SetChannelComponents.from_qsettings(settings) + return SetChannelComponents.from_qsettings(settings) - if not component == '*': - alter_comp = compclass.getCompPosition(component) - # alter_comp = str(alter_comp[0]) - - st_select = wfdata.select(component=component) - st_select += wfdata.select(component=alter_comp) - else: - st_select = wfdata - - if mapping: - plot_positions = self.calcPlotPositions(st_select, compclass) - - # list containing tuples of network, station, channel and plot position (for sorting) - nslc = [] - for plot_pos, trace in enumerate(st_select): - if not trace.stats.channel[-1] in ['Z', 'N', 'E', '1', '2', '3']: - print('Warning: Unrecognized channel {}'.format(trace.stats.channel)) - continue - nslc.append(trace.get_id()) - nslc.sort() - nslc.reverse() + def get_plot_streams(self, wfdata, wfdata_compare, component, compclass): + def get_wf_dict(data=Stream(), linecolor='k', offset=0., **plot_kwargs): + return dict(data=data, linecolor=linecolor, offset=offset, plot_kwargs=plot_kwargs) linecolor = (0., 0., 0., 1.) if not self.style else self.style['linecolor']['rgba_mpl'] + plot_streams = { + 'wfdata': get_wf_dict(linecolor=linecolor, linewidth=0.7), + 'wfdata_comp': get_wf_dict(offset=0.1, linecolor='b', alpha=0.7, linewidth=0.5) + } + if component != '*': + alter_comp = compclass.getCompPosition(component) + plot_streams['wfdata']['data'] = wfdata.select(component=component) + wfdata.select(component=alter_comp) + if wfdata_compare: + plot_streams['wfdata_comp']['data'] = wfdata_compare.select( + component=component) + wfdata_compare.select(component=alter_comp) + else: + plot_streams['wfdata']['data'] = wfdata + if wfdata_compare: + plot_streams['wfdata_comp']['data'] = wfdata_compare + + return plot_streams + + def get_sorted_nslc(self, st_main): + nslc = [trace.get_id() for trace in st_main if trace.stats.channel[-1] in ['Z', 'N', 'E', '1', '2', '3']] + nslc.sort(reverse=True) + return nslc + + def plot_traces(self, ax, plot_streams, nslc, wfstart, mapping, plot_positions, scaleToChannel, noiselevel, + scaleddata, nth_sample, verbosity): + nmax = 0 for n, seed_id in enumerate(nslc): network, station, location, channel = seed_id.split('.') - st = st_select.select(id=seed_id) - trace = st[0].copy() - if mapping: - n = plot_positions[trace.stats.channel] - if n > nmax: - nmax = n - if verbosity: - msg = 'plotting %s channel of station %s' % (channel, station) - print(msg) - stime = trace.stats.starttime - wfstart - time_ax = prepTimeAxis(stime, trace) - if time_ax is not None: - if scaleToChannel: - st_scale = wfdata.select(channel=scaleToChannel) - if st_scale: - tr = st_scale[0] - trace.detrend('constant') - trace.normalize(np.max(np.abs(tr.data)) * 2) - scaleddata = True - if not scaleddata: - trace.detrend('constant') - trace.normalize(np.max(np.abs(trace.data)) * 2) - - times = [time for index, time in enumerate(time_ax) if not index % nth_sample] - data = [datum + n for index, datum in enumerate(trace.data) if not index % nth_sample] - ax.axhline(n, color="0.5", lw=0.5) - ax.plot(times, data, color=linecolor, linewidth=0.7) - if noiselevel is not None: - for level in [-noiselevel[channel], noiselevel[channel]]: - ax.plot([time_ax[0], time_ax[-1]], - [n + level, n + level], - color=linecolor, - linestyle='dashed') + for wf_name, wf_dict in plot_streams.items(): + st_select = wf_dict.get('data') + if not st_select: + continue + trace = st_select.select(id=seed_id)[0].copy() + if mapping: + n = plot_positions[trace.stats.channel] + if n > nmax: + nmax = n + if verbosity: + print(f'plotting {channel} channel of station {station}') + time_ax = prep_time_axis(trace.stats.starttime - wfstart, trace) + self.plot_trace(ax, trace, time_ax, wf_dict, n, scaleToChannel, noiselevel, scaleddata, nth_sample) self.setPlotDict(n, seed_id) - if plot_additional and additional_channel: - compare_stream = wfdata.select(channel=additional_channel) - if compare_stream: - trace = compare_stream[0] - if scaleToChannel: - st_scale = wfdata.select(channel=scaleToChannel) - if st_scale: - tr = st_scale[0] - trace.detrend('constant') - trace.normalize(np.max(np.abs(tr.data)) * 2) - scaleddata = True - if not scaleddata: - trace.detrend('constant') - trace.normalize(np.max(np.abs(trace.data)) * 2) - time_ax = prepTimeAxis(stime, trace) - times = [time for index, time in enumerate(time_ax) if not index % nth_sample] - p_data = compare_stream[0].data - # #normalize - # p_max = max(abs(p_data)) - # p_data /= p_max - for index in range(3): - ax.plot(times, p_data, color='red', alpha=0.5, linewidth=0.7) - p_data += 1 + return nmax + def plot_trace(self, ax, trace, time_ax, wf_dict, n, scaleToChannel, noiselevel, scaleddata, nth_sample): + if time_ax is not None: + if scaleToChannel: + self.scale_trace(trace, scaleToChannel) + scaleddata = True + if not scaleddata: + trace.detrend('constant') + trace.normalize(np.max(np.abs(trace.data)) * 2) + offset = wf_dict.get('offset') + + times = [time for index, time in enumerate(time_ax) if not index % nth_sample] + data = [datum + n + offset for index, datum in enumerate(trace.data) if not index % nth_sample] + ax.axhline(n, color="0.5", lw=0.5) + ax.plot(times, data, color=wf_dict.get('linecolor'), **wf_dict.get('plot_kwargs')) + if noiselevel is not None: + self.plot_noise_level(ax, time_ax, noiselevel, channel, n, wf_dict.get('linecolor')) + + def scale_trace(self, trace, scaleToChannel): + st_scale = wfdata.select(channel=scaleToChannel) + if st_scale: + tr = st_scale[0] + trace.detrend('constant') + trace.normalize(np.max(np.abs(tr.data)) * 2) + + def plot_noise_level(self, ax, time_ax, noiselevel, channel, n, linecolor): + for level in [-noiselevel[channel], noiselevel[channel]]: + ax.plot([time_ax[0], time_ax[-1]], [n + level, n + level], color=linecolor, linestyle='dashed') + + def plot_additional_trace(self, ax, wfdata, additional_channel, scaleToChannel, scaleddata, nth_sample, wfstart): + compare_stream = wfdata.select(channel=additional_channel) + if compare_stream: + trace = compare_stream[0] + if scaleToChannel: + self.scale_trace(trace, scaleToChannel) + scaleddata = True + if not scaleddata: + trace.detrend('constant') + trace.normalize(np.max(np.abs(trace.data)) * 2) + time_ax = prep_time_axis(trace.stats.starttime - wfstart, trace) + self.plot_additional_data(ax, trace, time_ax, nth_sample) + + def plot_additional_data(self, ax, trace, time_ax, nth_sample): + times = [time for index, time in enumerate(time_ax) if not index % nth_sample] + p_data = trace.data + for index in range(3): + ax.plot(times, p_data, color='red', alpha=0.5, linewidth=0.7) + p_data += 1 + + def finalize_plot(self, ax, wfstart, wfend, nmax, zoomx, zoomy, iniPick, title, snr): if iniPick: - ax.vlines(iniPick, ax.get_ylim()[0], ax.get_ylim()[1], - colors='m', linestyles='dashed', - linewidth=2) - xlabel = 'seconds since {0}'.format(wfstart) + ax.vlines(iniPick, ax.get_ylim()[0], ax.get_ylim()[1], colors='m', linestyles='dashed', linewidth=2) + xlabel = f'seconds since {wfstart}' ylabel = '' self.updateWidget(xlabel, ylabel, title) self.setXLims(ax, [0, wfend - wfstart]) @@ -1484,15 +1500,14 @@ class PylotCanvas(FigureCanvas): if zoomy is not None: self.setYLims(ax, zoomy) if snr is not None: - if snr < 2: - warning = 'LOW SNR' - if snr < 1.5: - warning = 'VERY LOW SNR' - ax.text(0.1, 0.9, 'WARNING - {}'.format(warning), ha='center', va='center', transform=ax.transAxes, - color='red') - + self.plot_snr_warning(ax, snr) self.draw() + def plot_snr_warning(self, ax, snr): + if snr < 2: + warning = 'LOW SNR' if snr >= 1.5 else 'VERY LOW SNR' + ax.text(0.1, 0.9, f'WARNING - {warning}', ha='center', va='center', transform=ax.transAxes, color='red') + @staticmethod def getXLims(ax): return ax.get_xlim() @@ -1574,6 +1589,8 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog): self.events = events self.filepaths = [] self.file_extensions = [] + self.check_all_state = True + self.merge_strategy = None self.default_text = default_text self.label = label self.setButtons() @@ -1581,16 +1598,17 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog): self.connectSignals() self.showPaths() self.refreshSelectionBox() - self.refresh_timer = QTimer(self) - self.refresh_timer.timeout.connect(self.showPaths) - self.refresh_timer.start(10000) + # self.refresh_timer = QTimer(self) + # self.refresh_timer.timeout.connect(self.showPaths) + # self.refresh_timer.start(10000) self.resize(800, 450) - def setupUi(self): + ncol = 4 self.main_layout = QtWidgets.QVBoxLayout() self.header_layout = QtWidgets.QHBoxLayout() + self.footer_layout = QtWidgets.QHBoxLayout() # self.setLayout(self.main_layout) @@ -1604,11 +1622,24 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog): self.searchButton = QtWidgets.QPushButton('Search') self.searchButton.setVisible(False) + # check/uncheck button for table + self.checkAllButton = QtWidgets.QPushButton('Check/Uncheck all') + + # radiobutton for merge selection + self.mergeRadioButtonGroup = QtWidgets.QButtonGroup() + self.merge_button = QtWidgets.QRadioButton('Merge') + self.overwrite_button = QtWidgets.QRadioButton('Overwrite') + self.mergeRadioButtonGroup.addButton(self.merge_button) + self.mergeRadioButtonGroup.addButton(self.overwrite_button) + self.merge_button.setChecked(True) + self.merge_strategy = self.merge_button.text() + + # table self.tableWidget = QtWidgets.QTableWidget() tableWidget = self.tableWidget - tableWidget.setColumnCount(3) + tableWidget.setColumnCount(ncol) tableWidget.setRowCount(len(self.events)) - tableWidget.setHorizontalHeaderLabels(('Event ID', 'Filename', 'Last modified')) + tableWidget.setHorizontalHeaderLabels(('', 'Event ID', 'Filename', 'Last modified')) tableWidget.setEditTriggers(tableWidget.NoEditTriggers) tableWidget.setSortingEnabled(True) header = tableWidget.horizontalHeader() @@ -1621,9 +1652,17 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog): self.header_layout.addWidget(self.comboBox) self.header_layout.addWidget(self.searchButton) + self.footer_layout.addWidget(self.checkAllButton) + self.footer_layout.addWidget(self.statusText) + self.footer_layout.addWidget(self.merge_button) + self.footer_layout.addWidget(self.overwrite_button) + + self.footer_layout.setStretch(0, 0) + self.footer_layout.setStretch(1, 1) + self.main_layout.addLayout(self.header_layout) self.main_layout.addWidget(self.tableWidget) - self.main_layout.addWidget(self.statusText) + self.main_layout.addLayout(self.footer_layout) self.main_layout.addWidget(self._buttonbox) def showPaths(self): @@ -1632,23 +1671,23 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog): self.tableWidget.clearContents() for index, event in enumerate(self.events): filename = get_pylot_eventfile_with_extension(event, fext) - self.tableWidget.setItem(index, 0, QtWidgets.QTableWidgetItem(f'{event.pylot_id}')) + pf_selected_item = QtWidgets.QTableWidgetItem() + check_state = QtCore.Qt.Checked if filename else QtCore.Qt.Unchecked + pf_selected_item.setCheckState(check_state) + self.tableWidget.setItem(index, 0, pf_selected_item) + self.tableWidget.setItem(index, 1, QtWidgets.QTableWidgetItem(f'{event.pylot_id}')) if filename: self.filepaths.append(filename) ts = int(os.path.getmtime(filename)) # create QTableWidgetItems of filepath and last modification time fname_item = QtWidgets.QTableWidgetItem(f'{os.path.split(filename)[-1]}') + fname_item.setData(3, filename) ts_item = QtWidgets.QTableWidgetItem(f'{datetime.datetime.fromtimestamp(ts)}') - self.tableWidget.setItem(index, 1, fname_item) - self.tableWidget.setItem(index, 2, ts_item) + self.tableWidget.setItem(index, 2, fname_item) + self.tableWidget.setItem(index, 3, ts_item) - # TODO: Idea -> only refresh if table contents changed. Use selection to load only a subset of files - if len(self.filepaths) > 0: - status_text = f'Found {len(self.filepaths)} eventfiles. Do you want to load them?' - else: - status_text = 'Did not find any files for specified file mask.' - self.statusText.setText(status_text) + self.update_status() def refreshSelectionBox(self): fext = self.comboBox.currentText() @@ -1668,12 +1707,52 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog): self._buttonbox = QDialogButtonBox(QDialogButtonBox.Ok | QDialogButtonBox.Cancel) + def toggleCheckAll(self): + self.check_all_state = not self.check_all_state + self.checkAll(self.check_all_state) + + def checkAll(self, state): + state = QtCore.Qt.Checked if state else QtCore.Qt.Unchecked + for row_ind in range(self.tableWidget.rowCount()): + item = self.tableWidget.item(row_ind, 0) + item.setCheckState(state) + + def getChecked(self): + filepaths = [] + for row_ind in range(self.tableWidget.rowCount()): + item_check = self.tableWidget.item(row_ind, 0) + if item_check.checkState() == QtCore.Qt.Checked: + item_fname = self.tableWidget.item(row_ind, 2) + if item_fname: + filepath = item_fname.data(3) + filepaths.append(filepath) + return filepaths + + def update_status(self, row=None, col=None): + if col is not None and col != 0: + return + filepaths = self.getChecked() + if len(filepaths) > 0: + status_text = f"Found {len(filepaths)} eventfile{'s' if len(filepaths) > 1 else ''}. Do you want to load them?" + else: + status_text = 'Did not find/select any files for specified file mask.' + self.statusText.setText(status_text) + + def update_merge_strategy(self): + for button in (self.merge_button, self.overwrite_button): + if button.isChecked(): + self.merge_strategy = button.text() + def connectSignals(self): self._buttonbox.accepted.connect(self.accept) self._buttonbox.rejected.connect(self.reject) self.comboBox.editTextChanged.connect(self.showPaths) self.searchButton.clicked.connect(self.showPaths) - + self.checkAllButton.clicked.connect(self.toggleCheckAll) + self.checkAllButton.clicked.connect(self.update_status) + self.tableWidget.cellClicked.connect(self.update_status) + self.merge_button.clicked.connect(self.update_merge_strategy) + self.overwrite_button.clicked.connect(self.update_merge_strategy) class SingleTextLineDialog(QtWidgets.QDialog): @@ -1780,8 +1859,8 @@ class PhaseDefaults(QtWidgets.QDialog): class PickDlg(QDialog): update_picks = QtCore.Signal(dict) - def __init__(self, parent=None, data=None, station=None, network=None, location=None, picks=None, - autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None, + def __init__(self, parent=None, data=None, data_compare=None, station=None, network=None, location=None, picks=None, + autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None, show_comp_data=False, event=None, filteroptions=None, model=None, wftype=None): super(PickDlg, self).__init__(parent, Qt.Window) self.orig_parent = parent @@ -1790,6 +1869,7 @@ class PickDlg(QDialog): # initialize attributes self.parameter = parameter self._embedded = embedded + self.showCompData = show_comp_data self.station = station self.network = network self.location = location @@ -1828,11 +1908,32 @@ class PickDlg(QDialog): else: self.filteroptions = FILTERDEFAULTS self.pick_block = False + + # set attribute holding data + if data is None or not data: + try: + data = parent.get_data().get_wf_data().copy() + self.data = data.select(station=station) + except AttributeError as e: + errmsg = 'You either have to put in a data or an appropriate ' \ + 'parent (PyLoT MainWindow) object: {0}'.format(e) + raise Exception(errmsg) + else: + self.data = data + self.data_compare = data_compare + self.nextStation = QtWidgets.QCheckBox('Continue with next station ') # comparison channel - self.compareChannel = QtWidgets.QComboBox() - self.compareChannel.activated.connect(self.resetPlot) + self.referenceChannel = QtWidgets.QComboBox() + self.referenceChannel.activated.connect(self.resetPlot) + + # comparison channel + self.compareCB = QtWidgets.QCheckBox() + self.compareCB.setChecked(self.showCompData) + self.compareCB.clicked.connect(self.switchCompData) + self.compareCB.clicked.connect(self.resetPlot) + self.compareCB.setVisible(bool(self.data_compare)) # scale channel self.scaleChannel = QtWidgets.QComboBox() @@ -1845,19 +1946,7 @@ class PickDlg(QDialog): self.cur_xlim = None self.cur_ylim = None - # set attribute holding data - if data is None or not data: - try: - data = parent.get_data().getWFData().copy() - self.data = data.select(station=station) - except AttributeError as e: - errmsg = 'You either have to put in a data or an appropriate ' \ - 'parent (PyLoT MainWindow) object: {0}'.format(e) - raise Exception(errmsg) - else: - self.data = data - - self.stime, self.etime = full_range(self.getWFData()) + self.stime, self.etime = full_range(self.get_wf_data()) # initialize plotting widget self.multicompfig = PylotCanvas(parent=self, multicursor=True) @@ -1868,12 +1957,12 @@ class PickDlg(QDialog): self.setupUi() # fill compare and scale channels - self.compareChannel.addItem('-', None) + self.referenceChannel.addItem('-', None) self.scaleChannel.addItem('individual', None) - for trace in self.getWFData(): + for trace in self.get_wf_data(): channel = trace.stats.channel - self.compareChannel.addItem(channel, trace) + self.referenceChannel.addItem(channel, trace) if not channel[-1] in ['Z', 'N', 'E', '1', '2', '3']: print('Skipping unknown channel for scaling: {}'.format(channel)) continue @@ -1890,7 +1979,7 @@ class PickDlg(QDialog): if self.wftype is not None: title += ' | ({})'.format(self.wftype) - self.multicompfig.plotWFData(wfdata=self.getWFData(), + self.multicompfig.plotWFData(wfdata=self.get_wf_data(), wfdata_compare=self.get_wf_dataComp(), title=title) self.multicompfig.setZoomBorders2content() @@ -2066,8 +2155,11 @@ class PickDlg(QDialog): _dialtoolbar.addWidget(est_label) _dialtoolbar.addWidget(self.plot_arrivals_button) _dialtoolbar.addSeparator() - _dialtoolbar.addWidget(QtWidgets.QLabel('Compare to channel: ')) - _dialtoolbar.addWidget(self.compareChannel) + _dialtoolbar.addWidget(QtWidgets.QLabel('Plot reference channel: ')) + _dialtoolbar.addWidget(self.referenceChannel) + _dialtoolbar.addSeparator() + _dialtoolbar.addWidget(QtWidgets.QLabel('Compare: ')) + _dialtoolbar.addWidget(self.compareCB) _dialtoolbar.addSeparator() _dialtoolbar.addWidget(QtWidgets.QLabel('Scaling: ')) _dialtoolbar.addWidget(self.scaleChannel) @@ -2152,10 +2244,12 @@ class PickDlg(QDialog): station_id = trace.get_id() starttime = trace.stats.starttime station_coords = self.metadata.get_coordinates(station_id, starttime) + if not station_coords: + print('get_arrivals: No station coordinates found. Return!') + return origins = self.pylot_event.origins if phases == ['None', 'None']: - print("get_arrivals: Creation info (manual or auto) not available!") - print("Return!") + print("get_arrivals: Creation info (manual or auto) not available! Return!") return if origins: source_origin = origins[0] @@ -2267,7 +2361,7 @@ class PickDlg(QDialog): # create action and add to menu # phase name transferred using lambda function - slot = lambda phase=phase, phaseID=phaseID: phaseSelect[phaseID](phase) + slot = lambda ph=phase, phID=phaseID: phaseSelect[phID](ph) picksAction = createAction(parent=self, text=phase, slot=slot, shortcut=shortcut) @@ -2402,7 +2496,7 @@ class PickDlg(QDialog): def activatePicking(self): self.leave_rename_phase() self.renamePhaseAction.setEnabled(False) - self.compareChannel.setEnabled(False) + self.referenceChannel.setEnabled(False) self.scaleChannel.setEnabled(False) phase = self.currentPhase phaseID = self.getPhaseID(phase) @@ -2420,7 +2514,7 @@ class PickDlg(QDialog): if self.autoFilterAction.isChecked(): for filteraction in [self.filterActionP, self.filterActionS]: filteraction.setChecked(False) - self.filterWFData() + self.filter_wf_data() self.draw() else: self.draw() @@ -2434,7 +2528,7 @@ class PickDlg(QDialog): self.disconnectPressEvent() self.multicompfig.connectEvents() self.renamePhaseAction.setEnabled(True) - self.compareChannel.setEnabled(True) + self.referenceChannel.setEnabled(True) self.scaleChannel.setEnabled(True) self.connect_pick_delete() self.draw() @@ -2504,9 +2598,15 @@ class PickDlg(QDialog): def getGlobalLimits(self, ax, axis): return self.multicompfig.getGlobalLimits(ax, axis) - def getWFData(self): + def get_wf_data(self): return self.data + def get_wf_dataComp(self): + if self.showCompData: + return self.data_compare + else: + return Stream() + def selectWFData(self, channel): component = channel[-1].upper() wfdata = Stream() @@ -2516,17 +2616,17 @@ class PickDlg(QDialog): return tr if component == 'E' or component == 'N': - for trace in self.getWFData(): + for trace in self.get_wf_data(): trace = selectTrace(trace, 'NE') if trace: wfdata.append(trace) elif component == '1' or component == '2': - for trace in self.getWFData(): + for trace in self.get_wf_data(): trace = selectTrace(trace, '12') if trace: wfdata.append(trace) else: - wfdata = self.getWFData().select(component=component) + wfdata = self.get_wf_data().select(component=component) return wfdata def getPicks(self, picktype='manual'): @@ -2628,11 +2728,16 @@ class PickDlg(QDialog): stime = self.getStartTime() - # copy data for plotting - data = self.getWFData().copy() - data = self.getPickPhases(data, phase) - data.normalize() - if not data: + # copy wfdata for plotting + wfdata = self.get_wf_data().copy() + wfdata_comp = self.get_wf_dataComp().copy() + wfdata = self.getPickPhases(wfdata, phase) + wfdata_comp = self.getPickPhases(wfdata_comp, phase) + for wfd in [wfdata, wfdata_comp]: + if wfd: + wfd.normalize() + + if not wfdata: QtWidgets.QMessageBox.warning(self, 'No channel to plot', 'No channel to plot for phase: {}. ' 'Make sure to select the correct channels for P and S ' @@ -2640,14 +2745,16 @@ class PickDlg(QDialog): self.leave_picking_mode() return - # filter data and trace on which is picked prior to determination of SNR + # filter wfdata and trace on which is picked prior to determination of SNR filterphase = self.currentFilterPhase() if filterphase: filteroptions = self.getFilterOptions(filterphase).parseFilterOptions() try: - data.detrend('linear') - data.filter(**filteroptions) - # wfdata.filter(**filteroptions)# MP MP removed filtering of original data + for wfd in [wfdata, wfdata_comp]: + if wfd: + wfd.detrend('linear') + wfd.filter(**filteroptions) + # wfdata.filter(**filteroptions)# MP MP removed filtering of original wfdata except ValueError as e: self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Information, 'Denied', @@ -2657,8 +2764,8 @@ class PickDlg(QDialog): snr = [] noiselevels = {} # determine SNR and noiselevel - for trace in data.traces: - st = data.select(channel=trace.stats.channel) + for trace in wfdata.traces: + st = wfdata.select(channel=trace.stats.channel) stime_diff = trace.stats.starttime - stime result = getSNR(st, (noise_win, gap_win, signal_win), ini_pick - stime_diff) snr.append(result[0]) @@ -2669,23 +2776,25 @@ class PickDlg(QDialog): noiselevel = nfac noiselevels[trace.stats.channel] = noiselevel - # prepare plotting of data - for trace in data: - t = prepTimeAxis(trace.stats.starttime - stime, trace) - inoise = getnoisewin(t, ini_pick, noise_win, gap_win) - trace = demeanTrace(trace, inoise) - # upscale trace data in a way that each trace is vertically zoomed to noiselevel*factor - channel = trace.stats.channel - noiselevel = noiselevels[channel] - noiseScaleFactor = self.calcNoiseScaleFactor(noiselevel, zoomfactor=5.) - trace.data *= noiseScaleFactor - noiselevels[channel] *= noiseScaleFactor + # prepare plotting of wfdata + for wfd in [wfdata, wfdata_comp]: + if wfd: + for trace in wfd: + t = prep_time_axis(trace.stats.starttime - stime, trace) + inoise = getnoisewin(t, ini_pick, noise_win, gap_win) + trace = demeanTrace(trace, inoise) + # upscale trace wfdata in a way that each trace is vertically zoomed to noiselevel*factor + channel = trace.stats.channel + noiselevel = noiselevels[channel] + noiseScaleFactor = self.calcNoiseScaleFactor(noiselevel, zoomfactor=5.) + trace.data *= noiseScaleFactor + noiselevels[channel] *= noiseScaleFactor mean_snr = np.mean(snr) x_res = getResolutionWindow(mean_snr, parameter.get('extent')) xlims = [ini_pick - x_res, ini_pick + x_res] - ylims = list(np.array([-.5, .5]) + [0, len(data) - 1]) + ylims = list(np.array([-.5, .5]) + [0, len(wfdata) - 1]) title = self.getStation() + ' picking mode' title += ' | SNR: {}'.format(mean_snr) @@ -2693,9 +2802,10 @@ class PickDlg(QDialog): filtops_str = transformFilteroptions2String(filteroptions) title += ' | Filteroptions: {}'.format(filtops_str) - plot_additional = bool(self.compareChannel.currentText()) - additional_channel = self.compareChannel.currentText() - self.multicompfig.plotWFData(wfdata=data, + plot_additional = bool(self.referenceChannel.currentText()) + additional_channel = self.referenceChannel.currentText() + self.multicompfig.plotWFData(wfdata=wfdata, + wfdata_compare=wfdata_comp, title=title, zoomx=xlims, zoomy=ylims, @@ -2729,7 +2839,7 @@ class PickDlg(QDialog): filteroptions = None # copy and filter data for earliest and latest possible picks - wfdata = self.getWFData().copy().select(channel=channel) + wfdata = self.get_wf_data().copy().select(channel=channel) if filteroptions: try: wfdata.detrend('linear') @@ -2776,7 +2886,7 @@ class PickDlg(QDialog): minFMSNR = parameter.get('minFMSNR') quality = get_quality_class(spe, parameter.get('timeerrorsP')) if quality <= minFMweight and snr >= minFMSNR: - FM = fmpicker(self.getWFData().select(channel=channel).copy(), wfdata.copy(), parameter.get('fmpickwin'), + FM = fmpicker(self.get_wf_data().select(channel=channel).copy(), wfdata.copy(), parameter.get('fmpickwin'), pick - stime_diff) # save pick times for actual phase @@ -3068,7 +3178,7 @@ class PickDlg(QDialog): def togglePickBlocker(self): return not self.pick_block - def filterWFData(self, phase=None): + def filter_wf_data(self, phase=None): if not phase: phase = self.currentPhase if self.getPhaseID(phase) == 'P': @@ -3086,7 +3196,8 @@ class PickDlg(QDialog): self.cur_xlim = self.multicompfig.axes[0].get_xlim() self.cur_ylim = self.multicompfig.axes[0].get_ylim() # self.multicompfig.updateCurrentLimits() - data = self.getWFData().copy() + wfdata = self.get_wf_data().copy() + wfdata_comp = self.get_wf_dataComp().copy() title = self.getStation() if filter: filtoptions = None @@ -3094,19 +3205,22 @@ class PickDlg(QDialog): filtoptions = self.getFilterOptions(self.getPhaseID(phase), gui_filter=True).parseFilterOptions() if filtoptions is not None: - data.detrend('linear') - data.taper(0.02, type='cosine') - data.filter(**filtoptions) + for wfd in [wfdata, wfdata_comp]: + if wfd: + wfd.detrend('linear') + wfd.taper(0.02, type='cosine') + wfd.filter(**filtoptions) filtops_str = transformFilteroptions2String(filtoptions) title += ' | Filteroptions: {}'.format(filtops_str) if self.wftype is not None: title += ' | ({})'.format(self.wftype) - plot_additional = bool(self.compareChannel.currentText()) - additional_channel = self.compareChannel.currentText() + plot_additional = bool(self.referenceChannel.currentText()) + additional_channel = self.referenceChannel.currentText() scale_channel = self.scaleChannel.currentText() - self.multicompfig.plotWFData(wfdata=data, title=title, + self.multicompfig.plotWFData(wfdata=wfdata, wfdata_compare=wfdata_comp, + title=title, zoomx=self.getXLims(), zoomy=self.getYLims(), plot_additional=plot_additional, @@ -3120,14 +3234,14 @@ class PickDlg(QDialog): def filterP(self): self.filterActionS.setChecked(False) if self.filterActionP.isChecked(): - self.filterWFData('P') + self.filter_wf_data('P') else: self.refreshPlot() def filterS(self): self.filterActionP.setChecked(False) if self.filterActionS.isChecked(): - self.filterWFData('S') + self.filter_wf_data('S') else: self.refreshPlot() @@ -3179,11 +3293,14 @@ class PickDlg(QDialog): self.resetZoom() self.refreshPlot() + def switchCompData(self): + self.showCompData = self.compareCB.isChecked() + def refreshPlot(self): if self.autoFilterAction.isChecked(): self.filterActionP.setChecked(False) self.filterActionS.setChecked(False) - # data = self.getWFData().copy() + # data = self.get_wf_data().copy() # title = self.getStation() filter = False phase = None @@ -3683,8 +3800,8 @@ class TuneAutopicker(QWidget): fnames = self.station_ids[self.get_current_station_id()] if not fnames == self.fnames: self.fnames = fnames - self.data.setWFData(fnames) - wfdat = self.data.getWFData() # all available streams + self.data.set_wf_data(fnames) + wfdat = self.data.get_wf_data() # all available streams # remove possible underscores in station names # wfdat = remove_underscores(wfdat) # rotate misaligned stations to ZNE @@ -3789,12 +3906,14 @@ class TuneAutopicker(QWidget): network = None location = None - wfdata = self.data.getWFData() + wfdata = self.data.get_wf_data() + wfdata_comp = self.data.get_wf_dataComp() metadata = self.parent().metadata event = self.get_current_event() filteroptions = self.parent().filteroptions wftype = self.wftype if self.obspy_dmt else '' self.pickDlg = PickDlg(self.parent(), data=wfdata.select(station=station).copy(), + data_comp=wfdata_comp.select(station=station).copy(), station=station, network=network, location=location, parameter=self.parameter, picks=self.get_current_event_picks(station), @@ -3843,7 +3962,7 @@ class TuneAutopicker(QWidget): for plotitem in self._manual_pick_plots: self.clear_plotitem(plotitem) self._manual_pick_plots = [] - st = self.data.getWFData() + st = self.data.get_wf_data() tr = st.select(station=self.get_current_station())[0] starttime = tr.stats.starttime # create two lists with figure names and subindices (for subplots) to get the correct axes @@ -3879,7 +3998,7 @@ class TuneAutopicker(QWidget): self.plot_manual_pick_to_ax(ax=ax, picks=picks, phase='S', starttime=starttime, quality=qualitySpick) for canvas in self.parent().canvas_dict.values(): - canvas.draw() + canvas.draw_idle() def plot_manual_pick_to_ax(self, ax, picks, phase, starttime, quality): mpp = picks[phase]['mpp'] - starttime @@ -4794,8 +4913,8 @@ class InputsTab(PropTab): self.tstopBox = QSpinBox() for spinbox in [self.tstartBox, self.tstopBox]: spinbox.setRange(-99999, 99999) - self.tstartBox.setValue(float(settings.value('tstart')) if get_None(settings.value('tstart')) else 0) - self.tstopBox.setValue(float(settings.value('tstop')) if get_None(settings.value('tstop')) else 1e6) + self.tstartBox.setValue(float(settings.value('tstart')) if get_none(settings.value('tstart')) else 0) + self.tstopBox.setValue(float(settings.value('tstop')) if get_none(settings.value('tstop')) else 1e6) self.cuttimesLayout.addWidget(self.tstartBox, 10) self.cuttimesLayout.addWidget(QLabel('[s] and'), 0) self.cuttimesLayout.addWidget(self.tstopBox, 10) @@ -5791,7 +5910,7 @@ class ChooseWaveFormWindow(QWidget): #self.currentSpectro = self.traces[ # self.chooseBoxTraces.currentText()[3:]][self.chooseBoxComponent.currentText()].spectrogram(show=False, title=t) #self.currentSpectro.show() - applyFFT() + self.applyFFT() def applyFFT(self, trace): tra = self.traces[self.chooseBoxTraces.currentText()[3:]]['Z'] diff --git a/requirements.txt b/requirements.txt index b5e3fe81..a80fadd8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,12 +1,9 @@ # This file may be used to create an environment using: # $ conda create --name --file # platform: win-64 -cartopy=0.20.2 -matplotlib-base=3.3.4 -numpy=1.22.3 -obspy=1.3.0 -pyqtgraph=0.12.4 -pyside2=5.13.2 -python=3.8.12 -qt=5.12.9 -scipy=1.8.0 \ No newline at end of file +cartopy>=0.20.2 +numpy<2 +obspy>=1.3.0 +pyqtgraph>=0.12.4 +pyside2>=5.13.2 +scipy>=1.8.0 \ No newline at end of file