WIP: Simplify data structure #39

Draft
sebastianw wants to merge 31 commits from 38-simplify-data-structure into develop
7 changed files with 352 additions and 399 deletions
Showing only changes of commit 76d4ec290c - Show all commits

View File

@ -1773,7 +1773,7 @@ class MainWindow(QMainWindow):
def getStime(self):
if self.get_data():
return full_range(self.get_data().getWFData())[0]
return full_range(self.get_data().get_wf_data())[0]
def addActions(self, target, actions):
for action in actions:
@ -1983,7 +1983,7 @@ class MainWindow(QMainWindow):
tstart = None
tstop = None
self.data.setWFData(self.fnames,
self.data.set_wf_data(self.fnames,
self.fnames_comp,
checkRotated=True,
metadata=self.metadata,
@ -2035,7 +2035,7 @@ class MainWindow(QMainWindow):
def get_npts_to_plot(self):
if not hasattr(self.data, 'wfdata'):
return 0
return sum(trace.stats.npts for trace in self.data.getWFData())
return sum(trace.stats.npts for trace in self.data.get_wf_data())
def connectWFplotEvents(self):
'''
@ -2248,14 +2248,14 @@ class MainWindow(QMainWindow):
zne_text = {'Z': 'vertical', 'N': 'north-south', 'E': 'east-west'}
comp = self.getComponent()
title = 'section: {0} components'.format(zne_text[comp])
wfst = self.get_data().getWFData()
wfst = self.get_data().get_wf_data()
wfsyn = self.get_data().getAltWFdata()
if self.filterActionP.isChecked() and filter:
self.filterWaveformData(plot=False, phase='P')
elif self.filterActionS.isChecked() and filter:
self.filterWaveformData(plot=False, phase='S')
# wfst = self.get_data().getWFData().select(component=comp)
# wfst += self.get_data().getWFData().select(component=alter_comp)
# wfst = self.get_data().get_wf_data().select(component=comp)
# wfst += self.get_data().get_wf_data().select(component=alter_comp)
plotWidget = self.getPlotWidget()
self.adjustPlotHeight()
if get_bool(settings.value('large_dataset')) == True:
@ -2270,7 +2270,7 @@ class MainWindow(QMainWindow):
def adjustPlotHeight(self):
if self.pg:
return
height_need = len(self.data.getWFData()) * self.height_factor
height_need = len(self.data.get_wf_data()) * self.height_factor
plotWidget = self.getPlotWidget()
if self.tabs.widget(0).frameSize().height() < height_need:
plotWidget.setMinimumHeight(height_need)
@ -2290,24 +2290,24 @@ class MainWindow(QMainWindow):
self.plotWaveformDataThread()
def pushFilterWF(self, param_args):
self.get_data().filterWFData(param_args)
self.get_data().filter_wf_data(param_args)
def filterP(self):
self.filterActionS.setChecked(False)
if self.filterActionP.isChecked():
self.filterWaveformData(phase='P')
else:
self.resetWFData()
self.reset_wf_data()
def filterS(self):
self.filterActionP.setChecked(False)
if self.filterActionS.isChecked():
self.filterWaveformData(phase='S')
else:
self.resetWFData()
self.reset_wf_data()
def resetWFData(self):
self.get_data().resetWFData()
def reset_wf_data(self):
self.get_data().reset_wf_data()
self.plotWaveformDataThread()
def filterWaveformData(self, plot=True, phase=None):
@ -2326,11 +2326,11 @@ class MainWindow(QMainWindow):
kwargs = self.getFilterOptions()[phase].parseFilterOptions()
self.pushFilterWF(kwargs)
else:
self.get_data().resetWFData()
self.get_data().reset_wf_data()
elif self.filterActionP.isChecked() or self.filterActionS.isChecked():
self.adjustFilterOptions()
else:
self.get_data().resetWFData()
self.get_data().reset_wf_data()
if plot:
self.plotWaveformDataThread(filter=False)
@ -2531,10 +2531,10 @@ class MainWindow(QMainWindow):
show_comp_data=self.dataPlot.comp_checkbox.isChecked())
if self.filterActionP.isChecked():
pickDlg.currentPhase = "P"
pickDlg.filterWFData()
pickDlg.filter_wf_data()
elif self.filterActionS.isChecked():
pickDlg.currentPhase = "S"
pickDlg.filterWFData()
pickDlg.filter_wf_data()
pickDlg.nextStation.setChecked(self.nextStation)
pickDlg.nextStation.stateChanged.connect(self.toggle_next_station)
if pickDlg.exec_():
@ -3478,7 +3478,7 @@ class MainWindow(QMainWindow):
if not self.metadata:
return None
wf_copy = self.get_data().getWFData().copy()
wf_copy = self.get_data().get_wf_data().copy()
wf_select = Stream()
# restitute only picked traces

View File

@ -243,7 +243,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
pylot_event = Event(eventpath) # event should be path to event directory
data.setEvtData(pylot_event)
if fnames == 'None':
data.setWFData(glob.glob(os.path.join(datapath, event_datapath, '*')))
data.set_wf_data(glob.glob(os.path.join(datapath, event_datapath, '*')))
# the following is necessary because within
# multiple event processing no event ID is provided
# in autopylot.in
@ -258,7 +258,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
now.minute)
parameter.setParam(eventID=eventID)
else:
data.setWFData(fnames)
data.set_wf_data(fnames)
eventpath = events[0]
# now = datetime.datetime.now()
@ -268,7 +268,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
# now.hour,
# now.minute)
parameter.setParam(eventID=eventid)
wfdat = data.getWFData() # all available streams
wfdat = data.get_wf_data() # all available streams
if not station == 'all':
wfdat = wfdat.select(station=station)
if not wfdat:

View File

@ -9,11 +9,13 @@ from obspy import UTCDateTime
from pylot.core.io.event import EventData
from pylot.core.io.waveformdata import WaveformData
from pylot.core.util.dataprocessing import Metadata
@dataclass
class Data:
event_data: EventData = field(default_factory=EventData)
waveform_data: WaveformData = field(default_factory=WaveformData)
metadata: Metadata = field(default_factory=Metadata)
_parent: Union[None, 'QtWidgets.QWidget'] = None
def __init__(self, parent=None, evtdata=None):
@ -52,10 +54,17 @@ class Data:
self.waveform_data.dirty = True
def set_wf_data(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0):
return self.waveform_data.set_wf_data(fnames, fnames_alt, check_rotated, metadata, tstart, tstop)
return self.waveform_data.load_waveforms(fnames, fnames_alt, check_rotated, metadata, tstart, tstop)
def reset_wf_data(self):
self.waveform_data.reset_wf_data()
self.waveform_data.reset()
def get_wf_data(self):
return self.waveform_data.wfdata
def rotate_wf_data(self):
self.waveform_data.rotate_zne(self.metadata)
class GenericDataStructure(object):
"""

13
pylot/core/io/utils.py Normal file
View File

@ -0,0 +1,13 @@
import os
from typing import List
def validate_filenames(filenames: List[str]) -> List[str]:
"""
validate a list of filenames for file abundance
:param filenames: list of possible filenames
:type filenames: List[str]
:return: list of valid filenames
:rtype: List[str]
"""
return [fn for fn in filenames if os.path.isfile(fn)]

View File

@ -1,14 +1,13 @@
import logging
import os
from dataclasses import dataclass, field
from typing import Union, List
import numpy as np
from obspy import Stream, read
from obspy.io.sac import SacIOError
from obspy.signal.rotate import rotate2zne
from pylot.core.util.utils import full_range, get_stations
from pylot.core.io.utils import validate_filenames
from pylot.core.util.dataprocessing import Metadata
from pylot.core.util.utils import get_stations, check_for_nan, check4rotated
@dataclass
@ -18,26 +17,39 @@ class WaveformData:
wf_alt: Stream = field(default_factory=Stream)
dirty: bool = False
def set_wf_data(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0):
self.clear_data()
fnames = self.check_fname_exists(fnames)
fnames_alt = self.check_fname_exists(fnames_alt)
def load_waveforms(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0):
fn_list = validate_filenames(fnames)
if not fn_list:
logging.warning('No valid filenames given for loading waveforms')
else:
self.clear()
self.add_waveforms(fn_list)
if fnames:
self.append_wf_data(fnames)
if fnames_alt:
self.append_wf_data(fnames_alt, alternative=True)
self.wfdata, _ = self.check_for_gaps_and_merge(self.wfdata)
self.check_for_nan(self.wfdata)
if check_rotated and metadata:
self.wfdata = self.check4rotated(self.wfdata, metadata, verbosity=0)
self.trim_station_components(self.wfdata, trim_start=True, trim_end=False)
if fnames_alt is None:
pass
else:
alt_fn_list = validate_filenames(fnames_alt)
if not alt_fn_list:
logging.warning('No valid alternative filenames given for loading waveforms')
else:
self.add_waveforms(alt_fn_list, alternative=True)
if not fn_list and not alt_fn_list:
logging.error('No filenames or alternative filenames given for loading waveforms')
return False
self.merge()
self.replace_nan()
if not check_rotated or not metadata:
pass
else:
self.rotate_zne()
self.trim_station_traces()
self.wforiginal = self.wfdata.copy()
self.dirty = False
return True
return False
def append_wf_data(self, fnames: List[str], alternative: bool = False):
def add_waveforms(self, fnames: List[str], alternative: bool = False):
data_stream = self.wf_alt if alternative else self.wfdata
warnmsg = ''
for fname in set(fnames):
@ -55,189 +67,57 @@ class WaveformData:
warnmsg += f'{fname}\n{se}\n'
if warnmsg:
print(f'WARNING in appendWFData: unable to read waveform data\n{warnmsg}')
print(f'WARNING in add_waveforms: unable to read waveform data\n{warnmsg}')
def clear_data(self):
def clear(self):
self.wfdata = Stream()
self.wforiginal = None
self.wf_alt = Stream()
def reset_wf_data(self):
def reset(self):
"""
Resets the waveform data to its original state.
"""
if self.wforiginal:
self.wfdata = self.wforiginal.copy()
else:
self.wfdata = Stream()
self.dirty = False
def check_fname_exists(self, filenames: List[str]) -> List[str]:
return [fn for fn in filenames if os.path.isfile(fn)]
def check_for_gaps_and_merge(self, stream):
def merge(self):
"""
check for gaps in Stream and merge if gaps are found
:param stream: stream of seismic data
:type stream: `~obspy.core.stream.Stream`
:return: data stream, gaps returned from obspy get_gaps
:rtype: `~obspy.core.stream.Stream`
"""
gaps = stream.get_gaps()
gaps = self.wfdata.get_gaps()
if gaps:
merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
stream.merge(method=1)
print('Merged the following stations because of gaps:')
for merged_station in merged:
print(merged_station)
self.wfdata.merge(method=1)
logging.info('Merged the following stations because of gaps:')
for station in merged:
logging.info(station)
return stream, gaps
def check_for_nan(self, stream):
def replace_nan(self):
"""
Replace all NaNs in data with nan_value (in place)
:param stream: stream of seismic data
:type stream: `~obspy.core.stream.Stream`
:param nan_value: value which all NaNs are set to
:type nan_value: float, int
:return: None
Replace all NaNs in data with 0. (in place)
"""
if not stream:
return
for trace in stream:
np.nan_to_num(trace.data, copy=False, nan=0.)
self.wfdata = check_for_nan(self.wfdata)
def check4rotated(self, stream, metadata=None, verbosity=1):
def rotate_zne(self, metadata: Metadata = None):
"""
Check all traces in data. If a trace is not in ZNE rotation (last symbol of channel code is numeric) and the trace
Check all traces in stream for rotation. If a trace is not in ZNE rotation (last symbol of channel code is numeric) and the trace
is in the metadata with azimuth and dip, rotate it to classical ZNE orientation.
Rotating the traces requires them to be of the same length, so, all traces will be trimmed to a common length as a
side effect.
:param stream: stream object containing seismic traces
:type stream: `~obspy.core.stream.Stream`
:param metadata: tuple containing metadata type string and metadata parser object
:type metadata: (str, `~obspy.io.xseed.parser.Parser`)
:param verbosity: if 0 print no information at runtime
:type verbosity: int
:return: stream object with traditionally oriented traces (ZNE) for stations that had misaligned traces (123) before
:rtype: `~obspy.core.stream.Stream`
"""
def rotation_required(trace_ids):
"""
Derive if any rotation is required from the orientation code of the input.
self.wfdata = check4rotated(self.wfdata, metadata)
:param trace_ids: string identifier of waveform data trace
:type trace_ids: List(str)
:return: boolean representing if rotation is necessary for any of the traces
:rtype: bool
def trim_station_traces(self):
"""
orientations = [trace_id[-1] for trace_id in trace_ids]
return any([orientation.isnumeric() for orientation in orientations])
def rotate_components(wfs_in, metadata=None):
"""
Rotate components if orientation code is numeric (= non traditional orientation).
Azimut and dip are fetched from metadata. To be rotated, traces of a station have to be cut to the same length.
Returns unrotated traces of no metadata is provided
:param wfs_in: stream containing seismic traces of a station
:type wfs_in: `~obspy.core.stream.Stream`
:param metadata: tuple containing metadata type string and metadata parser object
:type metadata: (str, `~obspy.io.xseed.parser.Parser`)
:return: stream object with traditionally oriented traces (ZNE)
:rtype: `~obspy.core.stream.Stream`
trim data stream to common time window
"""
if len(wfs_in) < 3:
print(f"Stream {wfs_in=}, has not enough components to rotate.")
return wfs_in
# check if any traces in this station need to be rotated
trace_ids = [trace.id for trace in wfs_in]
if not rotation_required(trace_ids):
logging.debug(f"Stream does not need any rotation: Traces are {trace_ids=}")
return wfs_in
# check metadata quality
t_start = full_range(wfs_in)
try:
azimuths = []
dips = []
for tr_id in trace_ids:
azimuths.append(metadata.get_coordinates(tr_id, t_start)['azimuth'])
dips.append(metadata.get_coordinates(tr_id, t_start)['dip'])
except (KeyError, TypeError) as err:
logging.error(
f"{type(err)=} occurred: {err=} Rotating not possible, not all azimuth and dip information "
f"available in metadata. Stream remains unchanged.")
return wfs_in
except Exception as err:
print(f"Unexpected {err=}, {type(err)=}")
raise
# to rotate all traces must have same length, so trim them
wfs_out = self.trim_station_components(wfs_in, trim_start=True, trim_end=True)
try:
z, n, e = rotate2zne(wfs_out[0], azimuths[0], dips[0],
wfs_out[1], azimuths[1], dips[1],
wfs_out[2], azimuths[2], dips[2])
print('check4rotated: rotated trace {} to ZNE'.format(trace_ids))
# replace old data with rotated data, change the channel code to ZNE
z_index = dips.index(min(
dips)) # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90
# being vertical)
wfs_out[z_index].data = z
wfs_out[z_index].stats.channel = wfs_out[z_index].stats.channel[0:-1] + 'Z'
del trace_ids[z_index]
for trace_id in trace_ids:
coordinates = metadata.get_coordinates(trace_id, t_start)
dip, az = coordinates['dip'], coordinates['azimuth']
trace = wfs_out.select(id=trace_id)[0]
if az > 315 or az <= 45 or 135 < az <= 225:
trace.data = n
trace.stats.channel = trace.stats.channel[0:-1] + 'N'
elif 45 < az <= 135 or 225 < az <= 315:
trace.data = e
trace.stats.channel = trace.stats.channel[0:-1] + 'E'
except ValueError as err:
print(f"{err=} Rotation failed. Stream remains unchanged.")
return wfs_in
return wfs_out
if metadata is None:
if verbosity:
msg = 'Warning: could not rotate traces since no metadata was given\nset Inventory file!'
print(msg)
return stream
stations = get_stations(stream)
for station in stations: # loop through all stations and rotate data if neccessary
wf_station = stream.select(station=station)
rotate_components(wf_station, metadata)
return stream
def trim_station_components(stream, trim_start=True, trim_end=True):
"""
cut a stream so only the part common to all three traces is kept to avoid dealing with offsets
:param stream: stream of seismic data
:type stream: `~obspy.core.stream.Stream`
:param trim_start: trim start of stream
:type trim_start: bool
:param trim_end: trim end of stream
:type trim_end: bool
:return: data stream
:rtype: `~obspy.core.stream.Stream`
"""
starttime = {False: None}
endtime = {False: None}
stations = get_stations(stream)
print('trim_station_components: Will trim stream for trim_start: {} and for '
'trim_end: {}.'.format(trim_start, trim_end))
for station in stations:
wf_station = stream.select(station=station)
starttime[True] = max([trace.stats.starttime for trace in wf_station])
endtime[True] = min([trace.stats.endtime for trace in wf_station])
wf_station.trim(starttime=starttime[trim_start], endtime=endtime[trim_end])
return stream
for station in get_stations(self.wfdata):
station_traces = self.wfdata.select(station=station)
station_traces.trim(starttime=max([trace.stats.starttime for trace in station_traces]),
endtime=min([trace.stats.endtime for trace in station_traces]))

View File

@ -474,8 +474,8 @@ class Array_map(QtWidgets.QWidget):
transform=ccrs.PlateCarree(), label='deleted'))
def openPickDlg(self, ind):
wfdata = self._parent.get_data().getWFData()
wfdata_comp = self._parent.get_data().getWFDataComp()
wfdata = self._parent.get_data().get_wf_data()
wfdata_comp = self._parent.get_data().get_wf_dataComp()
for index in ind:
network, station = self._station_onpick_ids[index].split('.')[:2]
pyl_mw = self._parent

View File

@ -140,18 +140,6 @@ class LogWidget(QtWidgets.QWidget):
self.stderr.append(60 * '#' + '\n\n')
def getDataType(parent):
type = QInputDialog().getItem(parent, "Select phases type", "Type:",
["manual", "automatic"])
if type[0].startswith('auto'):
type = 'auto'
else:
type = type[0]
return type
def plot_pdf(_axes, x, y, annotation, bbox_props, xlabel=None, ylabel=None,
title=None):
# try method or data
@ -795,7 +783,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
def connect_signals(self):
self.qcombo_processed.activated.connect(self.parent().newWF)
self.syn_checkbox.clicked.connect(self.parent().newWF)
self.comp_checkbox.clicked.connect(self.parent().newWF)
def init_labels(self):
self.label_layout.addWidget(self.status_label)
@ -806,13 +794,13 @@ class WaveformWidgetPG(QtWidgets.QWidget):
# use widgets as placeholder, so that child widgets keep position when others are hidden
mid_layout = QHBoxLayout()
right_layout = QHBoxLayout()
mid_layout.addWidget(self.syn_checkbox)
mid_layout.addWidget(self.comp_checkbox)
right_layout.addWidget(self.qcombo_processed)
mid_widget.setLayout(mid_layout)
right_widget.setLayout(right_layout)
self.label_layout.addWidget(mid_widget)
self.label_layout.addWidget(right_widget)
self.syn_checkbox.setLayoutDirection(Qt.RightToLeft)
self.comp_checkbox.setLayoutDirection(Qt.RightToLeft)
self.label_layout.setStretch(0, 4)
self.label_layout.setStretch(1, 0)
self.label_layout.setStretch(2, 0)
@ -827,7 +815,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
label = QtWidgets.QLabel()
self.perm_labels.append(label)
self.qcombo_processed = QtWidgets.QComboBox()
self.syn_checkbox = QtWidgets.QCheckBox('synthetics')
self.comp_checkbox = QtWidgets.QCheckBox('Load comparison data')
self.addQCboxItem('processed', 'green')
self.addQCboxItem('raw', 'black')
# self.perm_qcbox_right.setAlignment(2)
@ -836,9 +824,11 @@ class WaveformWidgetPG(QtWidgets.QWidget):
def getPlotDict(self):
return self.plotdict
def activateObspyDMToptions(self, activate):
self.syn_checkbox.setVisible(activate)
self.qcombo_processed.setVisible(activate)
def activateObspyDMToptions(self, activate: bool) -> None:
self.qcombo_processed.setEnabled(activate)
def activateCompareOptions(self, activate: bool) -> None:
self.comp_checkbox.setEnabled(activate)
def setPermText(self, number, text=None, color='black'):
if not 0 <= number < len(self.perm_labels):
@ -961,7 +951,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
[time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else [])
trace.data = np.array(
[datum * gain + n for index, datum in enumerate(trace.data) if not index % nth_sample])
trace_syn.data = np.array([datum + n for index, datum in enumerate(trace_syn.data)
trace_syn.data = np.array([datum + n + shift_syn for index, datum in enumerate(trace_syn.data)
if not index % nth_sample] if st_syn else [])
plots.append((times, trace.data,
times_syn, trace_syn.data))
@ -1007,15 +997,6 @@ class WaveformWidgetPG(QtWidgets.QWidget):
time_ax = np.linspace(time_ax[0], time_ax[-1], num=len(data))
return data, time_ax
# def getAxes(self):
# return self.axes
# def getXLims(self):
# return self.getAxes().get_xlim()
# def getYLims(self):
# return self.getAxes().get_ylim()
def setXLims(self, lims):
vb = self.plotWidget.getPlotItem().getViewBox()
vb.setXRange(float(lims[0]), float(lims[1]), padding=0)
@ -1169,8 +1150,6 @@ class PylotCanvas(FigureCanvas):
break
if not ax_check: return
# self.updateCurrentLimits() #maybe put this down to else:
# calculate delta (relative values in axis)
old_x, old_y = self.press_rel
xdiff = gui_event.x - old_x
@ -1373,110 +1352,145 @@ class PylotCanvas(FigureCanvas):
plot_positions[channel] = plot_pos
return plot_positions
def plotWFData(self, wfdata, title=None, zoomx=None, zoomy=None,
def plotWFData(self, wfdata, wfdata_compare=None, title=None, zoomx=None, zoomy=None,
noiselevel=None, scaleddata=False, mapping=True,
component='*', nth_sample=1, iniPick=None, verbosity=0,
plot_additional=False, additional_channel=None, scaleToChannel=None,
snr=None):
ax = self.prepare_plot()
self.clearPlotDict()
wfstart, wfend = self.get_wf_range(wfdata)
compclass = self.get_comp_class()
plot_streams = self.get_plot_streams(wfdata, wfdata_compare, component, compclass)
st_main = plot_streams['wfdata']['data']
if mapping:
plot_positions = self.calcPlotPositions(st_main, compclass)
nslc = self.get_sorted_nslc(st_main)
nmax = self.plot_traces(ax, plot_streams, nslc, wfstart, mapping, plot_positions,
scaleToChannel, noiselevel, scaleddata, nth_sample, verbosity)
if plot_additional and additional_channel:
self.plot_additional_trace(ax, wfdata, additional_channel, scaleToChannel,
scaleddata, nth_sample, wfstart)
self.finalize_plot(ax, wfstart, wfend, nmax, zoomx, zoomy, iniPick, title, snr)
def prepare_plot(self):
ax = self.axes[0]
ax.cla()
return ax
self.clearPlotDict()
wfstart, wfend = full_range(wfdata)
nmax = 0
def get_wf_range(self, wfdata):
return full_range(wfdata)
def get_comp_class(self):
settings = QSettings()
compclass = SetChannelComponents.from_qsettings(settings)
return SetChannelComponents.from_qsettings(settings)
if not component == '*':
alter_comp = compclass.getCompPosition(component)
# alter_comp = str(alter_comp[0])
st_select = wfdata.select(component=component)
st_select += wfdata.select(component=alter_comp)
else:
st_select = wfdata
if mapping:
plot_positions = self.calcPlotPositions(st_select, compclass)
# list containing tuples of network, station, channel and plot position (for sorting)
nslc = []
for plot_pos, trace in enumerate(st_select):
if not trace.stats.channel[-1] in ['Z', 'N', 'E', '1', '2', '3']:
print('Warning: Unrecognized channel {}'.format(trace.stats.channel))
continue
nslc.append(trace.get_id())
nslc.sort()
nslc.reverse()
def get_plot_streams(self, wfdata, wfdata_compare, component, compclass):
def get_wf_dict(data=Stream(), linecolor='k', offset=0., **plot_kwargs):
return dict(data=data, linecolor=linecolor, offset=offset, plot_kwargs=plot_kwargs)
linecolor = (0., 0., 0., 1.) if not self.style else self.style['linecolor']['rgba_mpl']
plot_streams = {
'wfdata': get_wf_dict(linecolor=linecolor, linewidth=0.7),
'wfdata_comp': get_wf_dict(offset=0.1, linecolor='b', alpha=0.7, linewidth=0.5)
}
if component != '*':
alter_comp = compclass.getCompPosition(component)
plot_streams['wfdata']['data'] = wfdata.select(component=component) + wfdata.select(component=alter_comp)
if wfdata_compare:
plot_streams['wfdata_comp']['data'] = wfdata_compare.select(
component=component) + wfdata_compare.select(component=alter_comp)
else:
plot_streams['wfdata']['data'] = wfdata
if wfdata_compare:
plot_streams['wfdata_comp']['data'] = wfdata_compare
return plot_streams
def get_sorted_nslc(self, st_main):
nslc = [trace.get_id() for trace in st_main if trace.stats.channel[-1] in ['Z', 'N', 'E', '1', '2', '3']]
nslc.sort(reverse=True)
return nslc
def plot_traces(self, ax, plot_streams, nslc, wfstart, mapping, plot_positions, scaleToChannel, noiselevel,
scaleddata, nth_sample, verbosity):
nmax = 0
for n, seed_id in enumerate(nslc):
network, station, location, channel = seed_id.split('.')
st = st_select.select(id=seed_id)
trace = st[0].copy()
for wf_name, wf_dict in plot_streams.items():
st_select = wf_dict.get('data')
if not st_select:
continue
trace = st_select.select(id=seed_id)[0].copy()
if mapping:
n = plot_positions[trace.stats.channel]
if n > nmax:
nmax = n
if verbosity:
msg = 'plotting %s channel of station %s' % (channel, station)
print(msg)
stime = trace.stats.starttime - wfstart
time_ax = prep_time_axis(stime, trace)
print(f'plotting {channel} channel of station {station}')
time_ax = prep_time_axis(trace.stats.starttime - wfstart, trace)
self.plot_trace(ax, trace, time_ax, wf_dict, n, scaleToChannel, noiselevel, scaleddata, nth_sample)
self.setPlotDict(n, seed_id)
return nmax
def plot_trace(self, ax, trace, time_ax, wf_dict, n, scaleToChannel, noiselevel, scaleddata, nth_sample):
if time_ax is not None:
if scaleToChannel:
self.scale_trace(trace, scaleToChannel)
scaleddata = True
if not scaleddata:
trace.detrend('constant')
trace.normalize(np.max(np.abs(trace.data)) * 2)
offset = wf_dict.get('offset')
times = [time for index, time in enumerate(time_ax) if not index % nth_sample]
data = [datum + n + offset for index, datum in enumerate(trace.data) if not index % nth_sample]
ax.axhline(n, color="0.5", lw=0.5)
ax.plot(times, data, color=wf_dict.get('linecolor'), **wf_dict.get('plot_kwargs'))
if noiselevel is not None:
self.plot_noise_level(ax, time_ax, noiselevel, channel, n, wf_dict.get('linecolor'))
def scale_trace(self, trace, scaleToChannel):
st_scale = wfdata.select(channel=scaleToChannel)
if st_scale:
tr = st_scale[0]
trace.detrend('constant')
trace.normalize(np.max(np.abs(tr.data)) * 2)
scaleddata = True
if not scaleddata:
trace.detrend('constant')
trace.normalize(np.max(np.abs(trace.data)) * 2)
times = [time for index, time in enumerate(time_ax) if not index % nth_sample]
data = [datum + n for index, datum in enumerate(trace.data) if not index % nth_sample]
ax.axhline(n, color="0.5", lw=0.5)
ax.plot(times, data, color=linecolor, linewidth=0.7)
if noiselevel is not None:
def plot_noise_level(self, ax, time_ax, noiselevel, channel, n, linecolor):
for level in [-noiselevel[channel], noiselevel[channel]]:
ax.plot([time_ax[0], time_ax[-1]],
[n + level, n + level],
color=linecolor,
linestyle='dashed')
self.setPlotDict(n, seed_id)
if plot_additional and additional_channel:
ax.plot([time_ax[0], time_ax[-1]], [n + level, n + level], color=linecolor, linestyle='dashed')
def plot_additional_trace(self, ax, wfdata, additional_channel, scaleToChannel, scaleddata, nth_sample, wfstart):
compare_stream = wfdata.select(channel=additional_channel)
if compare_stream:
trace = compare_stream[0]
if scaleToChannel:
st_scale = wfdata.select(channel=scaleToChannel)
if st_scale:
tr = st_scale[0]
trace.detrend('constant')
trace.normalize(np.max(np.abs(tr.data)) * 2)
self.scale_trace(trace, scaleToChannel)
scaleddata = True
if not scaleddata:
trace.detrend('constant')
trace.normalize(np.max(np.abs(trace.data)) * 2)
time_ax = prep_time_axis(stime, trace)
time_ax = prep_time_axis(trace.stats.starttime - wfstart, trace)
self.plot_additional_data(ax, trace, time_ax, nth_sample)
def plot_additional_data(self, ax, trace, time_ax, nth_sample):
times = [time for index, time in enumerate(time_ax) if not index % nth_sample]
p_data = compare_stream[0].data
# #normalize
# p_max = max(abs(p_data))
# p_data /= p_max
p_data = trace.data
for index in range(3):
ax.plot(times, p_data, color='red', alpha=0.5, linewidth=0.7)
p_data += 1
def finalize_plot(self, ax, wfstart, wfend, nmax, zoomx, zoomy, iniPick, title, snr):
if iniPick:
ax.vlines(iniPick, ax.get_ylim()[0], ax.get_ylim()[1],
colors='m', linestyles='dashed',
linewidth=2)
xlabel = 'seconds since {0}'.format(wfstart)
ax.vlines(iniPick, ax.get_ylim()[0], ax.get_ylim()[1], colors='m', linestyles='dashed', linewidth=2)
xlabel = f'seconds since {wfstart}'
ylabel = ''
self.updateWidget(xlabel, ylabel, title)
self.setXLims(ax, [0, wfend - wfstart])
@ -1486,15 +1500,14 @@ class PylotCanvas(FigureCanvas):
if zoomy is not None:
self.setYLims(ax, zoomy)
if snr is not None:
if snr < 2:
warning = 'LOW SNR'
if snr < 1.5:
warning = 'VERY LOW SNR'
ax.text(0.1, 0.9, 'WARNING - {}'.format(warning), ha='center', va='center', transform=ax.transAxes,
color='red')
self.plot_snr_warning(ax, snr)
self.draw()
def plot_snr_warning(self, ax, snr):
if snr < 2:
warning = 'LOW SNR' if snr >= 1.5 else 'VERY LOW SNR'
ax.text(0.1, 0.9, f'WARNING - {warning}', ha='center', va='center', transform=ax.transAxes, color='red')
@staticmethod
def getXLims(ax):
return ax.get_xlim()
@ -1846,8 +1859,8 @@ class PhaseDefaults(QtWidgets.QDialog):
class PickDlg(QDialog):
update_picks = QtCore.Signal(dict)
def __init__(self, parent=None, data=None, station=None, network=None, location=None, picks=None,
autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None,
def __init__(self, parent=None, data=None, data_compare=None, station=None, network=None, location=None, picks=None,
autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None, show_comp_data=False,
event=None, filteroptions=None, model=None, wftype=None):
super(PickDlg, self).__init__(parent, Qt.Window)
self.orig_parent = parent
@ -1856,6 +1869,7 @@ class PickDlg(QDialog):
# initialize attributes
self.parameter = parameter
self._embedded = embedded
self.showCompData = show_comp_data
self.station = station
self.network = network
self.location = location
@ -1894,11 +1908,32 @@ class PickDlg(QDialog):
else:
self.filteroptions = FILTERDEFAULTS
self.pick_block = False
# set attribute holding data
if data is None or not data:
try:
data = parent.get_data().get_wf_data().copy()
self.data = data.select(station=station)
except AttributeError as e:
errmsg = 'You either have to put in a data or an appropriate ' \
'parent (PyLoT MainWindow) object: {0}'.format(e)
raise Exception(errmsg)
else:
self.data = data
self.data_compare = data_compare
self.nextStation = QtWidgets.QCheckBox('Continue with next station ')
# comparison channel
self.compareChannel = QtWidgets.QComboBox()
self.compareChannel.activated.connect(self.resetPlot)
self.referenceChannel = QtWidgets.QComboBox()
self.referenceChannel.activated.connect(self.resetPlot)
# comparison channel
self.compareCB = QtWidgets.QCheckBox()
self.compareCB.setChecked(self.showCompData)
self.compareCB.clicked.connect(self.switchCompData)
self.compareCB.clicked.connect(self.resetPlot)
self.compareCB.setVisible(bool(self.data_compare))
# scale channel
self.scaleChannel = QtWidgets.QComboBox()
@ -1911,19 +1946,7 @@ class PickDlg(QDialog):
self.cur_xlim = None
self.cur_ylim = None
# set attribute holding data
if data is None or not data:
try:
data = parent.get_data().getWFData().copy()
self.data = data.select(station=station)
except AttributeError as e:
errmsg = 'You either have to put in a data or an appropriate ' \
'parent (PyLoT MainWindow) object: {0}'.format(e)
raise Exception(errmsg)
else:
self.data = data
self.stime, self.etime = full_range(self.getWFData())
self.stime, self.etime = full_range(self.get_wf_data())
# initialize plotting widget
self.multicompfig = PylotCanvas(parent=self, multicursor=True)
@ -1934,12 +1957,12 @@ class PickDlg(QDialog):
self.setupUi()
# fill compare and scale channels
self.compareChannel.addItem('-', None)
self.referenceChannel.addItem('-', None)
self.scaleChannel.addItem('individual', None)
for trace in self.getWFData():
for trace in self.get_wf_data():
channel = trace.stats.channel
self.compareChannel.addItem(channel, trace)
self.referenceChannel.addItem(channel, trace)
if not channel[-1] in ['Z', 'N', 'E', '1', '2', '3']:
print('Skipping unknown channel for scaling: {}'.format(channel))
continue
@ -1956,7 +1979,7 @@ class PickDlg(QDialog):
if self.wftype is not None:
title += ' | ({})'.format(self.wftype)
self.multicompfig.plotWFData(wfdata=self.getWFData(),
self.multicompfig.plotWFData(wfdata=self.get_wf_data(), wfdata_compare=self.get_wf_dataComp(),
title=title)
self.multicompfig.setZoomBorders2content()
@ -2132,8 +2155,11 @@ class PickDlg(QDialog):
_dialtoolbar.addWidget(est_label)
_dialtoolbar.addWidget(self.plot_arrivals_button)
_dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtWidgets.QLabel('Compare to channel: '))
_dialtoolbar.addWidget(self.compareChannel)
_dialtoolbar.addWidget(QtWidgets.QLabel('Plot reference channel: '))
_dialtoolbar.addWidget(self.referenceChannel)
_dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtWidgets.QLabel('Compare: '))
_dialtoolbar.addWidget(self.compareCB)
_dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtWidgets.QLabel('Scaling: '))
_dialtoolbar.addWidget(self.scaleChannel)
@ -2470,7 +2496,7 @@ class PickDlg(QDialog):
def activatePicking(self):
self.leave_rename_phase()
self.renamePhaseAction.setEnabled(False)
self.compareChannel.setEnabled(False)
self.referenceChannel.setEnabled(False)
self.scaleChannel.setEnabled(False)
phase = self.currentPhase
phaseID = self.getPhaseID(phase)
@ -2488,7 +2514,7 @@ class PickDlg(QDialog):
if self.autoFilterAction.isChecked():
for filteraction in [self.filterActionP, self.filterActionS]:
filteraction.setChecked(False)
self.filterWFData()
self.filter_wf_data()
self.draw()
else:
self.draw()
@ -2502,7 +2528,7 @@ class PickDlg(QDialog):
self.disconnectPressEvent()
self.multicompfig.connectEvents()
self.renamePhaseAction.setEnabled(True)
self.compareChannel.setEnabled(True)
self.referenceChannel.setEnabled(True)
self.scaleChannel.setEnabled(True)
self.connect_pick_delete()
self.draw()
@ -2572,9 +2598,15 @@ class PickDlg(QDialog):
def getGlobalLimits(self, ax, axis):
return self.multicompfig.getGlobalLimits(ax, axis)
def getWFData(self):
def get_wf_data(self):
return self.data
def get_wf_dataComp(self):
if self.showCompData:
return self.data_compare
else:
return Stream()
def selectWFData(self, channel):
component = channel[-1].upper()
wfdata = Stream()
@ -2584,17 +2616,17 @@ class PickDlg(QDialog):
return tr
if component == 'E' or component == 'N':
for trace in self.getWFData():
for trace in self.get_wf_data():
trace = selectTrace(trace, 'NE')
if trace:
wfdata.append(trace)
elif component == '1' or component == '2':
for trace in self.getWFData():
for trace in self.get_wf_data():
trace = selectTrace(trace, '12')
if trace:
wfdata.append(trace)
else:
wfdata = self.getWFData().select(component=component)
wfdata = self.get_wf_data().select(component=component)
return wfdata
def getPicks(self, picktype='manual'):
@ -2696,11 +2728,16 @@ class PickDlg(QDialog):
stime = self.getStartTime()
# copy data for plotting
data = self.getWFData().copy()
data = self.getPickPhases(data, phase)
data.normalize()
if not data:
# copy wfdata for plotting
wfdata = self.get_wf_data().copy()
wfdata_comp = self.get_wf_dataComp().copy()
wfdata = self.getPickPhases(wfdata, phase)
wfdata_comp = self.getPickPhases(wfdata_comp, phase)
for wfd in [wfdata, wfdata_comp]:
if wfd:
wfd.normalize()
if not wfdata:
QtWidgets.QMessageBox.warning(self, 'No channel to plot',
'No channel to plot for phase: {}. '
'Make sure to select the correct channels for P and S '
@ -2708,14 +2745,16 @@ class PickDlg(QDialog):
self.leave_picking_mode()
return
# filter data and trace on which is picked prior to determination of SNR
# filter wfdata and trace on which is picked prior to determination of SNR
filterphase = self.currentFilterPhase()
if filterphase:
filteroptions = self.getFilterOptions(filterphase).parseFilterOptions()
try:
data.detrend('linear')
data.filter(**filteroptions)
# wfdata.filter(**filteroptions)# MP MP removed filtering of original data
for wfd in [wfdata, wfdata_comp]:
if wfd:
wfd.detrend('linear')
wfd.filter(**filteroptions)
# wfdata.filter(**filteroptions)# MP MP removed filtering of original wfdata
except ValueError as e:
self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Information,
'Denied',
@ -2725,8 +2764,8 @@ class PickDlg(QDialog):
snr = []
noiselevels = {}
# determine SNR and noiselevel
for trace in data.traces:
st = data.select(channel=trace.stats.channel)
for trace in wfdata.traces:
st = wfdata.select(channel=trace.stats.channel)
stime_diff = trace.stats.starttime - stime
result = getSNR(st, (noise_win, gap_win, signal_win), ini_pick - stime_diff)
snr.append(result[0])
@ -2737,12 +2776,14 @@ class PickDlg(QDialog):
noiselevel = nfac
noiselevels[trace.stats.channel] = noiselevel
# prepare plotting of data
for trace in data:
# prepare plotting of wfdata
for wfd in [wfdata, wfdata_comp]:
if wfd:
for trace in wfd:
t = prep_time_axis(trace.stats.starttime - stime, trace)
inoise = getnoisewin(t, ini_pick, noise_win, gap_win)
trace = demeanTrace(trace, inoise)
# upscale trace data in a way that each trace is vertically zoomed to noiselevel*factor
# upscale trace wfdata in a way that each trace is vertically zoomed to noiselevel*factor
channel = trace.stats.channel
noiselevel = noiselevels[channel]
noiseScaleFactor = self.calcNoiseScaleFactor(noiselevel, zoomfactor=5.)
@ -2753,7 +2794,7 @@ class PickDlg(QDialog):
x_res = getResolutionWindow(mean_snr, parameter.get('extent'))
xlims = [ini_pick - x_res, ini_pick + x_res]
ylims = list(np.array([-.5, .5]) + [0, len(data) - 1])
ylims = list(np.array([-.5, .5]) + [0, len(wfdata) - 1])
title = self.getStation() + ' picking mode'
title += ' | SNR: {}'.format(mean_snr)
@ -2761,9 +2802,10 @@ class PickDlg(QDialog):
filtops_str = transformFilteroptions2String(filteroptions)
title += ' | Filteroptions: {}'.format(filtops_str)
plot_additional = bool(self.compareChannel.currentText())
additional_channel = self.compareChannel.currentText()
self.multicompfig.plotWFData(wfdata=data,
plot_additional = bool(self.referenceChannel.currentText())
additional_channel = self.referenceChannel.currentText()
self.multicompfig.plotWFData(wfdata=wfdata,
wfdata_compare=wfdata_comp,
title=title,
zoomx=xlims,
zoomy=ylims,
@ -2797,7 +2839,7 @@ class PickDlg(QDialog):
filteroptions = None
# copy and filter data for earliest and latest possible picks
wfdata = self.getWFData().copy().select(channel=channel)
wfdata = self.get_wf_data().copy().select(channel=channel)
if filteroptions:
try:
wfdata.detrend('linear')
@ -2844,7 +2886,7 @@ class PickDlg(QDialog):
minFMSNR = parameter.get('minFMSNR')
quality = get_quality_class(spe, parameter.get('timeerrorsP'))
if quality <= minFMweight and snr >= minFMSNR:
FM = fmpicker(self.getWFData().select(channel=channel).copy(), wfdata.copy(), parameter.get('fmpickwin'),
FM = fmpicker(self.get_wf_data().select(channel=channel).copy(), wfdata.copy(), parameter.get('fmpickwin'),
pick - stime_diff)
# save pick times for actual phase
@ -3136,7 +3178,7 @@ class PickDlg(QDialog):
def togglePickBlocker(self):
return not self.pick_block
def filterWFData(self, phase=None):
def filter_wf_data(self, phase=None):
if not phase:
phase = self.currentPhase
if self.getPhaseID(phase) == 'P':
@ -3154,7 +3196,8 @@ class PickDlg(QDialog):
self.cur_xlim = self.multicompfig.axes[0].get_xlim()
self.cur_ylim = self.multicompfig.axes[0].get_ylim()
# self.multicompfig.updateCurrentLimits()
data = self.getWFData().copy()
wfdata = self.get_wf_data().copy()
wfdata_comp = self.get_wf_dataComp().copy()
title = self.getStation()
if filter:
filtoptions = None
@ -3162,19 +3205,22 @@ class PickDlg(QDialog):
filtoptions = self.getFilterOptions(self.getPhaseID(phase), gui_filter=True).parseFilterOptions()
if filtoptions is not None:
data.detrend('linear')
data.taper(0.02, type='cosine')
data.filter(**filtoptions)
for wfd in [wfdata, wfdata_comp]:
if wfd:
wfd.detrend('linear')
wfd.taper(0.02, type='cosine')
wfd.filter(**filtoptions)
filtops_str = transformFilteroptions2String(filtoptions)
title += ' | Filteroptions: {}'.format(filtops_str)
if self.wftype is not None:
title += ' | ({})'.format(self.wftype)
plot_additional = bool(self.compareChannel.currentText())
additional_channel = self.compareChannel.currentText()
plot_additional = bool(self.referenceChannel.currentText())
additional_channel = self.referenceChannel.currentText()
scale_channel = self.scaleChannel.currentText()
self.multicompfig.plotWFData(wfdata=data, title=title,
self.multicompfig.plotWFData(wfdata=wfdata, wfdata_compare=wfdata_comp,
title=title,
zoomx=self.getXLims(),
zoomy=self.getYLims(),
plot_additional=plot_additional,
@ -3188,14 +3234,14 @@ class PickDlg(QDialog):
def filterP(self):
self.filterActionS.setChecked(False)
if self.filterActionP.isChecked():
self.filterWFData('P')
self.filter_wf_data('P')
else:
self.refreshPlot()
def filterS(self):
self.filterActionP.setChecked(False)
if self.filterActionS.isChecked():
self.filterWFData('S')
self.filter_wf_data('S')
else:
self.refreshPlot()
@ -3247,11 +3293,14 @@ class PickDlg(QDialog):
self.resetZoom()
self.refreshPlot()
def switchCompData(self):
self.showCompData = self.compareCB.isChecked()
def refreshPlot(self):
if self.autoFilterAction.isChecked():
self.filterActionP.setChecked(False)
self.filterActionS.setChecked(False)
# data = self.getWFData().copy()
# data = self.get_wf_data().copy()
# title = self.getStation()
filter = False
phase = None
@ -3751,8 +3800,8 @@ class TuneAutopicker(QWidget):
fnames = self.station_ids[self.get_current_station_id()]
if not fnames == self.fnames:
self.fnames = fnames
self.data.setWFData(fnames)
wfdat = self.data.getWFData() # all available streams
self.data.set_wf_data(fnames)
wfdat = self.data.get_wf_data() # all available streams
# remove possible underscores in station names
# wfdat = remove_underscores(wfdat)
# rotate misaligned stations to ZNE
@ -3857,12 +3906,14 @@ class TuneAutopicker(QWidget):
network = None
location = None
wfdata = self.data.getWFData()
wfdata = self.data.get_wf_data()
wfdata_comp = self.data.get_wf_dataComp()
metadata = self.parent().metadata
event = self.get_current_event()
filteroptions = self.parent().filteroptions
wftype = self.wftype if self.obspy_dmt else ''
self.pickDlg = PickDlg(self.parent(), data=wfdata.select(station=station).copy(),
data_comp=wfdata_comp.select(station=station).copy(),
station=station, network=network,
location=location, parameter=self.parameter,
picks=self.get_current_event_picks(station),
@ -3911,7 +3962,7 @@ class TuneAutopicker(QWidget):
for plotitem in self._manual_pick_plots:
self.clear_plotitem(plotitem)
self._manual_pick_plots = []
st = self.data.getWFData()
st = self.data.get_wf_data()
tr = st.select(station=self.get_current_station())[0]
starttime = tr.stats.starttime
# create two lists with figure names and subindices (for subplots) to get the correct axes