diff --git a/.gitignore b/.gitignore index c99398e7..78b6b4e7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ *.pyc *~ pylot/RELEASE-VERSION +*.idea diff --git a/PyLoT.py b/PyLoT.py index 1c790609..60dee706 100755 --- a/PyLoT.py +++ b/PyLoT.py @@ -86,7 +86,7 @@ from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \ getDataType, ComparisonWidget, TuneAutopicker, PylotParaBox, AutoPickDlg, CanvasWidget, AutoPickWidget, \ CompareEventsWidget -from pylot.core.util.map_projection import map_projection +from pylot.core.util.array_map import Array_map from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.thread import Thread, Worker from pylot.core.util.version import get_git_version as _getVersionString @@ -138,6 +138,9 @@ class MainWindow(QMainWindow): self._eventChanged = [False, False] self.apd_local = None self.apd_sge = None + self.stations_highlighted = [] + self.obspy_dmt = False + self.nextStation = False self.poS_id = None self.ae_id = None @@ -615,6 +618,10 @@ class MainWindow(QMainWindow): # add scroll area used in case number of traces gets too high self.wf_scroll_area = QtGui.QScrollArea(self) + self.wf_scroll_area.setVisible(False) + self.no_data_label = QLabel('No Data') + self.no_data_label.setStyleSheet('color: red') + self.no_data_label.setAlignment(Qt.AlignCenter) # create central matplotlib figure canvas widget self.init_wfWidget() @@ -632,11 +639,17 @@ class MainWindow(QMainWindow): array_tab.setLayout(self.array_layout) events_tab.setLayout(self.events_layout) + # tighten up layouts inside tabs + for layout in [self.wf_layout, self.array_layout, self.events_layout]: + layout.setSpacing(0) + layout.setContentsMargins(0, 0, 0, 0) + # add tabs to main tab widget self.tabs.addTab(wf_tab, 'Waveform Plot') self.tabs.addTab(array_tab, 'Array Map') self.tabs.addTab(events_tab, 'Eventlist') + self.wf_layout.addWidget(self.no_data_label) self.wf_layout.addWidget(self.wf_scroll_area) self.wf_scroll_area.setWidgetResizable(True) self.init_array_tab() @@ -1217,6 +1230,8 @@ class MainWindow(QMainWindow): tv.verticalHeader().hide() tv.setSelectionBehavior(QtGui.QAbstractItemView.SelectRows) + current_event = self.get_current_event() + eventBox.setView(tv) eventBox.clear() model = eventBox.model() @@ -1249,6 +1264,17 @@ class MainWindow(QMainWindow): event_ref = event.isRefEvent() event_test = event.isTestEvent() + time = lat = lon = depth = mag = 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) == 1: + magnitude = event.magnitudes[0] + mag = magnitude.mag + # text = '{path:{plen}} | manual: [{p:3d}] | auto: [{a:3d}]' # text = text.format(path=event_path, # plen=plmax, @@ -1256,6 +1282,11 @@ class MainWindow(QMainWindow): # a=event_nautopicks) item_path = QtGui.QStandardItem('{path:{plen}}'.format(path=event_path, plen=plmax)) + item_time = QtGui.QStandardItem('{}'.format(time)) + item_lat = QtGui.QStandardItem('{}'.format(lat)) + item_lon = QtGui.QStandardItem('{}'.format(lon)) + item_depth = QtGui.QStandardItem('{}'.format(depth)) + item_mag = QtGui.QStandardItem('{}'.format(mag)) item_nmp = QtGui.QStandardItem(str(ma_count['manual'])) item_nmp.setIcon(self.manupicksicon_small) item_nap = QtGui.QStandardItem(str(ma_count['auto'])) @@ -1281,10 +1312,17 @@ class MainWindow(QMainWindow): # item.setFont(font) # item2.setForeground(QtGui.QColor('black')) # item2.setFont(font) - itemlist = [item_path, item_nmp, item_nap, item_ref, item_test, item_notes] - if event_test and select_events == 'ref': + itemlist = [item_path, item_time, item_lat, item_lon, item_depth, + item_mag, 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).strip(): message = ('Path missmatch creating eventbox.\n' @@ -1296,6 +1334,23 @@ class MainWindow(QMainWindow): eventBox.setCurrentIndex(index) self.refreshRefTestButtons() + def isEmpty(self, event_path): + wf_stat = {True: 'processed', + False: 'raw', + None: None} + + # self.data.processed is only None for PyLoT datastructure, else True or False + wf_dir = wf_stat[self.data.processed] + if wf_dir is not None: + wf_path = os.path.join(event_path, wf_dir) + if wf_dir is 'processed' and not os.path.exists(wf_path): + wf_path = os.path.join(event_path, 'raw') + else: + wf_path = event_path + if not os.path.exists(wf_path): + return True + return not bool(os.listdir(wf_path)) + def filename_from_action(self, action): if action.data() is None: filt = "Supported file formats" \ @@ -1373,7 +1428,10 @@ class MainWindow(QMainWindow): def exportAllEvents(self, outformats=['.xml']): for event in self.project.eventlist: self.get_data().setEvtData(event) - self.saveData(event, event.path, outformats) + try: + self.saveData(event, event.path, outformats) + except Exception as e: + print('WARNING! Could not save event {}. Reason: {}'.format(event.path, e)) def enableSaveEventAction(self): self.saveEventAction.setEnabled(True) @@ -1617,10 +1675,10 @@ class MainWindow(QMainWindow): if self._eventChanged[1]: self.refresh_array_map() if not plotted and self._eventChanged[0]: - # newWF(False) = load data without plotting + # newWF(plot=False) = load data without plotting self.newWF(plot=False) - def newWF(self, plot=True): + def newWF(self, event=None, plot=True): ''' Load new data and plot if necessary. ''' @@ -1634,6 +1692,7 @@ class MainWindow(QMainWindow): call modal plot thread method when finished. ''' if load: + self.prepareLoadWaveformData() self.wfd_thread = Thread(self, self.loadWaveformData, progressText='Reading data input...', pb_widget=self.mainProgressBarWidget) @@ -1646,6 +1705,16 @@ class MainWindow(QMainWindow): if plot and not load: self.plotWaveformDataThread() + def prepareLoadWaveformData(self): + self.fnames = self.getWFFnames_from_eventbox() + self.fnames_syn = [] + 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) + def loadWaveformData(self): ''' Load waveform data corresponding to current selected event. @@ -1657,14 +1726,48 @@ class MainWindow(QMainWindow): # ans = self.data.setWFData(self.getWFFnames()) # else: # ans = False - self.fnames = self.getWFFnames_from_eventbox() - eventpath = self.get_current_event_path() - basepath = eventpath.split(os.path.basename(eventpath))[0] - obspy_dmt = check_obspydmt_structure(basepath) + self.data.setWFData(self.fnames, + self.fnames_syn, checkRotated=True, - metadata=self.metadata, - obspy_dmt=obspy_dmt) + metadata=self.metadata) + + def prepareObspyDMT_data(self, eventpath): + qcbox_processed = self.dataPlot.qcombo_processed + qcheckb_syn = self.dataPlot.syn_checkbox + qcbox_processed.setEnabled(False) + qcheckb_syn.setEnabled(False) + for fpath in os.listdir(eventpath): + fpath = fpath.split('/')[-1] + 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 'processed' in fpath: + qcbox_processed.setEnabled(True) + if qcbox_processed.isEnabled(): + wftype = qcbox_processed.currentText() + else: + wftype = 'raw' + qcbox_processed.setCurrentIndex(qcbox_processed.findText(wftype)) + 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): """ @@ -1757,7 +1860,10 @@ class MainWindow(QMainWindow): def finish_pg_plot(self): self.getPlotWidget().updateWidget() - plots = self.wfp_thread.data + plots, gaps = self.wfp_thread.data + # do not show plot if no data are given + self.wf_scroll_area.setVisible(len(plots) > 0) + self.no_data_label.setVisible(not len(plots) > 0) for times, data, times_syn, data_syn in plots: self.dataPlot.plotWidget.getPlotItem().plot(times, data, pen=self.dataPlot.pen_linecolor) @@ -1765,8 +1871,7 @@ class MainWindow(QMainWindow): self.dataPlot.plotWidget.getPlotItem().plot(times_syn, data_syn, pen=self.dataPlot.pen_linecolor_syn) self.dataPlot.reinitMoveProxy() - self.dataPlot.plotWidget.showAxis('left') - self.dataPlot.plotWidget.showAxis('bottom') + self.highlight_stations() def finishWaveformDataPlot(self): self.comparable = self.checkEvents4comparison() @@ -1842,12 +1947,10 @@ class MainWindow(QMainWindow): comparable[event.pylot_id] = self.checkEvent4comparison(event) return comparable - def clearWaveformDataPlot(self): + def clearWaveformDataPlot(self, refresh_plot=False): self.disconnectWFplotEvents() if self.pg: self.dataPlot.plotWidget.getPlotItem().clear() - self.dataPlot.plotWidget.hideAxis('bottom') - self.dataPlot.plotWidget.hideAxis('left') else: for ax in self.dataPlot.axes: ax.cla() @@ -1863,6 +1966,9 @@ class MainWindow(QMainWindow): self.openEventAction.setEnabled(False) self.openEventsAutoAction.setEnabled(False) self.loadpilotevent.setEnabled(False) + if not refresh_plot: + self.wf_scroll_area.setVisible(False) + self.no_data_label.setVisible(True) self.disableSaveEventAction() self.draw() @@ -1871,7 +1977,7 @@ class MainWindow(QMainWindow): Open a modal thread to plot current waveform data. ''' self.check_plot_quantity() - self.clearWaveformDataPlot() + self.clearWaveformDataPlot(refresh_plot=True) self.wfp_thread = Thread(self, self.plotWaveformData, arg=filter, progressText='Plotting waveform data...', @@ -1904,9 +2010,10 @@ class MainWindow(QMainWindow): self.plot_method = 'fast' else: self.plot_method = 'normal' - plots = plotWidget.plotWFData(wfdata=wfst, wfsyn=wfsyn, title=title, mapping=False, component=comp, - nth_sample=int(nth_sample), method=self.plot_method) - return plots + rval = plotWidget.plotWFData(wfdata=wfst, wfsyn=wfsyn, title=title, mapping=False, component=comp, + nth_sample=int(nth_sample), method=self.plot_method) + plots, gaps = rval if rval else ([], []) + return plots, gaps def adjustPlotHeight(self): if self.pg: @@ -2135,10 +2242,12 @@ class MainWindow(QMainWindow): def pickOnStation(self, gui_event): if self.pg: - if not gui_event.button() == 1: + button = gui_event.button() + if not button in [1, 4]: return else: - if not gui_event.button == 1: + button = gui_event.button + if not button == 1: return if self.pg: @@ -2149,27 +2258,59 @@ class MainWindow(QMainWindow): wfID = self.getWFID(ycoord) if wfID is None: return - self.pickDialog(wfID) - def pickDialog(self, wfID, nextStation=False): - station = self.getStationName(wfID) network = self.getNetworkName(wfID) + station = self.getStationName(wfID) + if button == 1: + self.pickDialog(wfID, network, station) + elif button == 4: + self.toggle_station_color(wfID, network, station) + + def toggle_station_color(self, wfID, network, station): + black_pen = pg.mkPen((0, 0, 0)) + red_pen = pg.mkPen((200, 50, 50)) + line_item = self.dataPlot.plotWidget.getPlotItem().listDataItems()[wfID - 1] + current_pen = line_item.opts['pen'] + nwst = '{}.{}'.format(network, station) + if current_pen == black_pen: + line_item.setPen(red_pen) + if not nwst in self.stations_highlighted: + self.stations_highlighted.append(nwst) + else: + line_item.setPen(black_pen) + if nwst in self.stations_highlighted: + self.stations_highlighted.pop(self.stations_highlighted.index(nwst)) + + def highlight_stations(self): + for wfID, value in self.getPlotWidget().getPlotDict().items(): + station, channel, network = value + nwst = '{}.{}'.format(network, station) + if nwst in self.stations_highlighted: + self.toggle_station_color(wfID, network, station) + + def pickDialog(self, wfID, network=None, station=None): + if not network: + network = self.getNetworkName(wfID) if not station: + station = self.getStationName(wfID) + if not station or not network: return self.update_status('picking on station {0}'.format(station)) data = self.get_data().getOriginalWFData().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), station=station, network=network, picks=self.getPicksOnStation(station, 'manual'), autopicks=self.getPicksOnStation(station, 'auto'), metadata=self.metadata, event=event, - filteroptions=self.filteroptions) + filteroptions=self.filteroptions, wftype=wftype) if self.filterActionP.isChecked() or self.filterActionS.isChecked(): pickDlg.currentPhase = self.getSeismicPhase() pickDlg.filterWFData() - pickDlg.nextStation.setChecked(nextStation) + pickDlg.nextStation.setChecked(self.nextStation) + pickDlg.nextStation.stateChanged.connect(self.toggle_next_station) if pickDlg.exec_(): if pickDlg._dirty: self.setDirty(True) @@ -2179,13 +2320,13 @@ class MainWindow(QMainWindow): self.enableSaveEventAction() if replot1 or replot2: self.plotWaveformDataThread() - self.drawPicks() + self.draw() else: self.drawPicks(station) self.draw() - if pickDlg.nextStation.isChecked(): - self.pickDialog(wfID - 1, nextStation=pickDlg.nextStation.isChecked()) + if self.nextStation: + self.pickDialog(wfID - 1) else: self.update_status('picks discarded ({0})'.format(station)) if not self.get_loc_flag() and self.check4Loc(): @@ -2194,8 +2335,13 @@ class MainWindow(QMainWindow): elif self.get_loc_flag() and not self.check4Loc(): self.set_loc_flag(False) + def toggle_next_station(self, signal): + self.nextStation = bool(signal) + def addListItem(self, text): - self.listWidget.addItem(text) + textlist = text.split('\n') + for text in textlist: + self.listWidget.addItem(text) self.listWidget.scrollToBottom() def init_fig_dict(self): @@ -2254,7 +2400,7 @@ class MainWindow(QMainWindow): def tune_autopicker(self): ''' - Initiates TuneAutopicker widget use to interactively + Initiates TuneAutopicker widget used to interactively tune parameters for autopicking algorithm. ''' # figures and canvas have to be iniated from the main GUI @@ -2263,7 +2409,7 @@ class MainWindow(QMainWindow): self.init_fig_dict() #if not self.tap: # init TuneAutopicker object - self.tap = TuneAutopicker(self) + self.tap = TuneAutopicker(self, self.obspy_dmt) # first call of update to init tabs with empty canvas self.update_autopicker() # connect update signal of TuneAutopicker with update function @@ -2345,6 +2491,7 @@ class MainWindow(QMainWindow): # export current picks etc. self.exportAllEvents(['.xml']) + wfpath = self.dataPlot.qcombo_processed.currentText() if self.obspy_dmt else '' # define arguments for picker args = {'parameter': self._inputs, 'station': 'all', @@ -2353,8 +2500,8 @@ class MainWindow(QMainWindow): 'iplot': 0, 'fig_dict': None, 'fig_dict_wadatijack': self.fig_dict_wadatijack, - 'locflag': 0} - + 'savexml': False, + 'obspyDMT_wfpath': wfpath} # init pick thread self.mp_thread = QtCore.QThreadPool() self.mp_worker = Worker(autoPyLoT, args, redirect_stdout=True) @@ -2363,6 +2510,7 @@ class MainWindow(QMainWindow): self.mp_worker.signals.message.connect(self.addListItem) self.mp_worker.signals.result.connect(self.finalizeAutoPick) + self.mp_worker.signals.finished.connect(self.enableAutoPickWidget) self.mp_thread.start(self.mp_worker) @@ -2377,7 +2525,6 @@ class MainWindow(QMainWindow): self.apd_sge.show() def finalizeAutoPick(self, result): - self.apw.enable(True) if result: self.init_canvas_dict_wadatijack() for eventID in result.keys(): @@ -2393,6 +2540,9 @@ class MainWindow(QMainWindow): self.drawPicks(picktype='auto') self.draw() + def enableAutoPickWidget(self, event=None): + self.apw.enable(True) + def get_event_from_id(self, eventID): for event in self.project.eventlist: if event.pylot_id == eventID: @@ -2607,55 +2757,28 @@ class MainWindow(QMainWindow): """ if not self.okToContinue(): return - settings = QSettings() - # get location tool hook - loctool = settings.value("loc/tool", "nll") + parameter = self._inputs + loctool = 'nll' lt = locateTool[loctool] # get working directory - locroot = settings.value("{0}/rootPath".format(loctool), None) + locroot = parameter['nllocroot'] if locroot is None: self.PyLoTprefs() self.locate_event() - infile = settings.value("{0}/inputFile".format(loctool), None) + ctrfile = os.path.join(locroot, 'run', parameter['ctrfile']) - if not infile: - caption = 'Select {0} input file'.format(loctool) - filt = "Supported file formats" \ - " (*.in *.ini *.conf *.cfg)" - ans = QFileDialog().getOpenFileName(self, caption=caption, - filter=filt, dir=locroot) - if ans[0]: - infile = ans[0] - else: - QMessageBox.information(self, - self.tr('No infile selected'), - self.tr('Inputfile necessary for localization.')) - return - settings.setValue("{0}/inputFile".format(loctool), infile) - settings.sync() - if loctool == 'nll': - ttt = settings.value("{0}/travelTimeTables", None) - ok = False - if ttt is None: - while not ok: - text, ok = QInputDialog.getText(self, 'Pattern for travel time tables', - 'Base name of travel time tables', - echo=QLineEdit.Normal, - text="ttime") - ttt = text - - outfile = settings.value("{0}/outputFile".format(loctool), - os.path.split(os.tempnam())[-1]) + ttt = parameter['ttpatter'] + outfile = parameter['outpatter'] eventname = self.get_current_event_name() obsdir = os.path.join(self._inputs['rootpath'], self._inputs['datapath'], self._inputs['database'], eventname) self.saveData(event=self.get_current_event(), directory=obsdir, outformats='.obs') filename = 'PyLoT_' + eventname locpath = os.path.join(locroot, 'loc', filename) phasefile = os.path.join(obsdir, filename + '.obs') - lt.modify_inputs(infile, locroot, filename, phasefile, ttt) + lt.modify_inputs(ctrfile, locroot, filename, phasefile, ttt) try: - lt.locate(infile) + lt.locate(ctrfile) except RuntimeError as e: print(e.message) #finally: @@ -2721,10 +2844,7 @@ class MainWindow(QMainWindow): self.init_metadata() if not self.metadata: return - self.am_figure = Figure() - self.am_canvas = FigureCanvas(self.am_figure) - self.am_toolbar = NavigationToolbar(self.am_canvas, self) - self.array_map = map_projection(self) + self.array_map = Array_map(self) # self.array_map_thread() self.array_layout.addWidget(self.array_map) self.tabs.setCurrentIndex(index) @@ -2763,7 +2883,8 @@ class MainWindow(QMainWindow): lon = event.origins[0].longitude self.array_map.eventLoc = (lat, lon) if self.get_current_event(): - self.array_map.refresh_drawings(self.get_current_event().getPicks()) + self.array_map.refresh_drawings(self.get_current_event().getPicks(), + self.get_current_event().getAutopicks(),) self._eventChanged[1] = False def init_event_table(self, tabindex=2): @@ -2817,14 +2938,6 @@ class MainWindow(QMainWindow): event.addNotes(notes) self.fill_eventbox() - def set_background_color(items, color): - for item in items: - item.setBackground(color) - - def set_foreground_color(items, color): - for item in items: - item.setForeground(color) - current_event = self.get_current_event() # generate delete icon @@ -2894,7 +3007,7 @@ class MainWindow(QMainWindow): if hasattr(event, 'origins'): if event.origins: origin = event.origins[0] - item_time.setText(str(origin.time).split('.')[0]) + item_time.setText(str(origin.time + 0).split('.')[0]) # +0 as workaround in case time=0s item_lon.setText(str(origin.longitude)) item_lat.setText(str(origin.latitude)) item_depth.setText(str(origin.depth)) @@ -2928,11 +3041,7 @@ class MainWindow(QMainWindow): item_nmp, item_nap, item_ref, item_test, item_notes] self.project._table.append(column) - if index%2: - set_background_color(column, QtGui.QColor(*(245, 245, 245, 255))) - - if event == current_event: - set_foreground_color(column, QtGui.QColor(*(0, 143, 143, 255))) + self.setItemColor(column, index, event, current_event) # manipulate items item_ref.setBackground(self._ref_test_colors['ref']) @@ -2954,6 +3063,24 @@ class MainWindow(QMainWindow): self.events_layout.addWidget(self.event_table) self.tabs.setCurrentIndex(tabindex) + def setItemColor(self, item_list, index, event, current_event): + def set_background_color(items, color): + for item in items: + item.setBackground(color) + + def set_foreground_color(items, color): + for item in items: + item.setForeground(color) + + if index % 2: + set_background_color(item_list, QtGui.QColor(*(245, 245, 245, 255))) + + if self.isEmpty(event.path): + set_foreground_color(item_list, QtGui.QColor(*(180, 180, 180, 255))) + + if event == current_event: + set_background_color(item_list, QtGui.QColor(*(0, 143, 143, 255))) + def read_metadata_thread(self, fninv): self.rm_thread = Thread(self, read_metadata, arg=fninv, progressText='Reading metadata...', pb_widget=self.mainProgressBarWidget) @@ -3003,20 +3130,18 @@ class MainWindow(QMainWindow): settings.setValue("inventoryFile", self.project.inv_path) fninv = settings.value("inventoryFile", None) - if fninv and ask_default: + if (fninv and ask_default) and not new: ans = QMessageBox.question(self, self.tr("Use default metadata..."), self.tr( "Do you want to use the default value for metadata?\n({})".format(fninv)), QMessageBox.Yes | QMessageBox.No, QMessageBox.Yes) - if ans == QMessageBox.No: - if not set_inv(settings): - return None - elif ans == QMessageBox.Yes: + if ans == QMessageBox.Yes: self.read_metadata_thread(fninv) - if fninv and not ask_default: - self.read_metadata_thread(fninv) + return + set_inv(settings) + def calc_magnitude(self, type='ML'): self.init_metadata() if not self.metadata: @@ -3197,9 +3322,11 @@ class MainWindow(QMainWindow): self.fill_eventbox() self.getPlotWidget().draw() if self.plot_method == 'fast': - self.dataPlot.setPermText('MIN/MAX plot', color='red') + self.dataPlot.setPermText(1, ' MIN/MAX plot ', color='red') else: - self.dataPlot.setPermText() + self.dataPlot.setPermText(1) + self.dataPlot.setPermText(0, '| Number of traces: {} |'.format(len(self.getPlotWidget().getPlotDict()))) + def _setDirty(self): self.setDirty(True) diff --git a/autoPyLoT.py b/autoPyLoT.py index eb355609..5db95399 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -34,7 +34,7 @@ __version__ = _getVersionString() def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, eventid=None, savepath=None, - savexml=True, station='all', iplot=0, ncores=0): + savexml=True, station='all', iplot=0, ncores=0, obspyDMT_wfpath=False): """ Determine phase onsets automatically utilizing the automatic picking algorithms by Kueperkoch et al. 2010/2012. @@ -94,7 +94,6 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even fig_dict = None fig_dict_wadatijack = None - locflag = 1 if input_dict and isinstance(input_dict, dict): if 'parameter' in input_dict: parameter = input_dict['parameter'] @@ -110,10 +109,10 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even eventid = input_dict['eventid'] if 'iplot' in input_dict: iplot = input_dict['iplot'] - if 'locflag' in input_dict: - locflag = input_dict['locflag'] if 'savexml' in input_dict: savexml = input_dict['savexml'] + if 'obspyDMT_wfpath' in input_dict: + obspyDMT_wfpath = input_dict['obspyDMT_wfpath'] if not parameter: if inputfile: @@ -151,7 +150,8 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even datastructure.setExpandFields(exf) # check if default location routine NLLoc is available - if real_None(parameter['nllocbin']) and locflag: + if real_None(parameter['nllocbin']): + locflag = 1 # get NLLoc-root path nllocroot = parameter.get('nllocroot') # get path to NLLoc executable @@ -167,7 +167,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even ttpat = parameter.get('ttpatter') # pattern of NLLoc-output file nllocoutpatter = parameter.get('outpatter') - maxnumit = 3 # maximum number of iterations for re-picking + maxnumit = 2 # maximum number of iterations for re-picking else: locflag = 0 print(" !!! ") @@ -175,6 +175,14 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even print("!!No source parameter estimation possible!!") print(" !!! ") + wfpath_extension = '' + if obspyDMT_wfpath not in [None, False, 'False', '']: + wfpath_extension = obspyDMT_wfpath + print('Using obspyDMT structure. There will be no restitution, as pre-processed data are expected.') + if wfpath_extension != 'processed': + print('WARNING: Expecting wfpath_extension to be "processed" for' + ' pre-processed data but received "{}" instead!!!'.format(wfpath_extension)) + if not input_dict: # started in production mode datapath = datastructure.expandDataPath() @@ -219,6 +227,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even glocflag = locflag for eventpath in events: evID = os.path.split(eventpath)[-1] + event_datapath = os.path.join(eventpath, wfpath_extension) fext = '.xml' filename = os.path.join(eventpath, 'PyLoT_' + evID + fext) try: @@ -231,7 +240,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, eventpath, '*'))) + data.setWFData(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 @@ -268,12 +277,18 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even wfdat = check4doubled(wfdat) wfdat = trim_station_components(wfdat, trim_start=True, trim_end=False) metadata = read_metadata(parameter.get('invdir')) - # rotate stations to ZNE - wfdat = check4rotated(wfdat, metadata) + # TODO: (idea) read metadata from obspy_dmt database + # if not wfpath_extension: + # metadata = read_metadata(parameter.get('invdir')) + # else: + # metadata = None corr_dat = None - if locflag: - print("Restitute data ...") - corr_dat = restitute_data(wfdat.copy(), *metadata, ncores=ncores) + if metadata: + # rotate stations to ZNE + wfdat = check4rotated(wfdat, metadata) + if locflag: + print("Restitute data ...") + corr_dat = restitute_data(wfdat.copy(), *metadata, ncores=ncores) if not corr_dat and locflag: locflag = 2 print('Working on event %s. Stations: %s' % (eventpath, station)) @@ -331,7 +346,8 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even picks[stats]['P'].update(props) evt = moment_mag.updated_event() net_mw = moment_mag.net_magnitude() - print("Network moment magnitude: %4.1f" % net_mw.mag) + if net_mw is not None: + print("Network moment magnitude: %4.1f" % net_mw.mag) # calculate local (Richter) magntiude WAscaling = parameter.get('WAscaling') magscaling = parameter.get('magscaling') @@ -348,8 +364,15 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even evt = local_mag.updated_event(magscaling) net_ml = local_mag.net_magnitude(magscaling) print("Network local magnitude: %4.1f" % net_ml.mag) - print("Network local magnitude scaled with:") - print("%f * Ml + %f" % (magscaling[0], magscaling[1])) + if magscaling == None: + scaling = False + elif magscaling[0] != 0 and magscaling[1] != 0: + scaling = False + else: + scaling = True + if scaling: + print("Network local magnitude scaled with:") + print("%f * Ml + %f" % (magscaling[0], magscaling[1])) else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") @@ -406,7 +429,8 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even picks[stats]['P'].update(props) evt = moment_mag.updated_event() net_mw = moment_mag.net_magnitude() - print("Network moment magnitude: %4.1f" % net_mw.mag) + if net_mw is not None: + print("Network moment magnitude: %4.1f" % net_mw.mag) # calculate local (Richter) magntiude WAscaling = parameter.get('WAscaling') magscaling = parameter.get('magscaling') @@ -424,8 +448,15 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even evt = local_mag.updated_event(magscaling) net_ml = local_mag.net_magnitude(magscaling) print("Network local magnitude: %4.1f" % net_ml.mag) - print("Network local magnitude scaled with:") - print("%f * Ml + %f" % (magscaling[0], magscaling[1])) + if magscaling == None: + scaling = False + elif magscaling[0] != 0 and magscaling[1] != 0: + scaling = False + else: + scaling = True + if scaling: + print("Network local magnitude scaled with:") + print("%f * Ml + %f" % (magscaling[0], magscaling[1])) else: print("autoPyLoT: No NLLoc-location file available! Stop iteration!") locflag = 9 @@ -512,9 +543,12 @@ if __name__ == "__main__": parser.add_argument('-c', '-C', '--ncores', type=int, action='store', default=0, help='''optional, number of CPU cores used for parallel processing (default: all available(=0))''') + parser.add_argument('-dmt', '-DMT', '--obspy_dmt_wfpath', type=str, + action='store', default=False, + help='''optional, wftype (raw, processed) used for obspyDMT database structure''') cla = parser.parse_args() picks = autoPyLoT(inputfile=str(cla.inputfile), fnames=str(cla.fnames), eventid=str(cla.eventid), savepath=str(cla.spath), - ncores=cla.ncores, iplot=int(cla.iplot)) + ncores=cla.ncores, iplot=int(cla.iplot), obspyDMT_wfpath=str(cla.obspy_dmt_wfpath)) diff --git a/icons_rc_2.py b/icons_rc_2.py index bcf764af..9f36b183 100644 --- a/icons_rc_2.py +++ b/icons_rc_2.py @@ -7,7 +7,7 @@ # # WARNING! All changes made in this file will be lost! -from PyQt4 import QtCore +from PySide import QtCore qt_resource_data = "\ \x00\x00\x9e\x04\ diff --git a/icons_rc_3.py b/icons_rc_3.py index a437650d..963cf353 100644 --- a/icons_rc_3.py +++ b/icons_rc_3.py @@ -7,7 +7,7 @@ # # WARNING! All changes made in this file will be lost! -from PyQt4 import QtCore +from PySide import QtCore qt_resource_data = "\ \x00\x00\x9e\x04\ diff --git a/inputs/pylot_local.in b/inputs/pylot_local.in index fe019440..ccfeddd4 100644 --- a/inputs/pylot_local.in +++ b/inputs/pylot_local.in @@ -4,19 +4,19 @@ %Parameters are optimized for %extent data sets! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #main settings# - #rootpath# %project path - #datapath# %data path - #database# %name of data base - #eventID# %event ID for single event processing (* for all events found in database) - #invdir# %full path to inventory or dataless-seed file +/DATA/Insheim #rootpath# %project path +EVENT_DATA/LOCAL #datapath# %data path +2018.02_Insheim #database# %name of data base +e0006.038.18 #eventID# %event ID for single event processing (* for all events found in database) +/DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure# %choose data structure True #apverbose# %choose 'True' or 'False' for terminal output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #NLLoc settings# -None #nllocbin# %path to NLLoc executable -None #nllocroot# %root of NLLoc-processing directory -None #phasefile# %name of autoPyLoT-output phase file for NLLoc -None #ctrfile# %name of autoPyLoT-output control file for NLLoc +/home/ludger/NLLOC #nllocbin# %path to NLLoc executable +/home/ludger/NLLOC/Insheim #nllocroot# %root of NLLoc-processing directory +AUTOPHASES.obs #phasefile# %name of autoPyLoT-output phase file for NLLoc +Insheim_min1d032016_auto.in #ctrfile# %name of autoPyLoT-output control file for NLLoc ttime #ttpatter# %pattern of NLLoc ttimes from grid AUTOLOC_nlloc #outpatter# %pattern of NLLoc-output file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -27,31 +27,31 @@ AUTOLOC_nlloc #outpatter# %pattern of %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #settings local magnitude# 1.11 0.0009 -2.0 #WAscaling# %Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm] If zeros are set, original Richter magnitude is calculated! -1.0382 -0.447 #magscaling# %Scaling relation for derived local magnitude [a*Ml+b]. If zeros are set, no scaling of network magnitude is applied! +0.0 0.0 #magscaling# %Scaling relation for derived local magnitude [a*Ml+b]. If zeros are set, no scaling of network magnitude is applied! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #filter settings# -1.0 1.0 #minfreq# %Lower filter frequency [P, S] -10.0 10.0 #maxfreq# %Upper filter frequency [P, S] -2 2 #filter_order# %filter order [P, S] +2.0 2.0 #minfreq# %Lower filter frequency [P, S] +30.0 15.0 #maxfreq# %Upper filter frequency [P, S] +3 3 #filter_order# %filter order [P, S] bandpass bandpass #filter_type# %filter type (bandpass, bandstop, lowpass, highpass) [P, S] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #common settings picker# local #extent# %extent of array ("local", "regional" or "global") -15.0 #pstart# %start time [s] for calculating CF for P-picking -60.0 #pstop# %end time [s] for calculating CF for P-picking --1.0 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking +7.0 #pstart# %start time [s] for calculating CF for P-picking +16.0 #pstop# %end time [s] for calculating CF for P-picking +-0.5 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking 10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking -True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF +False #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF iasp91 #taup_model# %define TauPy model for traveltime estimation -2.0 10.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] -2.0 12.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] -2.0 8.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] -2.0 10.0 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz] +2.0 20.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] +2.0 30.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] +2.0 10.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] +2.0 15.0 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz] #special settings for calculating CF# %!!Edit the following only if you know what you are doing!!% #Z-component# HOS #algoP# %choose algorithm for P-onset determination (HOS, ARZ, or AR3) -7.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s] +4.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s] 4 #hosorder# %for HOS-picker, order of Higher Order Statistics 2 #Parorder# %for AR-picker, order of AR process of Z-component 1.2 #tdet1z# %for AR-picker, length of AR determination window [s] for Z-component, 1st pick @@ -59,12 +59,12 @@ HOS #algoP# %choose algo 0.6 #tdet2z# %for AR-picker, length of AR determination window [s] for Z-component, 2nd pick 0.2 #tpred2z# %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick 0.001 #addnoise# %add noise to seismogram for stable AR prediction -3.0 0.1 0.5 1.0 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s] +3.0 0.0 1.0 0.5 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s] 3.0 #pickwinP# %for initial AIC pick, length of P-pick window [s] 6.0 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick) -0.2 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s] +0.4 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s] 0.1 #tsmoothP# %for HOS/AR, take average of samples for smoothing CF [s] -0.001 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P) +0.4 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P) 1.3 #nfacP# %for HOS/AR, noise factor for noise level determination (P) #H-components# ARH #algoS# %choose algorithm for S-onset determination (ARH or AR3) @@ -75,7 +75,7 @@ ARH #algoS# %choose algo 4 #Sarorder# %for AR-picker, order of AR process of H-components 5.0 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H) 4.0 #pickwinS# %for initial AIC pick, length of S-pick window [s] -2.0 0.3 1.5 1.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s] +2.0 0.2 1.5 1.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s] 1.0 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s] 0.7 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S) 0.9 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S) @@ -85,16 +85,16 @@ ARH #algoS# %choose algo 2.0 #minFMSNR# %miniumum required SNR for first-motion determination 0.2 #fmpickwin# %pick window around P onset for calculating zero crossings #quality assessment# -0.02 0.04 0.08 0.16 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P -0.04 0.08 0.16 0.32 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S +0.04 0.08 0.16 0.32 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P +0.05 0.10 0.20 0.40 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S 0.8 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected 1.1 #minAICPSNR# %below this SNR the initial P pick is rejected 1.0 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected 1.0 #minsiglength# %length of signal part for which amplitudes must exceed noiselevel [s] -1.0 #noisefactor# %noiselevel*noisefactor=threshold -10.0 #minpercent# %required percentage of amplitudes exceeding threshold -1.5 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude -6.0 #mdttolerance# %maximum allowed deviation of P picks from median [s] +1.1 #noisefactor# %noiselevel*noisefactor=threshold +50.0 #minpercent# %required percentage of amplitudes exceeding threshold +1.1 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude +5.0 #mdttolerance# %maximum allowed deviation of P picks from median [s] 1.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram -5.0 #jackfactor# %pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor +2.0 #jackfactor# %pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 1983bf7c..08704f7b 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -6,6 +6,7 @@ Revised/extended summer 2017. :author: Ludger Küperkoch / MAGS2 EP3 working group """ + import matplotlib.pyplot as plt import numpy as np import obspy.core.event as ope @@ -121,7 +122,13 @@ class Magnitude(object): def net_magnitude(self, magscaling=None): if self: - if magscaling is not None and str(magscaling) is not '[0.0, 0.0]': + if magscaling == None: + scaling = False + elif magscaling[0] != 0 and magscaling[1] != 0: + scaling = False + else: + scaling = True + if scaling: # scaling necessary print("Scaling network magnitude ...") mag = ope.Magnitude( @@ -140,7 +147,6 @@ class Magnitude(object): station_count=len(self.magnitudes), azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap) return mag - return None class LocalMagnitude(Magnitude): @@ -219,7 +225,7 @@ class LocalMagnitude(Magnitude): sqH = np.sqrt(power_sum) # get time array - th = np.arange(0, len(sqH) * dt, dt) + th=np.arange(0, st[0].stats.npts/st[0].stats.sampling_rate, st[0].stats.delta) # get maximum peak within pick window iwin = getsignalwin(th, t0 - stime, self.calc_win) ii = min([iwin[len(iwin) - 1], len(th)]) diff --git a/pylot/core/io/data.py b/pylot/core/io/data.py index a7cb85d5..5e6db468 100644 --- a/pylot/core/io/data.py +++ b/pylot/core/io/data.py @@ -15,7 +15,7 @@ from pylot.core.util.event import Event from pylot.core.util.utils import fnConstructor, full_range, remove_underscores, check4gaps, check4doubled, \ check4rotated, trim_station_components import pylot.core.loc.velest as velest - +from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT class Data(object): """ @@ -43,7 +43,7 @@ class Data(object): elif isinstance(evtdata, dict): evt = readPILOTEvent(**evtdata) evtdata = evt - elif isinstance(evtdata, str): + elif type(evtdata) in [str, unicode]: try: cat = read_events(evtdata) if len(cat) is not 1: @@ -60,6 +60,8 @@ class Data(object): raise NotImplementedError('PILOT location information ' 'read support not yet ' 'implemeted.') + elif 'event.pkl' in evtdata: + evtdata = qml_from_obspyDMT(evtdata) else: raise e else: @@ -72,6 +74,7 @@ class Data(object): self.wforiginal = None self.cuttimes = None self.dirty = False + self.processed = None def __str__(self): return str(self.wfdata) @@ -368,7 +371,7 @@ class Data(object): data.filter(**kwargs) self.dirty = True - def setWFData(self, fnames, checkRotated=False, metadata=None, obspy_dmt=False): + def setWFData(self, fnames, fnames_syn=None, checkRotated=False, metadata=None): """ Clear current waveform data and set given waveform data :param fnames: waveform data names to append @@ -377,30 +380,31 @@ class Data(object): self.wfdata = Stream() self.wforiginal = None self.wfsyn = Stream() - wffnames = None - wffnames_syn = None - wfdir = 'processed' if 'processed' in [fname.split('/')[-1] for fname in fnames] else 'raw' - if obspy_dmt: - 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 wffnames is not None: - self.appendWFData(wffnames) - if wffnames_syn is not None: - self.appendWFData(wffnames_syn, synthetic=True) + # 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 doubled channels - check4gaps(self.wfdata) - check4doubled(self.wfdata) # check for stations with rotated components if checkRotated and metadata is not None: self.wfdata = check4rotated(self.wfdata, metadata, verbosity=0) @@ -441,7 +445,7 @@ class Data(object): except SacIOError as se: warnmsg += '{0}\n{1}\n'.format(fname, se) if warnmsg: - warnmsg = 'WARNING: unable to read\n' + warnmsg + warnmsg = 'WARNING in appendWFData: unable to read waveform data\n' + warnmsg print(warnmsg) def getWFData(self): diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index a5f93d6b..9e28a08a 100644 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -100,9 +100,9 @@ def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None, # quality control # median check and jackknife on P-onset times - jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, 1, fig_dict_wadatijack) + jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, iplot, fig_dict_wadatijack) # check S-P times (Wadati) - wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, 1, fig_dict_wadatijack) + wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, iplot, fig_dict_wadatijack) return wadationsets @@ -1057,8 +1057,64 @@ def nautopickstation(wfstream, pickparam, verbose=False, # read your pylot.in for details! plt_flag = 0 - # get picking parameter dictionaries - p_params, s_params, first_motion_params, signal_length_params = get_pickparams(pickparam) + # special parameters for P picking + algoP = pickparam.get('algoP') + pstart = pickparam.get('pstart') + pstop = pickparam.get('pstop') + thosmw = pickparam.get('tlta') + tsnrz = pickparam.get('tsnrz') + hosorder = pickparam.get('hosorder') + bpz1 = pickparam.get('bpz1') + bpz2 = pickparam.get('bpz2') + pickwinP = pickparam.get('pickwinP') + aictsmoothP = pickparam.get('aictsmooth') + tsmoothP = pickparam.get('tsmoothP') + ausP = pickparam.get('ausP') + nfacP = pickparam.get('nfacP') + tpred1z = pickparam.get('tpred1z') + tdet1z = pickparam.get('tdet1z') + tpred2z = pickparam.get('tpred2z') + tdet2z = pickparam.get('tdet2z') + Parorder = pickparam.get('Parorder') + addnoise = pickparam.get('addnoise') + Precalcwin = pickparam.get('Precalcwin') + minAICPslope = pickparam.get('minAICPslope') + minAICPSNR = pickparam.get('minAICPSNR') + timeerrorsP = pickparam.get('timeerrorsP') + # special parameters for S picking + algoS = pickparam.get('algoS') + sstart = pickparam.get('sstart') + sstop = pickparam.get('sstop') + use_taup = real_Bool(pickparam.get('use_taup')) + taup_model = pickparam.get('taup_model') + bph1 = pickparam.get('bph1') + bph2 = pickparam.get('bph2') + tsnrh = pickparam.get('tsnrh') + pickwinS = pickparam.get('pickwinS') + tpred1h = pickparam.get('tpred1h') + tdet1h = pickparam.get('tdet1h') + tpred2h = pickparam.get('tpred2h') + tdet2h = pickparam.get('tdet2h') + Sarorder = pickparam.get('Sarorder') + aictsmoothS = pickparam.get('aictsmoothS') + tsmoothS = pickparam.get('tsmoothS') + ausS = pickparam.get('ausS') + minAICSslope = pickparam.get('minAICSslope') + minAICSSNR = pickparam.get('minAICSSNR') + Srecalcwin = pickparam.get('Srecalcwin') + nfacS = pickparam.get('nfacS') + timeerrorsS = pickparam.get('timeerrorsS') + # parameters for first-motion determination + minFMSNR = pickparam.get('minFMSNR') + fmpickwin = pickparam.get('fmpickwin') + minfmweight = pickparam.get('minfmweight') + # parameters for checking signal length + minsiglength = pickparam.get('minsiglength') + minpercent = pickparam.get('minpercent') + nfacsl = pickparam.get('noisefactor') + # parameter to check for spuriously picked S onset + zfac = pickparam.get('zfac') + # path to inventory-, dataless- or resp-files # initialize output Pweight = 4 # weight for P onset @@ -1159,6 +1215,8 @@ def nautopickstation(wfstream, pickparam, verbose=False, if p_params['use_taup'] is True: Lc = np.inf print('autopickstation: use_taup flag active.') + if not metadata: + metadata = [None, None] if not metadata[1]: print('Warning: Could not use TauPy to estimate onsets as there are no metadata given.') else: @@ -1207,7 +1265,7 @@ def nautopickstation(wfstream, pickparam, verbose=False, return Ldiff = Lwf - abs(Lc) - if Ldiff < 0 or pstop <= pstart: + if Ldiff <= 0 or pstop <= pstart or pstop - pstart <= thosmw: msg = 'autopickstation: Cutting times are too large for actual ' \ 'waveform!\nUsing entire waveform instead!' if verbose: print(msg) @@ -1345,8 +1403,8 @@ def nautopickstation(wfstream, pickparam, verbose=False, elif p_params['algoP'] == 'ARZ': # calculate ARZ-CF using subclass ARZcf of class # CharcteristicFunction - cf2 = ARZcf(z_copy, cuttimes2, p_params['tpred1z'], p_params['Parorder'], p_params['tdet1z'], - p_params['addnoise']) # instance of ARZcf + cf2 = ARZcf(z_copy, cuttimes2, tpred2z, Parorder, tdet2z, + addnoise) # instance of ARZcf ############################################################## # get refined onset time from CF2 using class Picker assert isinstance(cf2, CharacteristicFunction), 'cf2 is not set ' \ @@ -1401,7 +1459,8 @@ def nautopickstation(wfstream, pickparam, verbose=False, msg = "autopickstation: P-weight: {0}, " \ "SNR: {1}, SNR[dB]: {2}, Polarity: {3}".format(Pweight, SNRP, SNRPdB, FM) print(msg) - msg = 'autopickstation: Refined P-Pick: {} s | P-Error: {} s'.format(mpickP, Perror) + msg = 'autopickstation: Refind P-Pick: {} s | P-Error: {} s'.format(zdat[0].stats.starttime \ + + mpickP, Perror) print(msg) Sflag = 1 @@ -1667,7 +1726,8 @@ def nautopickstation(wfstream, pickparam, verbose=False, lpickS = lpick[ipick] Serror = pickerr[ipick] - msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(mpickS, Serror) + msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(hdat[0].stats.starttime \ + + mpickS, Serror) print(msg) # get SNR diff --git a/pylot/core/pick/charfuns.py b/pylot/core/pick/charfuns.py index a62870eb..f3aa8a64 100644 --- a/pylot/core/pick/charfuns.py +++ b/pylot/core/pick/charfuns.py @@ -78,7 +78,7 @@ class CharacteristicFunction(object): return self.cut def setCut(self, cut): - self.cut = cut + self.cut = (int(cut[0]), int(cut[1])) def getTime1(self): return self.t1 diff --git a/pylot/core/pick/picker.py b/pylot/core/pick/picker.py index c10e2420..478c3ce2 100644 --- a/pylot/core/pick/picker.py +++ b/pylot/core/pick/picker.py @@ -233,12 +233,41 @@ class AICPicker(AutoPicker): # quality assessment using SNR and slope from CF if self.Pick is not None: - self.Data[0].data = check_counts_ms(self.Data[0].data) + # get noise window + inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1]) + # check, if these are counts or m/s, important for slope estimation! + # this is quick and dirty, better solution? #Todo wtf + if max(self.Data[0].data < 1e-3) and max(self.Data[0].data >= 1e-6): + self.Data[0].data = self.Data[0].data * 1000000. + elif max(self.Data[0].data < 1e-6): + self.Data[0].data = self.Data[0].data * 1e13 + # get signal window + isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2]) + if len(isignal) == 0: + return + ii = min([isignal[len(isignal) - 1], len(self.Tcf)]) + isignal = isignal[0:ii] + try: + self.Data[0].data[isignal] + except IndexError as e: + msg = "Time series out of bounds! {}".format(e) + print(msg) + return # calculate SNR from CF - # subtract beginning of Tcf to get sample index in self.Data, which - # is cut out off trace - self.SNR = getSNR(self.Data, self.TSNR, self.Pick-self.Tcf[0])[0] # TODO Check wether this yields similar results + self.SNR = max(abs(self.Data[0].data[isignal])) / \ + abs(np.mean(self.Data[0].data[inoise])) # calculate slope from CF after initial pick + # get slope window + tslope = self.TSNR[3] # slope determination window + tsafety = self.TSNR[1] # safety gap, AIC is usually a little bit too late + if tsafety >= 0: + islope = np.where((self.Tcf <= min([self.Pick + tslope + tsafety, self.Tcf[-1]])) \ + & (self.Tcf >= self.Pick)) # TODO: put this in a seperate function like getsignalwin + else: + islope = np.where((self.Tcf <= min([self.Pick + tslope, self.Tcf[-1]])) \ + & (self.Tcf >= self.Pick + tsafety)) # TODO: put this in a seperate function like getsignalwin + # find maximum within slope determination window + # 'cause slope should be calculated up to first local minimum only! try: self.slope, iislope, datafit = calcSlope(self.Data, self.aicsmooth, self.Tcf, self.Pick, self.TSNR) except Exception: @@ -264,6 +293,56 @@ class AICPicker(AutoPicker): pass plt.close(fig) return + if len(dataslope) < 2: + print('No or not enough data in slope window found!') + return + try: + imaxs, = argrelmax(dataslope) + imaxs.size + imax = imaxs[0] + except ValueError as e: + print(e, 'picker: argrelmax not working!') + imax = np.argmax(dataslope) + iislope = islope[0][0:imax + 1] + if len(iislope) < 2: + # calculate slope from initial onset to maximum of AIC function + print("AICPicker: Not enough data samples left for slope calculation!") + print("Calculating slope from initial onset to maximum of AIC function ...") + imax = np.argmax(aicsmooth[islope[0][0:-1]]) + if imax == 0: + print("AICPicker: Maximum for slope determination right at the beginning of the window!") + print("Choose longer slope determination window!") + if self.iplot > 1: + if self.fig == None or self.fig == 'None': + fig = plt.figure() + plt_flag = 1 + else: + fig = self.fig + ax = fig.add_subplot(111) + x = self.Data[0].data + ax.plot(self.Tcf, x / max(x), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data') + ax.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF') + ax.legend(loc=1) + ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime) + ax.set_yticks([]) + ax.set_title(self.Data[0].stats.station) + if plt_flag == 1: + fig.show() + try: input() + except SyntaxError: pass + plt.close(fig) + return + iislope = islope[0][0:imax+1] + dataslope = self.Data[0].data[iislope] + # calculate slope as polynomal fit of order 1 + xslope = np.arange(0, len(dataslope), 1) + P = np.polyfit(xslope, dataslope, 1) + datafit = np.polyval(P, xslope) + if datafit[0] >= datafit[-1]: + print('AICPicker: Negative slope, bad onset skipped!') + else: + self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0]) + else: self.SNR = -1 self.slope = -1 @@ -305,24 +384,23 @@ class AICPicker(AutoPicker): label='Slope Window') ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2, label='Slope') - ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, - self.SNR, self.slope)) + if self.slope is not None: + ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, + self.SNR, self.slope)) + else: + ax1.set_title('Station %s, SNR=%7.2f' % (self.Data[0].stats.station, self.SNR)) ax2.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime) ax2.set_ylabel('Counts') ax2.set_yticks([]) ax2.legend(loc=1) - if plt_flag == 1: - fig.show() - try: input() - except SyntaxError: pass - plt.close(fig) else: ax1.set_title(self.Data[0].stats.station) - if plt_flag == 1: - fig.show() - try: input() - except SyntaxError: pass - plt.close(fig) + + if plt_flag == 1: + fig.show() + try: input() + except SyntaxError: pass + plt.close(fig) if self.Pick == None: print('AICPicker: Could not find minimum, picking window too short?') diff --git a/pylot/core/util/array_map.py b/pylot/core/util/array_map.py new file mode 100644 index 00000000..de9a3bf3 --- /dev/null +++ b/pylot/core/util/array_map.py @@ -0,0 +1,458 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import matplotlib.pyplot as plt +import numpy as np +import obspy +from PySide import QtGui +from mpl_toolkits.basemap import Basemap +from matplotlib.figure import Figure +from mpl_toolkits.axes_grid1.inset_locator import inset_axes +from pylot.core.util.widgets import PickDlg, PylotCanvas +from scipy.interpolate import griddata + +plt.interactive(False) + + +class Array_map(QtGui.QWidget): + def __init__(self, parent, figure=None): + ''' + Create a map of the array. + :param parent: PyLoT Mainwindow class + :param figure: + ''' + QtGui.QWidget.__init__(self) + self._parent = parent + self.metadata_type = parent.metadata[0] + self.metadata = parent.metadata[1] + self.picks = None + self.picks_dict = None + self.autopicks_dict = None + self.eventLoc = None + self.figure = figure + self.init_graphics() + self.init_stations() + self.init_basemap(resolution='l') + self.init_map() + self._style = parent._style + # self.show() + + def init_map(self): + self.init_lat_lon_dimensions() + self.init_lat_lon_grid() + self.init_x_y_dimensions() + self.connectSignals() + self.draw_everything() + self.canvas.setZoomBorders2content() + + def onpick(self, event): + ind = event.ind + button = event.mouseevent.button + if ind == [] or not button == 1: + return + data = self._parent.get_data().getWFData() + for index in ind: + station = str(self._station_onpick_ids[index].split('.')[-1]) + try: + data = data.select(station=station) + if not data: + self._warn('No data for station {}'.format(station)) + return + pickDlg = PickDlg(self._parent, parameter=self._parent._inputs, + data=data, + station=station, + picks=self._parent.get_current_event().getPick(station), + autopicks=self._parent.get_current_event().getAutopick(station), + filteroptions=self._parent.filteroptions) + except Exception as e: + message = 'Could not generate Plot for station {st}.\n {er}'.format(st=station, er=e) + self._warn(message) + print(message, e) + return + pyl_mw = self._parent + try: + if pickDlg.exec_(): + pyl_mw.setDirty(True) + pyl_mw.update_status('picks accepted ({0})'.format(station)) + replot = pyl_mw.get_current_event().setPick(station, pickDlg.getPicks()) + self._refresh_drawings() + if replot: + pyl_mw.plotWaveformData() + pyl_mw.drawPicks() + pyl_mw.draw() + else: + pyl_mw.drawPicks(station) + pyl_mw.draw() + else: + pyl_mw.update_status('picks discarded ({0})'.format(station)) + except Exception as e: + message = 'Could not save picks for station {st}.\n{er}'.format(st=station, er=e) + self._warn(message) + print(message, e) + + def connectSignals(self): + self.comboBox_phase.currentIndexChanged.connect(self._refresh_drawings) + self.comboBox_am.currentIndexChanged.connect(self._refresh_drawings) + self.canvas.mpl_connect('motion_notify_event', self.mouse_moved) + #self.zoom_id = self.basemap.ax.figure.canvas.mpl_connect('scroll_event', self.zoom) + + def _from_dict(self, function, key): + return function(self.stations_dict.values(), key=lambda x: x[key])[key] + + def get_min_from_stations(self, key): + return self._from_dict(min, key) + + def get_max_from_stations(self, key): + return self._from_dict(max, key) + + def get_min_from_picks(self): + return min(self.picks_rel.values()) + + def get_max_from_picks(self): + return max(self.picks_rel.values()) + + def mouse_moved(self, event): + if not event.inaxes == self.main_ax: + return + x = event.xdata + y = event.ydata + lat, lon = self.basemap(x, y, inverse=True) + self.status_label.setText('Latitude: {}, Longitude: {}'.format(lat, lon)) + + def current_picks_dict(self): + picktype = self.comboBox_am.currentText() + auto_manu = {'auto': self.autopicks_dict, + 'manual': self.picks_dict} + return auto_manu[picktype] + + def init_graphics(self): + if not self.figure: + self.figure = Figure() + + self.status_label = QtGui.QLabel() + + self.main_ax = self.figure.add_subplot(111) + self.canvas = PylotCanvas(self.figure, parent=self._parent, multicursor=True, + panZoomX=False, panZoomY=False) + + self.main_box = QtGui.QVBoxLayout() + self.setLayout(self.main_box) + + self.top_row = QtGui.QHBoxLayout() + self.main_box.addLayout(self.top_row, 1) + + self.comboBox_phase = QtGui.QComboBox() + self.comboBox_phase.insertItem(0, 'P') + self.comboBox_phase.insertItem(1, 'S') + + self.comboBox_am = QtGui.QComboBox() + self.comboBox_am.insertItem(0, 'auto') + self.comboBox_am.insertItem(1, 'manual') + + self.top_row.addWidget(QtGui.QLabel('Select a phase: ')) + self.top_row.addWidget(self.comboBox_phase) + self.top_row.setStretch(1, 1) # set stretch of item 1 to 1 + self.top_row.addWidget(QtGui.QLabel('Pick type: ')) + self.top_row.addWidget(self.comboBox_am) + self.top_row.setStretch(3, 1) # set stretch of item 1 to 1 + + self.main_box.addWidget(self.canvas, 1) + self.main_box.addWidget(self.status_label, 0) + + + def init_stations(self): + def stat_info_from_parser(parser): + stations_dict = {} + for station in parser.stations: + station_name = station[0].station_call_letters + network_name = station[0].network_code + if not station_name in stations_dict.keys(): + st_id = network_name + '.' + station_name + stations_dict[st_id] = {'latitude': station[0].latitude, + 'longitude': station[0].longitude} + return stations_dict + + def stat_info_from_inventory(inventory): + stations_dict = {} + for network in inventory.networks: + for station in network.stations: + station_name = station.code + network_name = network_name.code + if not station_name in stations_dict.keys(): + st_id = network_name + '.' + station_name + stations_dict[st_id] = {'latitude': station[0].latitude, + 'longitude': station[0].longitude} + return stations_dict + + read_stat = {'xml': stat_info_from_inventory, + 'dless': stat_info_from_parser} + + self.stations_dict = read_stat[self.metadata_type](self.metadata) + self.latmin = self.get_min_from_stations('latitude') + self.lonmin = self.get_min_from_stations('longitude') + self.latmax = self.get_max_from_stations('latitude') + self.lonmax = self.get_max_from_stations('longitude') + + def init_picks(self): + def get_picks(station_dict): + picks = {} + phase = self.comboBox_phase.currentText() + for st_id in station_dict.keys(): + try: + station_name = st_id.split('.')[-1] + pick = self.current_picks_dict()[station_name][phase] + if pick['picker'] == 'auto': + if pick['weight'] > 3: + continue + picks[st_id] = pick['mpp'] + except KeyError: + continue + except Exception as e: + print('Cannot display pick for station {}. Reason: {}'.format(station_name, e)) + return picks + + def get_picks_rel(picks): + picks_rel = {} + picks_utc = [] + for pick in picks.values(): + if type(pick) is obspy.core.utcdatetime.UTCDateTime: + picks_utc.append(pick) + self._earliest_picktime = min(picks_utc) + for st_id, pick in picks.items(): + if type(pick) is obspy.core.utcdatetime.UTCDateTime: + pick -= self._earliest_picktime + picks_rel[st_id] = pick + return picks_rel + + self.picks = get_picks(self.stations_dict) + self.picks_rel = get_picks_rel(self.picks) + + def init_lat_lon_dimensions(self): + # init minimum and maximum lon and lat dimensions + self.londim = self.lonmax - self.lonmin + self.latdim = self.latmax - self.latmin + + def init_x_y_dimensions(self): + # transformation of lat/lon to ax coordinate system + for st_id, coords in self.stations_dict.items(): + lat, lon = coords['latitude'], coords['longitude'] + coords['x'], coords['y'] = self.basemap(lon, lat) + + self.xdim = self.get_max_from_stations('x') - self.get_min_from_stations('x') + self.ydim = self.get_max_from_stations('y') - self.get_min_from_stations('y') + + def init_basemap(self, resolution='l'): + # basemap = Basemap(projection=projection, resolution = resolution, ax=self.main_ax) + width = 5e6 + height = 2e6 + basemap = Basemap(projection='lcc', resolution=resolution, ax=self.main_ax, + width=width, height=height, + lat_0=(self.latmin + self.latmax) / 2., + lon_0=(self.lonmin + self.lonmax) / 2.) + + # basemap.fillcontinents(color=None, lake_color='aqua',zorder=1) + basemap.drawmapboundary(zorder=2) # fill_color='darkblue') + basemap.shadedrelief(zorder=3) + basemap.drawcountries(zorder=4) + basemap.drawstates(zorder=5) + basemap.drawcoastlines(zorder=6) + # labels = [left,right,top,bottom] + parallels = np.arange(-90, 90, 5.) + parallels_small = np.arange(-90, 90, 2.5) + basemap.drawparallels(parallels_small, linewidth=0.5, zorder=7) + basemap.drawparallels(parallels, zorder=7) + meridians = np.arange(-180, 180, 5.) + meridians_small = np.arange(-180, 180, 2.5) + basemap.drawmeridians(meridians_small, linewidth=0.5, zorder=7) + basemap.drawmeridians(meridians, zorder=7) + self.basemap = basemap + self.figure._tight = True + self.figure.tight_layout() + + def init_lat_lon_grid(self, nstep=250): + # create a regular grid to display colormap + lataxis = np.linspace(self.latmin, self.latmax, nstep) + lonaxis = np.linspace(self.lonmin, self.lonmax, nstep) + self.longrid, self.latgrid = np.meshgrid(lonaxis, lataxis) + + def init_picksgrid(self): + picks, lats, lons = self.get_picks_lat_lon() + try: + self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid), + method='linear') + except Exception as e: + self._warn('Could not init picksgrid: {}'.format(e)) + + def get_st_lat_lon_for_plot(self): + stations = [] + latitudes = [] + longitudes = [] + for st_id, coords in self.stations_dict.items(): + stations.append(st_id) + latitudes.append(coords['latitude']) + longitudes.append(coords['longitude']) + return stations, latitudes, longitudes + + def get_st_x_y_for_plot(self): + stations = [] + xs = [] + ys = [] + for st_id, coords in self.stations_dict.items(): + stations.append(st_id) + xs.append(coords['x']) + ys.append(coords['y']) + return stations, xs, ys + + def get_picks_lat_lon(self): + picks = [] + latitudes = [] + longitudes = [] + for st_id, pick in self.picks_rel.items(): + picks.append(pick) + latitudes.append(self.stations_dict[st_id]['latitude']) + longitudes.append(self.stations_dict[st_id]['longitude']) + return picks, latitudes, longitudes + + def draw_contour_filled(self, nlevel='50'): + levels = np.linspace(self.get_min_from_picks(), self.get_max_from_picks(), nlevel) + self.contourf = self.basemap.contourf(self.longrid, self.latgrid, self.picksgrid_active, + levels, latlon=True, zorder=9, alpha=0.5) + + def scatter_all_stations(self): + stations, lats, lons = self.get_st_lat_lon_for_plot() + self.sc = self.basemap.scatter(lons, lats, s=50, facecolor='none', latlon=True, + zorder=10, picker=True, edgecolor='m', label='Not Picked') + self.cid = self.canvas.mpl_connect('pick_event', self.onpick) + self._station_onpick_ids = stations + if self.eventLoc: + lats, lons = self.eventLoc + self.sc_event = self.basemap.scatter(lons, lats, s=100, facecolor='red', + latlon=True, zorder=11, label='Event (might be outside map region)') + + def scatter_picked_stations(self): + picks, lats, lons = self.get_picks_lat_lon() + # workaround because of an issue with latlon transformation of arrays with len <3 + if len(lons) <= 2 and len(lats) <= 2: + self.sc_picked = self.basemap.scatter(lons[0], lats[0], s=50, facecolor='white', + c=picks[0], latlon=True, zorder=11, label='Picked') + if len(lons) == 2 and len(lats) == 2: + self.sc_picked = self.basemap.scatter(lons[1], lats[1], s=50, facecolor='white', + c=picks[1], latlon=True, zorder=11) + else: + self.sc_picked = self.basemap.scatter(lons, lats, s=50, facecolor='white', + c=picks, latlon=True, zorder=11, label='Picked') + + def annotate_ax(self): + self.annotations = [] + stations, xs, ys = self.get_st_x_y_for_plot() + for st, x, y in zip(stations, xs, ys): + self.annotations.append(self.main_ax.annotate(' %s' % st, xy=(x, y), + fontsize='x-small', color='white', zorder=12)) + self.legend = self.main_ax.legend(loc=1) + self.legend.get_frame().set_facecolor((1, 1, 1, 0.75)) + + def add_cbar(self, label): + self.cbax_bg = inset_axes(self.main_ax, width="6%", height="75%", loc=5) + cbax = inset_axes(self.main_ax, width='2%', height='70%', loc=5) + cbar = self.main_ax.figure.colorbar(self.sc_picked, cax = cbax) + cbar.set_label(label) + cbax.yaxis.tick_left() + cbax.yaxis.set_label_position('left') + for spine in self.cbax_bg.spines.values(): + spine.set_visible(False) + self.cbax_bg.yaxis.set_ticks([]) + self.cbax_bg.xaxis.set_ticks([]) + self.cbax_bg.patch.set_facecolor((1, 1, 1, 0.75)) + return cbar + + def refresh_drawings(self, picks=None, autopicks=None): + self.picks_dict = picks + self.autopicks_dict = autopicks + self._refresh_drawings() + + def _refresh_drawings(self): + self.remove_drawings() + self.draw_everything() + + def draw_everything(self): + if self.picks_dict or self.autopicks_dict: + self.init_picks() + if len(self.picks) >= 3: + self.init_picksgrid() + self.draw_contour_filled() + self.scatter_all_stations() + if self.picks_dict or self.autopicks_dict: + self.scatter_picked_stations() + self.cbar = self.add_cbar(label='Time relative to first onset ({}) [s]'.format(self._earliest_picktime)) + self.comboBox_phase.setEnabled(True) + else: + self.comboBox_phase.setEnabled(False) + self.annotate_ax() + self.canvas.draw() + + def remove_drawings(self): + if hasattr(self, 'cbar'): + self.cbar.remove() + self.cbax_bg.remove() + del (self.cbar, self.cbax_bg) + if hasattr(self, 'sc_picked'): + self.sc_picked.remove() + del (self.sc_picked) + if hasattr(self, 'sc_event'): + self.sc_event.remove() + del (self.sc_event) + if hasattr(self, 'contourf'): + self.remove_contourf() + del (self.contourf) + if hasattr(self, 'cid'): + self.canvas.mpl_disconnect(self.cid) + del (self.cid) + try: + self.sc.remove() + except Exception as e: + print('Warning: could not remove station scatter plot.\nReason: {}'.format(e)) + try: + self.legend.remove() + except Exception as e: + print('Warning: could not remove legend. Reason: {}'.format(e)) + self.canvas.draw() + + def remove_contourf(self): + for item in self.contourf.collections: + item.remove() + + def remove_annotations(self): + for annotation in self.annotations: + annotation.remove() + + def zoom(self, event): + map = self.basemap + xlim = map.ax.get_xlim() + ylim = map.ax.get_ylim() + x, y = event.xdata, event.ydata + zoom = {'up': 1. / 2., + 'down': 2.} + + if not event.xdata or not event.ydata: + return + + if event.button in zoom: + factor = zoom[event.button] + xdiff = (xlim[1] - xlim[0]) * factor + xl = x - 0.5 * xdiff + xr = x + 0.5 * xdiff + ydiff = (ylim[1] - ylim[0]) * factor + yb = y - 0.5 * ydiff + yt = y + 0.5 * ydiff + + if xl < map.xmin or yb < map.ymin or xr > map.xmax or yt > map.ymax: + xl, xr = map.xmin, map.xmax + yb, yt = map.ymin, map.ymax + map.ax.set_xlim(xl, xr) + map.ax.set_ylim(yb, yt) + map.ax.figure.canvas.draw() + + def _warn(self, message): + self.qmb = QtGui.QMessageBox(QtGui.QMessageBox.Icon.Warning, + 'Warning', message) + self.qmb.show() diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index 6b952164..ee74a960 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -12,6 +12,186 @@ from pylot.core.util.utils import key_for_set_value, find_in_list, \ remove_underscores, gen_Pool +class Metadata(object): + def __init__(self, inventory=None): + self.inventories = [] + if os.path.isdir(inventory): + self.add_inventory(inventory) + if os.path.isfile(inventory): + self.add_inventory_file(inventory) + self.seed_ids = {} + self.inventory_files = {} + + + def add_inventory(self, path_to_inventory): + # add paths to list of inventories + assert (os.path.isdir(path_to_inventory)), '{} is no directory'.format(path_to_inventory) + if not path_to_inventory in self.inventories: + self.inventories.append(path_to_inventory) + + + def add_inventory_file(self, path_to_inventory_file): + assert (os.path.isfile(path_to_inventory_file)), '{} is no directory'.format(path_to_inventory_file) + self.add_inventory(os.path.split(path_to_inventory_file)[0]) + if not path_to_inventory_file in self.inventory_files.keys(): + self.read_single_file(path_to_inventory_file) + + + def remove_inventory(self, path_to_inventory): + if not path_to_inventory in self.inventories: + print('Path {} not in inventories list.'.format(path_to_inventory)) + return + self.inventories.remove(path_to_inventory) + + + def get_metadata(self, seed_id): + # get metadata for a specific seed_id, if not already read, try to read from inventories + if not seed_id in self.seed_ids.keys(): + self._read_inventory_data(seed_id) + # if seed id is not found read all inventories and try to find it there + if not seed_id in self.seed_ids.keys(): + print('No data found for seed id {}. Trying to find it in all known inventories...'.format(seed_id)) + self.read_all() + for inv_fname, metadata in self.inventory_files.items(): + # use get_coordinates to check for seed_id + if metadata['data'].get_coordinates(seed_id): + self.seed_ids[seed_id] = inv_fname + return metadata + print('Could not find metadata for station {}'.format(seed_id)) + return None + fname = self.seed_ids[seed_id] + return self.inventory_files[fname] + + + def read_all(self): + for inventory in self.inventories: + for inv_fname in os.listdir(inventory): + inv_fname = os.path.join(inventory, inv_fname) + if not self.read_single_file(inv_fname): + continue + + + def read_single_file(self, inv_fname): + if not inv_fname in self.inventory_files.keys(): + pass + else: + if not self.inventory_files[inv_fname]: + pass + else: + return + try: + invtype, robj = self._read_metadata_file(inv_fname) + if robj == None: + return + except Exception as e: + print('Could not read file {}'.format(inv_fname)) + return + self.inventory_files[inv_fname] = {'invtype': invtype, + 'data': robj} + return True + + + def get_coordinates(self, seed_id): + metadata = self.get_metadata(seed_id) + return metadata['data'].get_coordinates(seed_id) + + + def get_paz(self, seed_id, time=None): + metadata = self.get_metadata(seed_id) + if metadata['invtype'] in ['dless', 'dseed']: + return metadata['data'].get_paz(seed_id) + elif metadata['invtype'] in ['resp', 'xml']: + if not time: + print('Time needed to extract metadata from station inventory.') + return None + resp = metadata['data'].get_response(seed_id, time) + return resp.get_paz(seed_id) + + + def _read_inventory_data(self, seed_id=None): + for inventory in self.inventories: + if self._read_metadata_iterator(path_to_inventory=inventory, station_seed_id=seed_id): + return + + + def _read_metadata_iterator(self, path_to_inventory, station_seed_id): + ''' + search for metadata for a specific station iteratively + ''' + station, network, location, channel = station_seed_id.split('.') + fnames = glob.glob(os.path.join(path_to_inventory, '*' + station_seed_id + '*')) + if not fnames: + # search for station name in filename + fnames = glob.glob(os.path.join(path_to_inventory, '*' + station + '*')) + if not fnames: + # search for network name in filename + fnames = glob.glob(os.path.join(path_to_inventory, '*' + network + '*')) + if not fnames: + print('Could not find filenames matching station name, network name or seed id') + return + for fname in fnames: + if fname in self.inventory_files.keys(): + if self.inventory_files[fname]: + # file already read + continue + invtype, robj = self._read_metadata_file(os.path.join(path_to_inventory, fname)) + try: + robj.get_coordinates(station_seed_id) + self.inventory_files[fname] = {'invtype': invtype, + 'data': robj} + if station_seed_id in self.seed_ids.keys(): + print('WARNING: Overwriting metadata for station {}'.format(station_seed_id)) + self.seed_ids[station_seed_id] = fname + return True + except Exception as e: + continue + print('Could not find metadata for station_seed_id {} in path {}'.format(station_seed_id, path_to_inventory)) + + + def _read_metadata_file(self, path_to_inventory_filename): + ''' + function reading metadata files (either dataless seed, xml or resp) + :param path_to_inventory_filename: + :return: file type/ending, inventory object (Parser or Inventory) + ''' + # functions used to read metadata for different file endings (or file types) + read_functions = {'dless': self._read_dless, + 'dseed': self._read_dless, + 'xml': self._read_inventory_file, + 'resp': self._read_inventory_file} + file_ending = path_to_inventory_filename.split('.')[-1] + if file_ending in read_functions.keys(): + robj, exc = read_functions[file_ending](path_to_inventory_filename) + if exc is not None: + raise exc + return file_ending, robj + # in case file endings did not match the above keys, try and error + for file_type in ['dless', 'xml']: + robj, exc = read_functions[file_type](path_to_inventory_filename) + if exc is None: + return file_type, robj + return None, None + + + def _read_dless(self, path_to_inventory): + exc = None + try: + parser = Parser(path_to_inventory) + except Exception as exc: + parser = None + return parser, exc + + + def _read_inventory_file(self, path_to_inventory): + exc = None + try: + inv = read_inventory(path_to_inventory) + except Exception as exc: + inv = None + 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. @@ -196,6 +376,33 @@ def read_metadata(path_to_inventory): 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): tr, invtype, inobj, unit, force = input_tuple diff --git a/pylot/core/util/event.py b/pylot/core/util/event.py index 4ded8d50..a4f1e9cb 100644 --- a/pylot/core/util/event.py +++ b/pylot/core/util/event.py @@ -7,6 +7,7 @@ from obspy import UTCDateTime from obspy.core.event import Event as ObsPyEvent from obspy.core.event import Origin, ResourceIdentifier from pylot.core.io.phases import picks_from_picksdict +from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT class Event(ObsPyEvent): @@ -33,6 +34,7 @@ class Event(ObsPyEvent): self._testEvent = False self._refEvent = False self.get_notes() + self.get_obspy_event_info() def get_notes_path(self): """ @@ -43,6 +45,18 @@ class Event(ObsPyEvent): notesfile = os.path.join(self.path, 'notes.txt') return notesfile + def get_obspy_event_info(self): + infile_pickle = os.path.join(self.path, 'info/event.pkl') + if not os.path.isfile(infile_pickle): + return + try: + event_dmt = qml_from_obspyDMT(infile_pickle) + except Exception as e: + print('Could not get obspy event info: {}'.format(e)) + return + self.magnitudes = event_dmt.magnitudes + self.origins = event_dmt.origins + def get_notes(self): """ set self.note attribute to content of notes file diff --git a/pylot/core/util/map_projection.py b/pylot/core/util/map_projection.py deleted file mode 100644 index 66cfb5ce..00000000 --- a/pylot/core/util/map_projection.py +++ /dev/null @@ -1,377 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -import matplotlib.pyplot as plt -import numpy as np -import obspy -from PySide import QtGui -from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar -from mpl_toolkits.basemap import Basemap -from pylot.core.util.widgets import PickDlg -from scipy.interpolate import griddata - -plt.interactive(False) - - -class map_projection(QtGui.QWidget): - def __init__(self, parent, figure=None): - ''' - - :param: picked, can be False, auto, manual - :value: str - ''' - QtGui.QWidget.__init__(self) - self._parent = parent - self.metadata = parent.metadata - self.parser = parent.metadata[1] - self.picks = None - self.picks_dict = None - self.eventLoc = None - self.figure = figure - self.init_graphics() - self.init_stations() - self.init_basemap(resolution='l') - self.init_map() - self._style = parent._style - # self.show() - - def init_map(self): - self.init_lat_lon_dimensions() - self.init_lat_lon_grid() - self.init_x_y_dimensions() - self.connectSignals() - self.draw_everything() - - def onpick(self, event): - ind = event.ind - button = event.mouseevent.button - if ind == [] or not button == 1: - return - data = self._parent.get_data().getWFData() - for index in ind: - station = str(self.station_names[index].split('.')[-1]) - try: - pickDlg = PickDlg(self, parameter=self._parent._inputs, - data=data.select(station=station), - station=station, - picks=self._parent.get_current_event().getPick(station), - autopicks=self._parent.get_current_event().getAutopick(station)) - except Exception as e: - message = 'Could not generate Plot for station {st}.\n {er}'.format(st=station, er=e) - self._warn(message) - print(message, e) - return - pyl_mw = self._parent - try: - if pickDlg.exec_(): - pyl_mw.setDirty(True) - pyl_mw.update_status('picks accepted ({0})'.format(station)) - replot = pyl_mw.get_current_event().setPick(station, pickDlg.getPicks()) - self._refresh_drawings() - if replot: - pyl_mw.plotWaveformData() - pyl_mw.drawPicks() - pyl_mw.draw() - else: - pyl_mw.drawPicks(station) - pyl_mw.draw() - else: - pyl_mw.update_status('picks discarded ({0})'.format(station)) - except Exception as e: - message = 'Could not save picks for station {st}.\n{er}'.format(st=station, er=e) - self._warn(message) - print(message, e) - - def connectSignals(self): - self.comboBox_phase.currentIndexChanged.connect(self._refresh_drawings) - self.zoom_id = self.basemap.ax.figure.canvas.mpl_connect('scroll_event', self.zoom) - - def init_graphics(self): - if not self.figure: - if not hasattr(self._parent, 'am_figure'): - self.figure = plt.figure() - self.toolbar = NavigationToolbar(self.figure.canvas, self) - else: - self.figure = self._parent.am_figure - self.toolbar = self._parent.am_toolbar - - self.main_ax = self.figure.add_subplot(111) - self.canvas = self.figure.canvas - - self.main_box = QtGui.QVBoxLayout() - self.setLayout(self.main_box) - - self.top_row = QtGui.QHBoxLayout() - self.main_box.addLayout(self.top_row) - - self.comboBox_phase = QtGui.QComboBox() - self.comboBox_phase.insertItem(0, 'P') - self.comboBox_phase.insertItem(1, 'S') - - self.comboBox_am = QtGui.QComboBox() - self.comboBox_am.insertItem(0, 'auto') - self.comboBox_am.insertItem(1, 'manual') - - self.top_row.addWidget(QtGui.QLabel('Select a phase: ')) - self.top_row.addWidget(self.comboBox_phase) - self.top_row.setStretch(1, 1) # set stretch of item 1 to 1 - - self.main_box.addWidget(self.canvas) - self.main_box.addWidget(self.toolbar) - - def init_stations(self): - def get_station_names_lat_lon(parser): - station_names = [] - lat = [] - lon = [] - for station in parser.stations: - station_name = station[0].station_call_letters - network = station[0].network_code - if not station_name in station_names: - station_names.append(network + '.' + station_name) - lat.append(station[0].latitude) - lon.append(station[0].longitude) - return station_names, lat, lon - - station_names, lat, lon = get_station_names_lat_lon(self.parser) - self.station_names = station_names - self.lat = lat - self.lon = lon - - def init_picks(self): - phase = self.comboBox_phase.currentText() - - def get_picks(station_names): - picks = [] - for station in station_names: - try: - station = station.split('.')[-1] - picks.append(self.picks_dict[station][phase]['mpp']) - except: - picks.append(np.nan) - return picks - - def get_picks_rel(picks): - picks_rel = [] - picks_utc = [] - for pick in picks: - if type(pick) is obspy.core.utcdatetime.UTCDateTime: - picks_utc.append(pick) - minp = min(picks_utc) - for pick in picks: - if type(pick) is obspy.core.utcdatetime.UTCDateTime: - pick -= minp - picks_rel.append(pick) - return picks_rel - - self.picks = get_picks(self.station_names) - self.picks_rel = get_picks_rel(self.picks) - - def init_picks_active(self): - def remove_nan_picks(picks): - picks_no_nan = [] - for pick in picks: - if not np.isnan(pick): - picks_no_nan.append(pick) - return picks_no_nan - - self.picks_no_nan = remove_nan_picks(self.picks_rel) - - def init_stations_active(self): - def remove_nan_lat_lon(picks, lat, lon): - lat_no_nan = [] - lon_no_nan = [] - for index, pick in enumerate(picks): - if not np.isnan(pick): - lat_no_nan.append(lat[index]) - lon_no_nan.append(lon[index]) - return lat_no_nan, lon_no_nan - - self.lat_no_nan, self.lon_no_nan = remove_nan_lat_lon(self.picks_rel, self.lat, self.lon) - - def init_lat_lon_dimensions(self): - def get_lon_lat_dim(lon, lat): - londim = max(lon) - min(lon) - latdim = max(lat) - min(lat) - return londim, latdim - - self.londim, self.latdim = get_lon_lat_dim(self.lon, self.lat) - - def init_x_y_dimensions(self): - def get_x_y_dim(x, y): - xdim = max(x) - min(x) - ydim = max(y) - min(y) - return xdim, ydim - - self.x, self.y = self.basemap(self.lon, self.lat) - self.xdim, self.ydim = get_x_y_dim(self.x, self.y) - - def init_basemap(self, resolution='l'): - # basemap = Basemap(projection=projection, resolution = resolution, ax=self.main_ax) - basemap = Basemap(projection='lcc', resolution=resolution, ax=self.main_ax, - width=5e6, height=2e6, - lat_0=(min(self.lat) + max(self.lat)) / 2., - lon_0=(min(self.lon) + max(self.lon)) / 2.) - - # basemap.fillcontinents(color=None, lake_color='aqua',zorder=1) - basemap.drawmapboundary(zorder=2) # fill_color='darkblue') - basemap.shadedrelief(zorder=3) - basemap.drawcountries(zorder=4) - basemap.drawstates(zorder=5) - basemap.drawcoastlines(zorder=6) - self.basemap = basemap - self.figure.tight_layout() - - def init_lat_lon_grid(self): - def get_lat_lon_axis(lat, lon): - steplat = (max(lat) - min(lat)) / 250 - steplon = (max(lon) - min(lon)) / 250 - - lataxis = np.arange(min(lat), max(lat), steplat) - lonaxis = np.arange(min(lon), max(lon), steplon) - return lataxis, lonaxis - - def get_lat_lon_grid(lataxis, lonaxis): - longrid, latgrid = np.meshgrid(lonaxis, lataxis) - return latgrid, longrid - - self.lataxis, self.lonaxis = get_lat_lon_axis(self.lat, self.lon) - self.latgrid, self.longrid = get_lat_lon_grid(self.lataxis, self.lonaxis) - - def init_picksgrid(self): - self.picksgrid_no_nan = griddata((self.lat_no_nan, self.lon_no_nan), - self.picks_no_nan, (self.latgrid, self.longrid), - method='linear') ################## - - def draw_contour_filled(self, nlevel='50'): - levels = np.linspace(min(self.picks_no_nan), max(self.picks_no_nan), nlevel) - self.contourf = self.basemap.contourf(self.longrid, self.latgrid, self.picksgrid_no_nan, - levels, latlon=True, zorder=9, alpha=0.5) - - def scatter_all_stations(self): - self.sc = self.basemap.scatter(self.lon, self.lat, s=50, facecolor='none', latlon=True, - zorder=10, picker=True, edgecolor='m', label='Not Picked') - self.cid = self.canvas.mpl_connect('pick_event', self.onpick) - if self.eventLoc: - lat, lon = self.eventLoc - self.sc_event = self.basemap.scatter(lon, lat, s=100, facecolor='red', - latlon=True, zorder=11, label='Event (might be outside map region)') - - def scatter_picked_stations(self): - lon = self.lon_no_nan - lat = self.lat_no_nan - - # workaround because of an issue with latlon transformation of arrays with len <3 - if len(lon) <= 2 and len(lat) <= 2: - self.sc_picked = self.basemap.scatter(lon[0], lat[0], s=50, facecolor='white', - c=self.picks_no_nan[0], latlon=True, zorder=11, label='Picked') - if len(lon) == 2 and len(lat) == 2: - self.sc_picked = self.basemap.scatter(lon[1], lat[1], s=50, facecolor='white', - c=self.picks_no_nan[1], latlon=True, zorder=11) - else: - self.sc_picked = self.basemap.scatter(lon, lat, s=50, facecolor='white', - c=self.picks_no_nan, latlon=True, zorder=11, label='Picked') - - def annotate_ax(self): - self.annotations = [] - for index, name in enumerate(self.station_names): - self.annotations.append(self.main_ax.annotate(' %s' % name, xy=(self.x[index], self.y[index]), - fontsize='x-small', color='white', zorder=12)) - self.legend = self.main_ax.legend(loc=1) - - def add_cbar(self, label): - cbar = self.main_ax.figure.colorbar(self.sc_picked, fraction=0.025) - cbar.set_label(label) - return cbar - - def refresh_drawings(self, picks=None): - self.picks_dict = picks - self._refresh_drawings() - - def _refresh_drawings(self): - self.remove_drawings() - self.draw_everything() - - def draw_everything(self): - if self.picks_dict: - self.init_picks() - self.init_picks_active() - self.init_stations_active() - if len(self.picks_no_nan) >= 3: - self.init_picksgrid() - self.draw_contour_filled() - self.scatter_all_stations() - if self.picks_dict: - self.scatter_picked_stations() - self.cbar = self.add_cbar(label='Time relative to first onset [s]') - self.comboBox_phase.setEnabled(True) - else: - self.comboBox_phase.setEnabled(False) - self.annotate_ax() - self.canvas.draw() - - def remove_drawings(self): - if hasattr(self, 'sc_picked'): - self.sc_picked.remove() - del (self.sc_picked) - if hasattr(self, 'sc_event'): - self.sc_event.remove() - del (self.sc_event) - if hasattr(self, 'cbar'): - self.cbar.remove() - del (self.cbar) - if hasattr(self, 'contourf'): - self.remove_contourf() - del (self.contourf) - if hasattr(self, 'cid'): - self.canvas.mpl_disconnect(self.cid) - del (self.cid) - try: - self.sc.remove() - except Exception as e: - print('Warning: could not remove station scatter plot.\nReason: {}'.format(e)) - try: - self.legend.remove() - except Exception as e: - print('Warning: could not remove legend. Reason: {}'.format(e)) - self.canvas.draw() - - def remove_contourf(self): - for item in self.contourf.collections: - item.remove() - - def remove_annotations(self): - for annotation in self.annotations: - annotation.remove() - - def zoom(self, event): - map = self.basemap - xlim = map.ax.get_xlim() - ylim = map.ax.get_ylim() - x, y = event.xdata, event.ydata - zoom = {'up': 1. / 2., - 'down': 2.} - - if not event.xdata or not event.ydata: - return - - if event.button in zoom: - factor = zoom[event.button] - xdiff = (xlim[1] - xlim[0]) * factor - xl = x - 0.5 * xdiff - xr = x + 0.5 * xdiff - ydiff = (ylim[1] - ylim[0]) * factor - yb = y - 0.5 * ydiff - yt = y + 0.5 * ydiff - - if xl < map.xmin or yb < map.ymin or xr > map.xmax or yt > map.ymax: - xl, xr = map.xmin, map.xmax - yb, yt = map.ymin, map.ymax - map.ax.set_xlim(xl, xr) - map.ax.set_ylim(yb, yt) - map.ax.figure.canvas.draw() - - def _warn(self, message): - self.qmb = QtGui.QMessageBox(QtGui.QMessageBox.Icon.Warning, - 'Warning', message) - self.qmb.show() diff --git a/pylot/core/util/obspyDMT_interface.py b/pylot/core/util/obspyDMT_interface.py index 661e4f5b..47d5c687 100644 --- a/pylot/core/util/obspyDMT_interface.py +++ b/pylot/core/util/obspyDMT_interface.py @@ -25,4 +25,20 @@ def check_obspydmt_eventfolder(folder): except Exception as e: return False, e -check_obspydmt_eventfolder('20110311_054623.a') \ No newline at end of file +def qml_from_obspyDMT(path): + import pickle + from obspy.core.event import Event, Magnitude, Origin + + if not os.path.exists(path): + return IOError('Could not find Event at {}'.format(path)) + infile = open(path, 'rb') + event_dmt = pickle.load(infile) + ev = Event(resource_id=event_dmt['event_id']) + origin = Origin(resource_id=event_dmt['origin_id'], time=event_dmt['datetime'], longitude=event_dmt['longitude'], + latitude=event_dmt['latitude'], depth=event_dmt['depth']) + mag = Magnitude(mag=event_dmt['magnitude'], magnitude_type=event_dmt['magnitude_type'], + origin_id=event_dmt['origin_id']) + ev.magnitudes.append(mag) + ev.origins.append(origin) + return ev + diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index cf03abac..ed556a9b 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -16,6 +16,9 @@ import time import numpy as np +import matplotlib +matplotlib.use('QT4Agg') + from matplotlib.figure import Figure try: @@ -26,6 +29,7 @@ from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT from matplotlib.widgets import MultiCursor from matplotlib.tight_layout import get_renderer, get_subplotspec_list, get_tight_layout_figure from scipy.signal import argrelmin, argrelmax +from obspy import read from PySide import QtCore, QtGui from PySide.QtGui import QAction, QApplication, QCheckBox, QComboBox, \ @@ -447,23 +451,25 @@ class WaveformWidgetPG(QtGui.QWidget): self.orig_parent = parent # attribute plotdict is a dictionary connecting position and a name self.plotdict = dict() + # init labels + self.xlabel = None + self.ylabel = None + self.title = None # create plot self.main_layout = QtGui.QVBoxLayout() self.label_layout = QtGui.QHBoxLayout() - self.status_label = QtGui.QLabel() - self.perm_label = QtGui.QLabel() - self.setLayout(self.main_layout) + self.add_labels() + self.connect_signals() self.plotWidget = self.pg.PlotWidget(self.parent(), title=title) self.main_layout.addWidget(self.plotWidget) self.main_layout.addLayout(self.label_layout) - self.label_layout.addWidget(self.status_label) - self.label_layout.addWidget(self.perm_label) + self.init_labels() + self.activateObspyDMToptions(False) self.plotWidget.showGrid(x=False, y=True, alpha=0.3) - self.plotWidget.hideAxis('bottom') - self.plotWidget.hideAxis('left') self.wfstart, self.wfend = 0, 0 self.pen_multicursor = self.pg.mkPen(self.parent()._style['multicursor']['rgba']) self.pen_linecolor = self.pg.mkPen(self.parent()._style['linecolor']['rgba']) + self.pen_linecolor_highlight = self.pg.mkPen((255, 100, 100, 255)) self.pen_linecolor_syn = self.pg.mkPen((100, 0, 255, 255)) self.reinitMoveProxy() self._proxy = self.pg.SignalProxy(self.plotWidget.scene().sigMouseMoved, rateLimit=60, slot=self.mouseMoved) @@ -489,12 +495,59 @@ class WaveformWidgetPG(QtGui.QWidget): self.vLine.setPos(mousePoint.x()) self.hLine.setPos(mousePoint.y()) + def connect_signals(self): + self.qcombo_processed.activated.connect(self.parent().newWF) + self.syn_checkbox.clicked.connect(self.parent().newWF) + + def init_labels(self): + self.label_layout.addWidget(self.status_label) + for label in self.perm_labels: + self.label_layout.addWidget(label) + self.label_layout.addWidget(self.syn_checkbox) + self.label_layout.addWidget(self.qcombo_processed) + self.syn_checkbox.setLayoutDirection(Qt.RightToLeft) + self.label_layout.setStretch(0, 4) + self.label_layout.setStretch(1, 0) + self.label_layout.setStretch(2, 0) + self.label_layout.setStretch(3, 0) + self.label_layout.setStretch(4, 3) + self.label_layout.setStretch(5, 1) + + def add_labels(self): + self.status_label = QtGui.QLabel() + self.perm_labels = [] + for index in range(3): + label = QtGui.QLabel() + self.perm_labels.append(label) + self.qcombo_processed = QtGui.QComboBox() + self.syn_checkbox = QtGui.QCheckBox('synthetics') + self.addQCboxItem('processed', 'green') + self.addQCboxItem('raw', 'black') + #self.perm_qcbox_right.setAlignment(2) + self.setLayout(self.main_layout) + def getPlotDict(self): return self.plotdict - def setPermText(self, text=None, color='black'): - self.perm_label.setText(text) - self.perm_label.setStyleSheet('color: {}'.format(color)) + def activateObspyDMToptions(self, activate): + self.syn_checkbox.setVisible(activate) + self.qcombo_processed.setVisible(activate) + + def setPermText(self, number, text=None, color='black'): + if not 0 <= number < len(self.perm_labels): + raise ValueError('No label for number {}'.format(number)) + self.perm_labels[number].setText(text) + self.perm_labels[number].setStyleSheet('color: {}'.format(color)) + + def addQCboxItem(self, text=None, color='black'): + item = QtGui.QStandardItem(text) + model = self.qcombo_processed.model() + model.appendRow(item) + item.setForeground(QtGui.QColor('{}'.format(color))) + + def setQCboxItem(self, text): + index = self.qcombo_processed.findText(text) + self.qcombo_processed.setCurrentIndex(index) def setPlotDict(self, key, value): self.plotdict[key] = value @@ -529,6 +582,14 @@ class WaveformWidgetPG(QtGui.QWidget): else: st_select = wfdata + gaps = st_select.get_gaps() + if gaps: + merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps] + st_select.merge() + print('Merged the following stations because of gaps:') + for merged_station in merged: + print(merged_station) + # list containing tuples of network, station, channel (for sorting) nsc = [] for trace in st_select: @@ -552,6 +613,8 @@ class WaveformWidgetPG(QtGui.QWidget): st_syn = wfsyn.select(network=network, station=station, channel=channel) if st_syn: trace_syn = st_syn[0].copy() + else: + trace_syn = Trace() if mapping: comp = channel[-1] n = compclass.getPlotPosition(str(comp)) @@ -568,9 +631,11 @@ class WaveformWidgetPG(QtGui.QWidget): time_ax_syn = prepTimeAxis(stime_syn, trace_syn) if method == 'fast': - trace.data, time_ax = self.minMax(trace, time_ax) + trace.data, time_ax = self.minMax(trace, time_ax) + if trace_syn: + trace_syn.data, time_ax_syn = self.minMax(trace_syn, time_ax_syn) - if time_ax not in [None, []]: + if len(time_ax) > 0: if not scaleddata: trace.detrend('constant') trace.normalize(np.max(np.abs(trace.data)) * 2) @@ -581,23 +646,23 @@ class WaveformWidgetPG(QtGui.QWidget): times = np.array([time for index, time in enumerate(time_ax) if not index % nth_sample]) times_syn = np.array([time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else []) trace.data = np.array([datum + n for index, datum in enumerate(trace.data) if not index % nth_sample]) - trace.data_syn = np.array([datum + n for index, datum in enumerate(trace.data_syn) + trace_syn.data = np.array([datum + n 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.data_syn)) + times_syn, trace_syn.data)) self.setPlotDict(n, (station, channel, network)) self.xlabel = 'seconds since {0}'.format(self.wfstart) self.ylabel = '' self.setXLims([0, self.wfend - self.wfstart]) self.setYLims([0.5, nmax + 0.5]) - return plots + return plots, gaps def minMax(self, trace, time_ax): ''' create min/max array for fast plotting (approach based on obspy __plot_min_max function) :returns data, time_ax ''' - npixel = self.width() + npixel = self.orig_parent.width() ndata = len(trace.data) pts_per_pixel = ndata/npixel if pts_per_pixel < 2: @@ -1013,6 +1078,14 @@ class PylotCanvas(FigureCanvas): if mapping: plot_positions = self.calcPlotPositions(st_select, compclass) + gaps = st_select.get_gaps() + if gaps: + merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps] + st_select.merge() + print('Merged the following stations because of gaps:') + for merged_station in merged: + print(merged_station) + # list containing tuples of network, station, channel and plot position (for sorting) nsc = [] for plot_pos, trace in enumerate(st_select): @@ -1245,9 +1318,10 @@ class PickDlg(QDialog): def __init__(self, parent=None, data=None, station=None, network=None, picks=None, autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None, - event=None, filteroptions=None, model='iasp91'): + event=None, filteroptions=None, model='iasp91', wftype=None): super(PickDlg, self).__init__(parent, 1) self.orig_parent = parent + self.setAttribute(Qt.WA_DeleteOnClose) # initialize attributes self.parameter = parameter @@ -1256,6 +1330,7 @@ class PickDlg(QDialog): self.network = network self.rotate = rotate self.metadata = metadata + self.wftype = wftype self.pylot_event = event self.components = 'ZNE' self.currentPhase = None @@ -1305,7 +1380,7 @@ class PickDlg(QDialog): self.cur_ylim = None # set attribute holding data - if data is None: + if data is None or not data: try: data = parent.get_data().getWFData().copy() self.data = data.select(station=station) @@ -1328,7 +1403,7 @@ class PickDlg(QDialog): # fill compare and scale channels self.compareChannel.addItem('-', None) - self.scaleChannel.addItem('normalized', None) + self.scaleChannel.addItem('individual', None) for trace in self.getWFData(): channel = trace.stats.channel @@ -1345,8 +1420,12 @@ class PickDlg(QDialog): actionS.setChecked(self.getChannelSettingsS(channel)) # plot data + title = self.getStation() + if self.wftype is not None: + title += ' | ({})'.format(self.wftype) + self.multicompfig.plotWFData(wfdata=self.getWFData(), - title=self.getStation()) + title=title) self.multicompfig.setZoomBorders2content() @@ -1380,7 +1459,6 @@ class PickDlg(QDialog): self.setWindowTitle('Pickwindow on station: {}'.format(self.getStation())) self.setWindowState(QtCore.Qt.WindowMaximized) - self.deleteLater() def setupUi(self): menuBar = QtGui.QMenuBar(self) @@ -1526,7 +1604,7 @@ class PickDlg(QDialog): _dialtoolbar.addWidget(QtGui.QLabel('Compare to channel: ')) _dialtoolbar.addWidget(self.compareChannel) _dialtoolbar.addSeparator() - _dialtoolbar.addWidget(QtGui.QLabel('Scale by: ')) + _dialtoolbar.addWidget(QtGui.QLabel('Scaling: ')) _dialtoolbar.addWidget(self.scaleChannel) # layout the innermost widget @@ -1649,7 +1727,7 @@ class PickDlg(QDialog): self.arrivalsText.append(ax.text(time_rel, ylims[0], arrival.name, color='0.5')) def drawArrivalsText(self): - return self.drawArrivals(True) + return self.drawArrivals(textOnly=True) def refreshArrivalsText(self, event=None): self.removeArrivalsText() @@ -2514,6 +2592,9 @@ class PickDlg(QDialog): 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() scale_channel = self.scaleChannel.currentText() @@ -2601,6 +2682,7 @@ class PickDlg(QDialog): phase = 'S' filter = True self.plotWFData(phase=phase, filter=filter) + self.drawArrivals() def resetZoom(self): ax = self.multicompfig.axes[0] @@ -2864,7 +2946,7 @@ class AutoPickWidget(MultiEventWidget): def reinitEvents2plot(self): for eventID, eventDict in self.events2plot.items(): for widget_key, widget in eventDict.items(): - widget.setParent(None) + del(widget) self.events2plot = {} self.eventbox.clear() self.refresh_plot_tabs() @@ -2946,13 +3028,15 @@ class TuneAutopicker(QWidget): :type: PyLoT Mainwindow ''' - def __init__(self, parent): + def __init__(self, parent, obspy_dmt=False): QtGui.QWidget.__init__(self, parent, 1) self._style = parent._style self.setWindowTitle('PyLoT - Tune Autopicker') self.parameter = self.parent()._inputs self.fig_dict = self.parent().fig_dict self.data = Data() + self.obspy_dmt = obspy_dmt + self.wftype = None self.pdlg_widget = None self.pylot_picks = None self.init_main_layouts() @@ -3009,8 +3093,45 @@ class TuneAutopicker(QWidget): self.eventBox.activated.connect(self.fill_tabs) self.stationBox.activated.connect(self.fill_tabs) + def catch_station_ids(self): + self.station_ids = {} + eventpath = self.get_current_event_fp() + self.wftype = 'processed' if self.obspy_dmt else '' + wf_path = os.path.join(eventpath, self.wftype) + if not os.path.exists(wf_path) and self.obspy_dmt: + self.wftype = 'raw' + wf_path = os.path.join(eventpath, self.wftype) + for filename in os.listdir(wf_path): + filename = os.path.join(eventpath, self.wftype, filename) + try: + st = read(filename, headonly=True) + except Exception as e: + print('Warning: Could not read file {} as a stream object: {}'.format(filename, e)) + continue + for trace in st: + network = trace.stats.network + station = trace.stats.station + location = trace.stats.location + station_id = '{}.{}.{}'.format(network, station, location) + if not station_id in self.station_ids: + self.station_ids[station_id] = [] + self.station_ids[station_id].append(filename) + def fill_stationbox(self): - fnames = self.parent().getWFFnames_from_eventbox(eventbox=self.eventBox) + self.stationBox.clear() + model = self.stationBox.model() + + self.catch_station_ids() + st_ids_list = list(self.station_ids.keys()) + st_ids_list.sort() + for station_id in st_ids_list: + item = QtGui.QStandardItem(station_id) + if station_id.split('.')[1] in self.get_current_event().pylot_picks: + item.setBackground(self.parent()._ref_test_colors['ref']) + model.appendRow(item) + + def load_wf_data(self): + fnames = self.station_ids[self.get_current_station_id()] self.data.setWFData(fnames) wfdat = self.data.getWFData() # all available streams # remove possible underscores in station names @@ -3022,21 +3143,6 @@ class TuneAutopicker(QWidget): wfdat = check4rotated(wfdat, self.parent().metadata, verbosity=0) # trim station components to same start value trim_station_components(wfdat, trim_start=True, trim_end=False) - self.stationBox.clear() - stations = [] - for trace in self.data.getWFData(): - station = trace.stats.station - network = trace.stats.network - ns_tup = (str(network), str(station)) - if not ns_tup in stations: - stations.append(ns_tup) - stations.sort() - model = self.stationBox.model() - for network, station in stations: - item = QtGui.QStandardItem(network + '.' + station) - if station in self.get_current_event().pylot_picks: - item.setBackground(self.parent()._ref_test_colors['ref']) - model.appendRow(item) def init_figure_tabs(self): self.figure_tabs = QtGui.QTabWidget() @@ -3093,10 +3199,14 @@ class TuneAutopicker(QWidget): def get_current_event_autopicks(self, station): event = self.get_current_event() if event.pylot_autopicks: - return event.pylot_autopicks[station] + if station in event.pylot_autopicks: + return event.pylot_autopicks[station] def get_current_station(self): - return str(self.stationBox.currentText()).split('.')[-1] + return str(self.stationBox.currentText()).split('.')[1] + + def get_current_station_id(self): + return str(self.stationBox.currentText()) def gen_tab_widget(self, name, canvas): widget = QtGui.QWidget() @@ -3111,17 +3221,19 @@ class TuneAutopicker(QWidget): self.pdlg_widget.setParent(None) self.pdlg_widget = None return + self.load_wf_data() station = self.get_current_station() - data = self.data.getWFData() + wfdata = self.data.getWFData() metadata = self.parent().metadata event = self.get_current_event() filteroptions = self.parent().filteroptions - self.pickDlg = PickDlg(self.parent(), data=data.select(station=station), + wftype = self.wftype if self.obspy_dmt else '' + self.pickDlg = PickDlg(self.parent(), data=wfdata.select(station=station).copy(), station=station, parameter=self.parameter, picks=self.get_current_event_picks(station), autopicks=self.get_current_event_autopicks(station), metadata=metadata, event=event, filteroptions=filteroptions, - embedded=True) + embedded=True, wftype=wftype) self.pickDlg.update_picks.connect(self.picks_from_pickdlg) self.pickDlg.update_picks.connect(self.fill_eventbox) self.pickDlg.update_picks.connect(self.fill_stationbox) @@ -3301,14 +3413,17 @@ class TuneAutopicker(QWidget): if not station: self._warn('No station selected') return + wfpath = self.wftype if self.obspy_dmt else '' args = {'parameter': self.parameter, 'station': station, 'fnames': 'None', 'eventid': [self.get_current_event_fp()], 'iplot': 2, 'fig_dict': self.fig_dict, - 'locflag': 0, - 'savexml': False} + 'savexml': False, + 'obspyDMT_wfpath': wfpath} + event = self.get_current_event() + self.parent().saveData(event, event.path, '.xml') for key in self.fig_dict.keys(): if not key == 'plot_style': self.fig_dict[key].clear() @@ -3364,8 +3479,7 @@ class TuneAutopicker(QWidget): if hasattr(self, 'pdlg_widget'): if self.pdlg_widget: self.pdlg_widget.setParent(None) - # TODO: removing widget by parent deletion raises exception when activating stationbox: - # RuntimeError: Internal C++ object (PylotCanvas) already deleted. + del(self.pdlg_widget) if hasattr(self, 'overview'): self.overview.setParent(None) if hasattr(self, 'p_tabs'): diff --git a/pylot/styles/bright.qss b/pylot/styles/bright.qss index 71e5c5a3..5f3de6a7 100644 --- a/pylot/styles/bright.qss +++ b/pylot/styles/bright.qss @@ -179,9 +179,9 @@ background-color:transparent; QTabWidget::pane{ background-color:rgba(0, 0, 0, 255); -border-style:solid; -border-color:rgba(245, 245, 245, 255); -border-width:1px; +margin: 0px, 0px, 0px, 0px; +padding: 0px; +border-width:0px; } QTabWidget::tab{ diff --git a/pylot/styles/dark.qss b/pylot/styles/dark.qss index fc866e20..0ca4cd46 100644 --- a/pylot/styles/dark.qss +++ b/pylot/styles/dark.qss @@ -178,9 +178,9 @@ background-color:transparent; QTabWidget::pane{ background-color:rgba(70, 70, 80, 255); -border-style:solid; -border-color:rgba(70, 70, 80, 255); -border-width:1px; +margin: 0px, 0px, 0px, 0px; +padding: 0px; +border-width:0px; } QTabWidget::tab{