marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

WIP: Simplify data structure #39

Draft
sebastianw wants to merge 31 commits from 38-simplify-data-structure into develop
pull from: 38-simplify-data-structure
merge into: marcel:develop
Conversation 1 Commits 31 Files Changed 25 +1608 -2532
25 changed files with 1608 additions and 2532 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

806
PyLoT.py
View File

File diff suppressed because it is too large Load Diff

24
autoPyLoT.py
View File

@ -28,7 +28,7 @@ from pylot.core.util.dataprocessing import restitute_data, Metadata
from pylot.core.util.defaults import SEPARATOR from pylot.core.util.defaults import SEPARATOR
from pylot.core.util.event import Event from pylot.core.util.event import Event
from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.structure import DATASTRUCTURE
from pylot.core.util.utils import get_None, trim_station_components, check4gapsAndRemove, check4doubled, \ from pylot.core.util.utils import get_none, trim_station_components, check4gapsAndRemove, check4doubled, \
check4rotated check4rotated
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
@ -91,9 +91,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
sp=sp_info) sp=sp_info)
print(splash) print(splash)
parameter = get_None(parameter) parameter = get_none(parameter)
inputfile = get_None(inputfile) inputfile = get_none(inputfile)
eventid = get_None(eventid) eventid = get_none(eventid)
fig_dict = None fig_dict = None
fig_dict_wadatijack = None fig_dict_wadatijack = None
@ -119,13 +119,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
obspyDMT_wfpath = input_dict['obspyDMT_wfpath'] obspyDMT_wfpath = input_dict['obspyDMT_wfpath']
if not parameter: if not parameter:
if inputfile: if not inputfile:
print('Using default input parameter')
parameter = PylotParameter(inputfile) parameter = PylotParameter(inputfile)
# iplot = parameter['iplot']
else:
infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
print('Using default input file {}'.format(infile))
parameter = PylotParameter(infile)
else: else:
if not type(parameter) == PylotParameter: if not type(parameter) == PylotParameter:
print('Wrong input type for parameter: {}'.format(type(parameter))) print('Wrong input type for parameter: {}'.format(type(parameter)))
@ -154,7 +150,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
datastructure.setExpandFields(exf) datastructure.setExpandFields(exf)
# check if default location routine NLLoc is available and all stations are used # check if default location routine NLLoc is available and all stations are used
if get_None(parameter['nllocbin']) and station == 'all': if get_none(parameter['nllocbin']) and station == 'all':
locflag = 1 locflag = 1
# get NLLoc-root path # get NLLoc-root path
nllocroot = parameter.get('nllocroot') nllocroot = parameter.get('nllocroot')
@ -247,7 +243,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
pylot_event = Event(eventpath) # event should be path to event directory pylot_event = Event(eventpath) # event should be path to event directory
data.setEvtData(pylot_event) data.setEvtData(pylot_event)
if fnames == 'None': 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 # the following is necessary because within
# multiple event processing no event ID is provided # multiple event processing no event ID is provided
# in autopylot.in # in autopylot.in
@ -262,7 +258,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
now.minute) now.minute)
parameter.setParam(eventID=eventID) parameter.setParam(eventID=eventID)
else: else:
data.setWFData(fnames) data.set_wf_data(fnames)
eventpath = events[0] eventpath = events[0]
# now = datetime.datetime.now() # now = datetime.datetime.now()
@ -272,7 +268,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
# now.hour, # now.hour,
# now.minute) # now.minute)
parameter.setParam(eventID=eventid) parameter.setParam(eventID=eventid)
wfdat = data.getWFData() # all available streams wfdat = data.get_wf_data() # all available streams
if not station == 'all': if not station == 'all':
wfdat = wfdat.select(station=station) wfdat = wfdat.select(station=station)
if not wfdat: if not wfdat:

687
pylot/core/io/data.py
View File

@ -1,658 +1,69 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import copy
import os import os
from dataclasses import dataclass, field
from typing import List, Union
from PySide2.QtWidgets import QMessageBox from obspy import UTCDateTime
from obspy import read_events
from obspy.core import read, Stream, UTCDateTime
from obspy.core.event import Event as ObsPyEvent
from obspy.io.sac import SacIOError
import pylot.core.loc.focmec as focmec from pylot.core.io.event import EventData
import pylot.core.loc.hypodd as hypodd from pylot.core.io.waveformdata import WaveformData
import pylot.core.loc.velest as velest from pylot.core.util.dataprocessing import Metadata
from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \
picksdict_from_pilot, merge_picks, PylotParameter
from pylot.core.util.errors import FormatError, OverwriteError
from pylot.core.util.event import Event
from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
from pylot.core.util.utils import fnConstructor, full_range, check4rotated, \
check_for_gaps_and_merge, trim_station_components, check_for_nan
@dataclass
class Data(object): class Data:
""" event_data: EventData = field(default_factory=EventData)
Data container with attributes wfdata holding ~obspy.core.stream. waveform_data: WaveformData = field(default_factory=WaveformData)
metadata: Metadata = field(default_factory=Metadata)
:type parent: PySide2.QtWidgets.QWidget object, optional _parent: Union[None, 'QtWidgets.QWidget'] = None
:param parent: A PySide2.QtWidgets.QWidget object utilized when
called by a GUI to display a PySide2.QtWidgets.QMessageBox instead of printing
to standard out.
:type evtdata: ~obspy.core.event.Event object, optional
:param evtdata ~obspy.core.event.Event object containing all derived or
loaded event. Container object holding, e.g. phase arrivals, etc.
"""
def __init__(self, parent=None, evtdata=None): def __init__(self, parent=None, evtdata=None):
self._parent = parent self._parent = parent
if self.getParent(): self.event_data = EventData(evtdata)
self.comp = parent.getComponent() self.waveform_data = WaveformData()
else:
self.comp = 'Z'
self.wfdata = Stream()
self._new = False
if isinstance(evtdata, ObsPyEvent) or isinstance(evtdata, Event):
pass
elif isinstance(evtdata, dict):
evt = readPILOTEvent(**evtdata)
evtdata = evt
elif isinstance(evtdata, str):
try:
cat = read_events(evtdata)
if len(cat) != 1:
raise ValueError('ambiguous event information for file: '
'{file}'.format(file=evtdata))
evtdata = cat[0]
except TypeError as e:
if 'Unknown format for file' in e.message:
if 'PHASES' in evtdata:
picks = picksdict_from_pilot(evtdata)
evtdata = ObsPyEvent()
evtdata.picks = picks_from_picksdict(picks)
elif 'LOC' in evtdata:
raise NotImplementedError('PILOT location information '
'read support not yet '
'implemented.')
elif 'event.pkl' in evtdata:
evtdata = qml_from_obspyDMT(evtdata)
else:
raise e
else:
raise e
else: # create an empty Event object
self.setNew()
evtdata = ObsPyEvent()
evtdata.picks = []
self.evtdata = evtdata
self.wforiginal = None
self.cuttimes = None
self.dirty = False
self.processed = None
def __str__(self): def __str__(self):
return str(self.wfdata) return str(self.waveform_data.wfdata)
def __add__(self, other): def __add__(self, other):
assert isinstance(other, Data), "operands must be of same type 'Data'" if not isinstance(other, Data):
rs_id = self.get_evt_data().get('resource_id') raise TypeError("Operands must be of type 'Data'")
rs_id_other = other.get_evt_data().get('resource_id') if self.event_data.is_new() and other.event_data.is_new():
if other.isNew() and not self.isNew():
picks_to_add = other.get_evt_data().picks
old_picks = self.get_evt_data().picks
wf_ids_old = [pick.waveform_id for pick in old_picks]
for new_pick in picks_to_add:
wf_id = new_pick.waveform_id
if wf_id in wf_ids_old:
for old_pick in old_picks:
comparison = [old_pick.waveform_id == new_pick.waveform_id,
old_pick.phase_hint == new_pick.phase_hint,
old_pick.method_id == new_pick.method_id]
if all(comparison):
del (old_pick)
old_picks.append(new_pick)
elif not other.isNew() and self.isNew():
new = other + self
self.evtdata = new.get_evt_data()
elif self.isNew() and other.isNew():
pass pass
elif rs_id == rs_id_other: elif other.event_data.is_new():
other.setNew() new_picks = other.event_data.evtdata.picks
old_picks = self.event_data.evtdata.picks
old_picks.extend([pick for pick in new_picks if pick not in old_picks])
elif self.event_data.is_new():
return other + self
elif self.event_data.get_id() == other.event_data.get_id():
other.event_data.set_new()
return self + other return self + other
else: else:
raise ValueError("both Data objects have differing " raise ValueError("Both Data objects have differing unique Event identifiers")
"unique Event identifiers")
return self return self
def getPicksStr(self): def get_parent(self):
"""
Return picks in event data
:return: picks seperated by newlines
:rtype: str
"""
picks_str = ''
for pick in self.get_evt_data().picks:
picks_str += str(pick) + '\n'
return picks_str
def getParent(self):
"""
Get PySide.QtGui.QWidget parent object
"""
return self._parent return self._parent
def isNew(self): def filter_wf_data(self, **kwargs):
return self._new self.waveform_data.wfdata.detrend('linear')
self.waveform_data.wfdata.taper(0.02, type='cosine')
self.waveform_data.wfdata.filter(**kwargs)
self.waveform_data.dirty = True
def setNew(self): def set_wf_data(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0):
self._new = True return self.waveform_data.load_waveforms(fnames, fnames_alt, check_rotated, metadata, tstart, tstop)
def getCutTimes(self): def reset_wf_data(self):
""" self.waveform_data.reset()
Returns earliest start and latest end of all waveform data
:return: minimum start time and maximum end time as a tuple
:rtype: (UTCDateTime, UTCDateTime)
"""
if self.cuttimes is None:
self.updateCutTimes()
return self.cuttimes
def updateCutTimes(self): def get_wf_data(self):
""" return self.waveform_data.wfdata
Update cuttimes to contain earliest start and latest end time
of all waveform data
:rtype: None
"""
self.cuttimes = full_range(self.getWFData())
def getEventFileName(self): def rotate_wf_data(self):
ID = self.getID() self.waveform_data.rotate_zne(self.metadata)
# handle forbidden filenames especially on windows systems
return fnConstructor(str(ID))
def checkEvent(self, event, fcheck, forceOverwrite=False):
"""
Check information in supplied event and own event and replace with own
information if no other information are given or forced by forceOverwrite
:param event: Event that supplies information for comparison
:type event: pylot.core.util.event.Event
:param fcheck: check and delete existing information
can be a str or a list of strings of ['manual', 'auto', 'origin', 'magnitude']
:type fcheck: str, [str]
:param forceOverwrite: Set to true to force overwrite own information. If false,
supplied information from event is only used if there is no own information in that
category (given in fcheck: manual, auto, origin, magnitude)
:type forceOverwrite: bool
:return:
:rtype: None
"""
if 'origin' in fcheck:
self.replaceOrigin(event, forceOverwrite)
if 'magnitude' in fcheck:
self.replaceMagnitude(event, forceOverwrite)
if 'auto' in fcheck:
self.replacePicks(event, 'auto')
if 'manual' in fcheck:
self.replacePicks(event, 'manual')
def replaceOrigin(self, event, forceOverwrite=False):
"""
Replace own origin with the one supplied in event if own origin is not
existing or forced by forceOverwrite = True
:param event: Event that supplies information for comparison
:type event: pylot.core.util.event.Event
:param forceOverwrite: always replace own information with supplied one if true
:type forceOverwrite: bool
:return:
:rtype: None
"""
if self.get_evt_data().origins or forceOverwrite:
if event.origins:
print("Found origin, replace it by new origin.")
event.origins = self.get_evt_data().origins
def replaceMagnitude(self, event, forceOverwrite=False):
"""
Replace own magnitude with the one supplied in event if own magnitude is not
existing or forced by forceOverwrite = True
:param event: Event that supplies information for comparison
:type event: pylot.core.util.event.Event
:param forceOverwrite: always replace own information with supplied one if true
:type forceOverwrite: bool
:return:
:rtype: None
"""
if self.get_evt_data().magnitudes or forceOverwrite:
if event.magnitudes:
print("Found magnitude, replace it by new magnitude")
event.magnitudes = self.get_evt_data().magnitudes
def replacePicks(self, event, picktype):
"""
Replace picks in event with own picks
:param event: Event that supplies information for comparison
:type event: pylot.core.util.event.Event
:param picktype: 'auto' or 'manual' picks
:type picktype: str
:return:
:rtype: None
"""
checkflag = 1
picks = event.picks
# remove existing picks
for j, pick in reversed(list(enumerate(picks))):
try:
if picktype in str(pick.method_id.id):
picks.pop(j)
checkflag = 2
except AttributeError as e:
msg = '{}'.format(e)
print(e)
checkflag = 0
if checkflag > 0:
if checkflag == 1:
print("Write new %s picks to catalog." % picktype)
if checkflag == 2:
print("Found %s pick(s), remove them and append new picks to catalog." % picktype)
# append new picks
for pick in self.get_evt_data().picks:
if picktype in str(pick.method_id.id):
picks.append(pick)
def exportEvent(self, fnout, fnext='.xml', fcheck='auto', upperErrors=None):
"""
Export event to file
:param fnout: basename of file
:param fnext: file extensions xml, cnv, obs, focmec, or/and pha
:param fcheck: check and delete existing information
can be a str or a list of strings of ['manual', 'auto', 'origin', 'magnitude']
"""
from pylot.core.util.defaults import OUTPUTFORMATS
if not type(fcheck) == list:
fcheck = [fcheck]
try:
evtformat = OUTPUTFORMATS[fnext]
except KeyError as e:
errmsg = '{0}; selected file extension {1} not ' \
'supported'.format(e, fnext)
raise FormatError(errmsg)
if hasattr(self.get_evt_data(), 'notes'):
try:
with open(os.path.join(os.path.dirname(fnout), 'notes.txt'), 'w') as notes_file:
notes_file.write(self.get_evt_data().notes)
except Exception as e:
print('Warning: Could not save notes.txt: ', str(e))
# check for already existing xml-file
if fnext == '.xml':
if os.path.isfile(fnout + fnext):
print("xml-file already exists! Check content ...")
cat = read_events(fnout + fnext)
if len(cat) > 1:
raise IOError('Ambigious event information in file {}'.format(fnout + fnext))
if len(cat) < 1:
raise IOError('No event information in file {}'.format(fnout + fnext))
event = cat[0]
if not event.resource_id == self.get_evt_data().resource_id:
QMessageBox.warning(self, 'Warning', 'Different resource IDs!')
return
self.checkEvent(event, fcheck)
self.setEvtData(event)
self.get_evt_data().write(fnout + fnext, format=evtformat)
# try exporting event
else:
evtdata_org = self.get_evt_data()
picks = evtdata_org.picks
eventpath = evtdata_org.path
picks_copy = copy.deepcopy(picks)
evtdata_copy = Event(eventpath)
evtdata_copy.picks = picks_copy
# check for stations picked automatically as well as manually
# Prefer manual picks!
for i in range(len(picks)):
if picks[i].method_id == 'manual':
mstation = picks[i].waveform_id.station_code
mstation_ext = mstation + '_'
for k in range(len(picks_copy)):
if ((picks_copy[k].waveform_id.station_code == mstation) or
(picks_copy[k].waveform_id.station_code == mstation_ext)) and \
(picks_copy[k].method_id == 'auto'):
del picks_copy[k]
break
lendiff = len(picks) - len(picks_copy)
if lendiff != 0:
print("Manual as well as automatic picks available. Prefered the {} manual ones!".format(lendiff))
no_uncertainties_p = []
no_uncertainties_s = []
if upperErrors:
# check for pick uncertainties exceeding adjusted upper errors
# Picks with larger uncertainties will not be saved in output file!
for j in range(len(picks)):
for i in range(len(picks_copy)):
if picks_copy[i].phase_hint[0] == 'P':
# Skipping pick if no upper_uncertainty is found and warning user
if picks_copy[i].time_errors['upper_uncertainty'] is None:
#print("{1} P-Pick of station {0} does not have upper_uncertainty and cant be checked".format(
# picks_copy[i].waveform_id.station_code,
# picks_copy[i].method_id))
if not picks_copy[i].waveform_id.station_code in no_uncertainties_p:
no_uncertainties_p.append(picks_copy[i].waveform_id.station_code)
continue
#print ("checking for upper_uncertainty")
if (picks_copy[i].time_errors['uncertainty'] is None) or \
(picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[0]):
print("Uncertainty exceeds or equal adjusted upper time error!")
print("Adjusted uncertainty: {}".format(upperErrors[0]))
print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty']))
print("{1} P-Pick of station {0} will not be saved in outputfile".format(
picks_copy[i].waveform_id.station_code,
picks_copy[i].method_id))
del picks_copy[i]
break
if picks_copy[i].phase_hint[0] == 'S':
# Skipping pick if no upper_uncertainty is found and warning user
if picks_copy[i].time_errors['upper_uncertainty'] is None:
#print("{1} S-Pick of station {0} does not have upper_uncertainty and cant be checked".format(
#picks_copy[i].waveform_id.station_code,
#picks_copy[i].method_id))
if not picks_copy[i].waveform_id.station_code in no_uncertainties_s:
no_uncertainties_s.append(picks_copy[i].waveform_id.station_code)
continue
if (picks_copy[i].time_errors['uncertainty'] is None) or \
(picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[1]):
print("Uncertainty exceeds or equal adjusted upper time error!")
print("Adjusted uncertainty: {}".format(upperErrors[1]))
print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty']))
print("{1} S-Pick of station {0} will not be saved in outputfile".format(
picks_copy[i].waveform_id.station_code,
picks_copy[i].method_id))
del picks_copy[i]
break
for s in no_uncertainties_p:
print("P-Pick of station {0} does not have upper_uncertainty and cant be checked".format(s))
for s in no_uncertainties_s:
print("S-Pick of station {0} does not have upper_uncertainty and cant be checked".format(s))
if fnext == '.obs':
try:
evtdata_copy.write(fnout + fnext, format=evtformat)
# write header afterwards
evid = str(evtdata_org.resource_id).split('/')[1]
header = '# EQEVENT: Label: EQ%s Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n' % evid
nllocfile = open(fnout + fnext)
l = nllocfile.readlines()
# Adding A0/Generic Amplitude to .obs file
# l2 = []
# for li in l:
# for amp in evtdata_org.amplitudes:
# if amp.waveform_id.station_code == li[0:5].strip():
# li = li[0:64] + '{:0.2e}'.format(amp.generic_amplitude) + li[73:-1] + '\n'
# l2.append(li)
# l = l2
nllocfile.close()
l.insert(0, header)
nllocfile = open(fnout + fnext, 'w')
nllocfile.write("".join(l))
nllocfile.close()
except KeyError as e:
raise KeyError('''{0} export format
not implemented: {1}'''.format(evtformat, e))
if fnext == '.cnv':
try:
velest.export(picks_copy, fnout + fnext, eventinfo=self.get_evt_data())
except KeyError as e:
raise KeyError('''{0} export format
not implemented: {1}'''.format(evtformat, e))
if fnext == '_focmec.in':
try:
infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
print('Using default input file {}'.format(infile))
parameter = PylotParameter(infile)
focmec.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data())
except KeyError as e:
raise KeyError('''{0} export format
not implemented: {1}'''.format(evtformat, e))
if fnext == '.pha':
try:
infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
print('Using default input file {}'.format(infile))
parameter = PylotParameter(infile)
hypodd.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data())
except KeyError as e:
raise KeyError('''{0} export format
not implemented: {1}'''.format(evtformat, e))
def getComp(self):
"""
Get component (ZNE)
"""
return self.comp
def getID(self):
"""
Get unique resource id
"""
try:
return self.evtdata.get('resource_id').id
except:
return None
def filterWFData(self, kwargs):
"""
Filter waveform data
:param kwargs: arguments to pass through to filter function
"""
data = self.getWFData()
data.detrend('linear')
data.taper(0.02, type='cosine')
data.filter(**kwargs)
self.dirty = True
def setWFData(self, fnames, fnames_syn=None, checkRotated=False, metadata=None, tstart=0, tstop=0):
"""
Clear current waveform data and set given waveform data
:param fnames: waveform data names to append
:type fnames: list
"""
def check_fname_exists(filenames: list) -> list:
if filenames:
filenames = [fn for fn in filenames if os.path.isfile(fn)]
return filenames
self.wfdata = Stream()
self.wforiginal = None
self.wfsyn = Stream()
if tstart == tstop:
tstart = tstop = None
self.tstart = tstart
self.tstop = tstop
fnames = check_fname_exists(fnames)
fnames_syn = check_fname_exists(fnames_syn)
# if obspy_dmt:
# wfdir = 'raw'
# self.processed = False
# for fname in fnames:
# if fname.endswith('processed'):
# wfdir = 'processed'
# self.processed = True
# break
# for fpath in fnames:
# if fpath.endswith(wfdir):
# wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
# if 'syngine' in fpath.split('/')[-1]:
# wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
# else:
# wffnames = fnames
if fnames is not None:
self.appendWFData(fnames)
if fnames_syn is not None:
self.appendWFData(fnames_syn, synthetic=True)
else:
return False
# various pre-processing steps:
# remove possible underscores in station names
# self.wfdata = remove_underscores(self.wfdata)
# check for gaps and merge
self.wfdata, _ = check_for_gaps_and_merge(self.wfdata)
# check for nans
check_for_nan(self.wfdata)
# check for stations with rotated components
if checkRotated and metadata is not None:
self.wfdata = check4rotated(self.wfdata, metadata, verbosity=0)
# trim station components to same start value
trim_station_components(self.wfdata, trim_start=True, trim_end=False)
# make a copy of original data
self.wforiginal = self.getWFData().copy()
self.dirty = False
return True
def appendWFData(self, fnames, synthetic=False):
"""
Read waveform data from fnames and append it to current wf data
:param fnames: waveform data to append
:type fnames: list
"""
assert isinstance(fnames, list), "input parameter 'fnames' is " \
"supposed to be of type 'list' " \
"but is actually" \
" {0}".format(type(fnames))
if self.dirty:
self.resetWFData()
real_or_syn_data = {True: self.wfsyn,
False: self.wfdata}
warnmsg = ''
for fname in set(fnames):
try:
real_or_syn_data[synthetic] += read(fname, starttime=self.tstart, endtime=self.tstop)
except TypeError:
try:
real_or_syn_data[synthetic] += read(fname, format='GSE2', starttime=self.tstart, endtime=self.tstop)
except Exception as e:
try:
real_or_syn_data[synthetic] += read(fname, format='SEGY', starttime=self.tstart,
endtime=self.tstop)
except Exception as e:
warnmsg += '{0}\n{1}\n'.format(fname, e)
except SacIOError as se:
warnmsg += '{0}\n{1}\n'.format(fname, se)
if warnmsg:
warnmsg = 'WARNING in appendWFData: unable to read waveform data\n' + warnmsg
print(warnmsg)
def getWFData(self):
return self.wfdata
def getOriginalWFData(self):
return self.wforiginal
def getSynWFData(self):
return self.wfsyn
def resetWFData(self):
"""
Set waveform data to original waveform data
"""
if self.getOriginalWFData():
self.wfdata = self.getOriginalWFData().copy()
else:
self.wfdata = Stream()
self.dirty = False
def resetPicks(self):
"""
Clear all picks from event
"""
self.get_evt_data().picks = []
def get_evt_data(self):
return self.evtdata
def setEvtData(self, event):
self.evtdata = event
def applyEVTData(self, data, typ='pick'):
"""
Either takes an `obspy.core.event.Event` object and applies all new
information on the event to the actual data if typ is 'event or
creates ObsPy pick objects and append it to the picks list from the
PyLoT dictionary contain all picks if type is pick
:param data: data to apply, either picks or complete event
:type data:
:param typ: which event data to apply, 'pick' or 'event'
:type typ: str
:param authority_id: (currently unused)
:type: str
:raise OverwriteError:
"""
def applyPicks(picks):
"""
Creates ObsPy pick objects and append it to the picks list from the
PyLoT dictionary contain all picks.
:param picks:
:raise OverwriteError: raises an OverwriteError if the picks list is
not empty. The GUI will then ask for a decision.
"""
# firstonset = find_firstonset(picks)
# check for automatic picks
print("Writing phases to ObsPy-quakeml file")
for key in picks:
if not picks[key].get('P'):
continue
if picks[key]['P']['picker'] == 'auto':
print("Existing auto-picks will be overwritten in pick-dictionary!")
picks = picks_from_picksdict(picks)
break
else:
if self.get_evt_data().picks:
raise OverwriteError('Existing picks would be overwritten!')
else:
picks = picks_from_picksdict(picks)
break
self.get_evt_data().picks = picks
# if 'smi:local' in self.getID() and firstonset:
# fonset_str = firstonset.strftime('%Y_%m_%d_%H_%M_%S')
# ID = ResourceIdentifier('event/' + fonset_str)
# ID.convertIDToQuakeMLURI(authority_id=authority_id)
# self.get_evt_data().resource_id = ID
def applyEvent(event):
"""
takes an `obspy.core.event.Event` object and applies all new
information on the event to the actual data
:param event:
"""
if event is None:
print("applyEvent: Received None")
return
if self.isNew():
self.setEvtData(event)
else:
# prevent overwriting original pick information
event_old = self.get_evt_data()
if not event_old.resource_id == event.resource_id:
print("WARNING: Missmatch in event resource id's: {} and {}".format(
event_old.resource_id,
event.resource_id))
else:
picks = copy.deepcopy(event_old.picks)
event = merge_picks(event, picks)
# apply event information from location
event_old.update(event)
applydata = {'pick': applyPicks,
'event': applyEvent}
applydata[typ](data)
self._new = False
class GenericDataStructure(object): class GenericDataStructure(object):
@ -837,22 +248,6 @@ class PilotDataStructure(GenericDataStructure):
self.setExpandFields(['root', 'database']) self.setExpandFields(['root', 'database'])
class ObspyDMTdataStructure(GenericDataStructure):
"""
Object containing the data access information for the old PILOT data
structure.
"""
def __init__(self, **fields):
if not fields:
fields = {'database': '',
'root': ''}
GenericDataStructure.__init__(self, **fields)
self.setExpandFields(['root', 'database'])
class SeiscompDataStructure(GenericDataStructure): class SeiscompDataStructure(GenericDataStructure):
""" """
Dictionary containing the data access information for an SDS data archive: Dictionary containing the data access information for an SDS data archive:

106
pylot/core/io/event.py Normal file
View File

@ -0,0 +1,106 @@
import copy
from dataclasses import dataclass
from typing import Union
from obspy import read_events
from obspy.core.event import Event as ObsPyEvent
@dataclass
class EventData:
evtdata: Union[ObsPyEvent, None] = None
_new: bool = False
def __init__(self, evtdata=None):
self.set_event_data(evtdata)
def set_event_data(self, evtdata):
if isinstance(evtdata, ObsPyEvent):
self.evtdata = evtdata
elif isinstance(evtdata, dict):
self.evtdata = self.read_pilot_event(**evtdata)
elif isinstance(evtdata, str):
self.evtdata = self.read_event_file(evtdata)
else:
self.set_new()
self.evtdata = ObsPyEvent(picks=[])
def read_event_file(self, evtdata: str) -> ObsPyEvent:
try:
cat = read_events(evtdata)
if len(cat) != 1:
raise ValueError(f'ambiguous event information for file: {evtdata}')
return cat[0]
except TypeError as e:
self.handle_event_file_error(e, evtdata)
def handle_event_file_error(self, e: TypeError, evtdata: str):
if 'Unknown format for file' in str(e):
if 'PHASES' in evtdata:
picks = self.picksdict_from_pilot(evtdata)
evtdata = ObsPyEvent(picks=self.picks_from_picksdict(picks))
elif 'LOC' in evtdata:
raise NotImplementedError('PILOT location information read support not yet implemented.')
elif 'event.pkl' in evtdata:
evtdata = self.qml_from_obspy_dmt(evtdata)
else:
raise e
else:
raise e
def set_new(self):
self._new = True
def is_new(self) -> bool:
return self._new
def get_picks_str(self) -> str:
return '\n'.join(str(pick) for pick in self.evtdata.picks)
def replace_origin(self, event: ObsPyEvent, force_overwrite: bool = False):
if self.evtdata.origins or force_overwrite:
event.origins = self.evtdata.origins
def replace_magnitude(self, event: ObsPyEvent, force_overwrite: bool = False):
if self.evtdata.magnitudes or force_overwrite:
event.magnitudes = self.evtdata.magnitudes
def replace_picks(self, event: ObsPyEvent, picktype: str):
checkflag = 1
picks = event.picks
for j, pick in reversed(list(enumerate(picks))):
if picktype in str(pick.method_id.id):
picks.pop(j)
checkflag = 2
if checkflag > 0:
for pick in self.evtdata.picks:
if picktype in str(pick.method_id.id):
picks.append(pick)
def get_id(self) -> Union[str, None]:
try:
return self.evtdata.resource_id.id
except:
return None
def apply_event_data(self, data, typ='pick'):
if typ == 'pick':
self.apply_picks(data)
elif typ == 'event':
self.apply_event(data)
def apply_picks(self, picks):
self.evtdata.picks = picks
def apply_event(self, event: ObsPyEvent):
if self.is_new():
self.evtdata = event
else:
old_event = self.evtdata
if old_event.resource_id == event.resource_id:
picks = copy.deepcopy(old_event.picks)
event = self.merge_picks(event, picks)
old_event.update(event)
else:
print(f"WARNING: Mismatch in event resource id's: {old_event.resource_id} and {event.resource_id}")

6
pylot/core/io/location.py
View File

@ -1,7 +1,7 @@
from obspy import UTCDateTime from obspy import UTCDateTime
from obspy.core import event as ope from obspy.core import event as ope
from pylot.core.util.utils import getLogin, getHash from pylot.core.util.utils import get_login, get_hash
def create_amplitude(pickID, amp, unit, category, cinfo): def create_amplitude(pickID, amp, unit, category, cinfo):
@ -61,7 +61,7 @@ def create_creation_info(agency_id=None, creation_time=None, author=None):
:return: :return:
''' '''
if author is None: if author is None:
author = getLogin() author = get_login()
if creation_time is None: if creation_time is None:
creation_time = UTCDateTime() creation_time = UTCDateTime()
return ope.CreationInfo(agency_id=agency_id, author=author, return ope.CreationInfo(agency_id=agency_id, author=author,
@ -210,7 +210,7 @@ def create_resourceID(timetohash, restype, authority_id=None, hrstr=None):
''' '''
assert isinstance(timetohash, UTCDateTime), "'timetohash' is not an ObsPy" \ assert isinstance(timetohash, UTCDateTime), "'timetohash' is not an ObsPy" \
"UTCDateTime object" "UTCDateTime object"
hid = getHash(timetohash) hid = get_hash(timetohash)
if hrstr is None: if hrstr is None:
resID = ope.ResourceIdentifier(restype + '/' + hid[0:6]) resID = ope.ResourceIdentifier(restype + '/' + hid[0:6])
else: else:

848
pylot/core/io/phases.py
View File

@ -1,6 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import glob import glob
import logging
import os import os
import warnings import warnings
@ -16,29 +17,10 @@ from pylot.core.io.inputs import PylotParameter
from pylot.core.io.location import create_event, \ from pylot.core.io.location import create_event, \
create_magnitude create_magnitude
from pylot.core.pick.utils import select_for_phase, get_quality_class from pylot.core.pick.utils import select_for_phase, get_quality_class
from pylot.core.util.utils import getOwner, full_range, four_digits, transformFilterString4Export, \ from pylot.core.util.utils import get_owner, full_range, four_digits, transformFilterString4Export, \
backtransformFilterString, loopIdentifyPhase, identifyPhase backtransformFilterString, loopIdentifyPhase, identifyPhase
def add_amplitudes(event, amplitudes):
amplitude_list = []
for pick in event.picks:
try:
a0 = amplitudes[pick.waveform_id.station_code]
amplitude = ope.Amplitude(generic_amplitude=a0 * 1e-3)
amplitude.unit = 'm'
amplitude.category = 'point'
amplitude.waveform_id = pick.waveform_id
amplitude.magnitude_hint = 'ML'
amplitude.pick_id = pick.resource_id
amplitude.type = 'AML'
amplitude_list.append(amplitude)
except KeyError:
continue
event.amplitudes = amplitude_list
return event
def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs): def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
""" """
readPILOTEvent - function readPILOTEvent - function
@ -58,7 +40,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
if phasfn is not None and os.path.isfile(phasfn): if phasfn is not None and os.path.isfile(phasfn):
phases = sio.loadmat(phasfn) phases = sio.loadmat(phasfn)
phasctime = UTCDateTime(os.path.getmtime(phasfn)) phasctime = UTCDateTime(os.path.getmtime(phasfn))
phasauthor = getOwner(phasfn) phasauthor = get_owner(phasfn)
else: else:
phases = None phases = None
phasctime = None phasctime = None
@ -66,7 +48,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
if locfn is not None and os.path.isfile(locfn): if locfn is not None and os.path.isfile(locfn):
loc = sio.loadmat(locfn) loc = sio.loadmat(locfn)
locctime = UTCDateTime(os.path.getmtime(locfn)) locctime = UTCDateTime(os.path.getmtime(locfn))
locauthor = getOwner(locfn) locauthor = get_owner(locfn)
else: else:
loc = None loc = None
locctime = None locctime = None
@ -192,32 +174,7 @@ def convert_pilot_times(time_array):
return UTCDateTime(*times) return UTCDateTime(*times)
def picksdict_from_obs(fn): def picksdict_from_picks(evt, parameter=None):
"""
create pick dictionary from obs file
:param fn: filename
:type fn:
:return:
:rtype:
"""
picks = dict()
station_name = str()
for line in open(fn, 'r'):
if line.startswith('#'):
continue
else:
phase_line = line.split()
if not station_name == phase_line[0]:
phase = dict()
station_name = phase_line[0]
phase_name = phase_line[4].upper()
pick = UTCDateTime(phase_line[6] + phase_line[7] + phase_line[8])
phase[phase_name] = dict(mpp=pick, fm=phase_line[5])
picks[station_name] = phase
return picks
def picksdict_from_picks(evt):
""" """
Takes an Event object and return the pick dictionary commonly used within Takes an Event object and return the pick dictionary commonly used within
PyLoT PyLoT
@ -230,6 +187,7 @@ def picksdict_from_picks(evt):
'auto': {} 'auto': {}
} }
for pick in evt.picks: for pick in evt.picks:
errors = None
phase = {} phase = {}
station = pick.waveform_id.station_code station = pick.waveform_id.station_code
if pick.waveform_id.channel_code is None: if pick.waveform_id.channel_code is None:
@ -273,33 +231,29 @@ def picksdict_from_picks(evt):
phase['epp'] = epp phase['epp'] = epp
phase['lpp'] = lpp phase['lpp'] = lpp
phase['spe'] = spe phase['spe'] = spe
try: weight = phase.get('weight')
phase['weight'] = weight if not weight:
except: if not parameter:
# get onset weight from uncertainty logging.warning('Using ')
infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in') logging.warning('Using default input parameter')
print('Using default input file {}'.format(infile)) parameter = PylotParameter()
parameter = PylotParameter(infile)
pick.phase_hint = identifyPhase(pick.phase_hint) pick.phase_hint = identifyPhase(pick.phase_hint)
if pick.phase_hint == 'P': if pick.phase_hint == 'P':
errors = parameter['timeerrorsP'] errors = parameter['timeerrorsP']
elif pick.phase_hint == 'S': elif pick.phase_hint == 'S':
errors = parameter['timeerrorsS'] errors = parameter['timeerrorsS']
if errors:
weight = get_quality_class(spe, errors) weight = get_quality_class(spe, errors)
phase['weight'] = weight phase['weight'] = weight
phase['channel'] = channel phase['channel'] = channel
phase['network'] = network phase['network'] = network
phase['picker'] = pick_method phase['picker'] = pick_method
try:
if pick.polarity == 'positive': if pick.polarity == 'positive':
phase['fm'] = 'U' phase['fm'] = 'U'
elif pick.polarity == 'negative': elif pick.polarity == 'negative':
phase['fm'] = 'D' phase['fm'] = 'D'
else: else:
phase['fm'] = 'N' phase['fm'] = 'N'
except:
print("No FM info available!")
phase['fm'] = 'N'
phase['filter_id'] = filter_id if filter_id is not None else '' phase['filter_id'] = filter_id if filter_id is not None else ''
onsets[pick.phase_hint] = phase.copy() onsets[pick.phase_hint] = phase.copy()
@ -375,636 +329,228 @@ def picks_from_picksdict(picks, creation_info=None):
return picks_list return picks_list
def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0): def write_phases(arrivals, fformat, filename, parameter=None, eventinfo=None):
# TODO: change root to datapath
db_root = os.path.join(root_dir, db_dir)
evt_list = glob.glob1(db_root, 'e????.???.??')
for evt in evt_list:
if verbosity > 0:
print('Reassessing event {0}'.format(evt))
reassess_pilot_event(root_dir, db_dir, evt, out_dir, fn_param, verbosity)
def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None, verbosity=0):
from obspy import read
from pylot.core.io.inputs import PylotParameter
from pylot.core.pick.utils import earllatepicker
# TODO: change root to datapath
default = PylotParameter(fn_param, verbosity)
search_base = os.path.join(root_dir, db_dir, event_id)
phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat'))
if not phases_file:
return
if verbosity > 1:
print('Opening PILOT phases file: {fn}'.format(fn=phases_file[0]))
picks_dict = picksdict_from_pilot(phases_file[0])
if verbosity > 0:
print('Dictionary read from PHASES.mat:\n{0}'.format(picks_dict))
datacheck = list()
info = None
for station in picks_dict.keys():
fn_pattern = os.path.join(search_base, '{0}*'.format(station))
try:
st = read(fn_pattern)
except TypeError as e:
if 'Unknown format for file' in e.message:
try:
st = read(fn_pattern, format='GSE2')
except ValueError as e:
if e.message == 'second must be in 0..59':
info = 'A known Error was raised. Please find the list of corrupted files and double-check these files.'
datacheck.append(fn_pattern + ' (time info)\n')
continue
else:
raise ValueError(e.message)
except Exception as e:
if 'No file matching file pattern:' in e.message:
if verbosity > 0:
warnings.warn('no waveform data found for station {station}'.format(station=station),
RuntimeWarning)
datacheck.append(fn_pattern + ' (no data)\n')
continue
else:
raise e
else:
raise e
for phase in picks_dict[station].keys():
try:
mpp = picks_dict[station][phase]['mpp']
except KeyError as e:
print(e.message, station)
continue
sel_st = select_for_phase(st, phase)
if not sel_st:
msg = 'no waveform data found for station {station}'.format(station=station)
warnings.warn(msg, RuntimeWarning)
continue
stime, etime = full_range(sel_st)
rel_pick = mpp - stime
epp, lpp, spe = earllatepicker(sel_st,
default.get('nfac{0}'.format(phase)),
default.get('tsnrz' if phase == 'P' else 'tsnrh'),
Pick1=rel_pick,
iplot=0,
verbosity=0)
if epp is None or lpp is None:
continue
epp = stime + epp
lpp = stime + lpp
min_diff = 3 * st[0].stats.delta
if lpp - mpp < min_diff:
lpp = mpp + min_diff
if mpp - epp < min_diff:
epp = mpp - min_diff
picks_dict[station][phase] = dict(epp=epp, mpp=mpp, lpp=lpp, spe=spe)
if datacheck:
if info:
if verbosity > 0:
print(info + ': {0}'.format(search_base))
fncheck = open(os.path.join(search_base, 'datacheck_list'), 'w')
fncheck.writelines(datacheck)
fncheck.close()
del datacheck
# create Event object for export
evt = ope.Event(resource_id=event_id)
evt.picks = picks_from_picksdict(picks_dict)
# write phase information to file
if not out_dir:
fnout_prefix = os.path.join(root_dir, db_dir, event_id, 'PyLoT_{0}.'.format(event_id))
else:
out_dir = os.path.join(out_dir, db_dir)
if not os.path.isdir(out_dir):
os.makedirs(out_dir)
fnout_prefix = os.path.join(out_dir, 'PyLoT_{0}.'.format(event_id))
evt.write(fnout_prefix + 'xml', format='QUAKEML')
def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
""" """
Function of methods to write phases to the following standard file Writes earthquake phase data to different file formats.
formats used for locating earthquakes:
HYPO71, NLLoc, VELEST, HYPOSAT, FOCMEC, and hypoDD :param arrivals: Dictionary containing phase information (station ID, phase, first motion, weight, etc.)
:param arrivals:dictionary containing all phase information including
station ID, phase, first motion, weight (uncertainty), ...
:type arrivals: dict :type arrivals: dict
:param fformat: File format to write to (e.g., 'NLLoc', 'HYPO71', 'HYPOSAT', 'VELEST', 'HYPODD', 'FOCMEC')
:param fformat: chosen file format (location routine),
choose between NLLoc, HYPO71, HYPOSAT, VELEST,
HYPOINVERSE, FOCMEC, and hypoDD
:type fformat: str :type fformat: str
:param filename: Path and name of the output phase file
:param filename: full path and name of phase file :type filename: str
:type filename: string :param parameter: Additional parameters for writing the phase data
:param parameter: all input information
:type parameter: object :type parameter: object
:param eventinfo: Event information needed for specific formats like VELEST, FOCMEC, and HASH
:param eventinfo: optional, needed for VELEST-cnv file :type eventinfo: obspy.core.event.Event
and FOCMEC- and HASH-input files
:type eventinfo: `obspy.core.event.Event` object
""" """
if fformat == 'NLLoc':
print("Writing phases to %s for NLLoc" % filename)
fid = open("%s" % filename, 'w')
# write header
fid.write('# EQEVENT: %s Label: EQ%s Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n' %
(parameter.get('database'), parameter.get('eventID')))
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
for key in arrivals:
# P onsets
if 'P' in arrivals[key]:
try:
fm = arrivals[key]['P']['fm']
except KeyError as e:
print(e)
fm = None
if fm is None:
fm = '?'
onset = arrivals[key]['P']['mpp']
year = onset.year
month = onset.month
day = onset.day
hh = onset.hour
mm = onset.minute
ss = onset.second
ms = onset.microsecond
ss_ms = ss + ms / 1000000.0
pweight = 1 # use pick
try:
if arrivals[key]['P']['weight'] >= 4:
pweight = 0 # do not use pick
print("Station {}: Uncertain pick, do not use it!".format(key))
except KeyError as e:
print(e.message + '; no weight set during processing')
fid.write('%s ? ? ? P %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 %d \n' % (key,
fm,
year,
month,
day,
hh,
mm,
ss_ms,
pweight))
# S onsets
if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None:
fm = '?'
onset = arrivals[key]['S']['mpp']
year = onset.year
month = onset.month
day = onset.day
hh = onset.hour
mm = onset.minute
ss = onset.second
ms = onset.microsecond
ss_ms = ss + ms / 1000000.0
sweight = 1 # use pick
try:
if arrivals[key]['S']['weight'] >= 4:
sweight = 0 # do not use pick
except KeyError as e:
print(str(e) + '; no weight set during processing')
Ao = arrivals[key]['S']['Ao'] # peak-to-peak amplitude
if Ao == None:
Ao = 0.0
# fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 %d \n' % (key,
fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 %9.2f 0 0 %d \n' % (key,
fm,
year,
month,
day,
hh,
mm,
ss_ms,
Ao,
sweight))
fid.close() def write_nlloc():
elif fformat == 'HYPO71': with open(filename, 'w') as fid:
print("Writing phases to %s for HYPO71" % filename) fid.write('# EQEVENT: {} Label: EQ{} Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n'.format(
fid = open("%s" % filename, 'w') parameter.get('database'), parameter.get('eventID')))
# write header for key, value in arrivals.items():
fid.write(' %s\n' % for phase in ['P', 'S']:
parameter.get('eventID')) if phase in value:
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere fm = value[phase].get('fm', '?')
for key in arrivals: onset = value[phase]['mpp']
if arrivals[key]['P']['weight'] < 4: ss_ms = onset.second + onset.microsecond / 1000000.0
stat = key weight = 1 if value[phase].get('weight', 0) < 4 else 0
if len(stat) > 4: # HYPO71 handles only 4-string station IDs amp = value[phase].get('Ao', 0.0) if phase == 'S' else ''
stat = stat[1:5] fid.write('{} ? ? ? {} {}{}{} {}{} {:7.4f} GAU 0 {} 0 0 {}\n'.format(
Ponset = arrivals[key]['P']['mpp'] key, phase, fm, onset.year, onset.month, onset.day, onset.hour, onset.minute, ss_ms, amp,
Sonset = arrivals[key]['S']['mpp'] weight))
pweight = arrivals[key]['P']['weight']
sweight = arrivals[key]['S']['weight']
fm = arrivals[key]['P']['fm']
if fm is None:
fm = '-'
Ao = arrivals[key]['S']['Ao']
if Ao is None:
Ao = ''
else:
Ao = str('%7.2f' % Ao)
year = Ponset.year
if year >= 2000:
year = year - 2000
else:
year = year - 1900
month = Ponset.month
day = Ponset.day
hh = Ponset.hour
mm = Ponset.minute
ss = Ponset.second
ms = Ponset.microsecond
ss_ms = ss + ms / 1000000.0
if pweight < 2:
pstr = 'I'
elif pweight >= 2:
pstr = 'E'
if arrivals[key]['S']['weight'] < 4:
Sss = Sonset.second
Sms = Sonset.microsecond
Sss_ms = Sss + Sms / 1000000.0
Sss_ms = str('%5.02f' % Sss_ms)
if sweight < 2:
sstr = 'I'
elif sweight >= 2:
sstr = 'E'
fid.write('%-4s%sP%s%d %02d%02d%02d%02d%02d%5.2f %s%sS %d %s\n' % (stat,
pstr,
fm,
pweight,
year,
month,
day,
hh,
mm,
ss_ms,
Sss_ms,
sstr,
sweight,
Ao))
else:
fid.write('%-4s%sP%s%d %02d%02d%02d%02d%02d%5.2f %s\n' % (stat,
pstr,
fm,
pweight,
year,
month,
day,
hh,
mm,
ss_ms,
Ao))
fid.close() def write_hypo71():
with open(filename, 'w') as fid:
elif fformat == 'HYPOSAT': fid.write(
print("Writing phases to %s for HYPOSAT" % filename) ' {}\n'.format(parameter.get('eventID')))
fid = open("%s" % filename, 'w') for key, value in arrivals.items():
# write header if value['P'].get('weight', 0) < 4:
fid.write('%s, event %s \n' % (parameter.get('database'), parameter.get('eventID'))) stat = key[:4]
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere Ponset = value['P']['mpp']
for key in arrivals: Sonset = value.get('S', {}).get('mpp')
# P onsets pweight = value['P'].get('weight', 0)
if 'P' in arrivals[key] and arrivals[key]['P']['mpp'] is not None: sweight = value.get('S', {}).get('weight', 0)
if arrivals[key]['P']['weight'] < 4: fm = value['P'].get('fm', '-')
Ponset = arrivals[key]['P']['mpp'] Ao = value.get('S', {}).get('Ao', '')
pyear = Ponset.year year = Ponset.year - 2000 if Ponset.year >= 2000 else Ponset.year - 1900
pmonth = Ponset.month ss_ms = Ponset.second + Ponset.microsecond / 1000000.0
pday = Ponset.day if Sonset:
phh = Ponset.hour Sss_ms = Sonset.second + Sonset.microsecond / 1000000.0
pmm = Ponset.minute fid.write('{}P{}{}{} {}{}{}{}{} {:5.2f} {}{}S {} {}\n'.format(
pss = Ponset.second stat, 'I' if pweight < 2 else 'E', fm, pweight, year, Ponset.month, Ponset.day,
pms = Ponset.microsecond Ponset.hour, Ponset.minute, ss_ms, Sss_ms, 'I' if sweight < 2 else 'E', sweight, Ao))
Pss = pss + pms / 1000000.0
# use symmetrized picking error as std
# (read the HYPOSAT manual)
pstd = arrivals[key]['P']['spe']
if pstd is None:
errorsP = parameter.get('timeerrorsP')
if arrivals[key]['P']['weight'] == 0:
pstd = errorsP[0]
elif arrivals[key]['P']['weight'] == 1:
pstd = errorsP[1]
elif arrivals[key]['P']['weight'] == 2:
pstd = errorsP[2]
elif arrivals[key]['P']['weight'] == 3:
psrd = errorsP[3]
else: else:
pstd = errorsP[4] fid.write('{}P{}{}{} {}{}{}{}{} {:5.2f} {}\n'.format(
fid.write('%-5s P1 %4.0f %02d %02d %02d %02d %05.02f %5.3f -999. 0.00 -999. 0.00\n' stat, 'I' if pweight < 2 else 'E', fm, pweight, year, Ponset.month, Ponset.day,
% (key, pyear, pmonth, pday, phh, pmm, Pss, pstd)) Ponset.hour, Ponset.minute, ss_ms, Ao))
# S onsets
if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None:
if arrivals[key]['S']['weight'] < 4:
Sonset = arrivals[key]['S']['mpp']
syear = Sonset.year
smonth = Sonset.month
sday = Sonset.day
shh = Sonset.hour
smm = Sonset.minute
sss = Sonset.second
sms = Sonset.microsecond
Sss = sss + sms / 1000000.0
sstd = arrivals[key]['S']['spe']
if pstd is None:
errorsS = parameter.get('timeerrorsS')
if arrivals[key]['S']['weight'] == 0:
pstd = errorsS[0]
elif arrivals[key]['S']['weight'] == 1:
pstd = errorsS[1]
elif arrivals[key]['S']['weight'] == 2:
pstd = errorsS[2]
elif arrivals[key]['S']['weight'] == 3:
psrd = errorsS[3]
else:
pstd = errorsP[4]
fid.write('%-5s S1 %4.0f %02d %02d %02d %02d %05.02f %5.3f -999. 0.00 -999. 0.00\n'
% (key, syear, smonth, sday, shh, smm, Sss, sstd))
fid.close()
elif fformat == 'VELEST': def write_hyposat():
print("Writing phases to %s for VELEST" % filename) with open(filename, 'w') as fid:
fid = open("%s" % filename, 'w') fid.write('{}, event {} \n'.format(parameter.get('database'), parameter.get('eventID')))
# get informations needed in cnv-file for key, value in arrivals.items():
# check, whether latitude is N or S and longitude is E or W for phase in ['P', 'S']:
try: if phase in value and value[phase].get('weight', 0) < 4:
eventsource = eventinfo.origins[0] onset = value[phase]['mpp']
except: ss_ms = onset.second + onset.microsecond / 1000000.0
std = value[phase].get('spe', parameter.get('timeerrorsP')[value[phase].get('weight', 0)])
fid.write(
'{:<5} {}1 {:4} {:02} {:02} {:02} {:02} {:05.02f} {:5.3f} -999. 0.00 -999. 0.00\n'.format(
key, phase, onset.year, onset.month, onset.day, onset.hour, onset.minute, ss_ms, std))
def write_velest():
if not eventinfo:
print("No source origin calculated yet, thus no cnv-file creation possible!") print("No source origin calculated yet, thus no cnv-file creation possible!")
return return
if eventsource['latitude'] < 0: with open(filename, 'w') as fid:
cns = 'S' origin = eventinfo.origins[0]
else: lat_dir = 'S' if origin.latitude < 0 else 'N'
cns = 'N' lon_dir = 'W' if origin.longitude < 0 else 'E'
if eventsource['longitude'] < 0: year = origin.time.year - 2000 if origin.time.year >= 2000 else origin.time.year - 1900
cew = 'W' fid.write(
else: '{}{}{} {}{} {} {:05.2f} {:7.4f}{} {:8.4f}{} {:7.2f} {:6.2f} {:02.0f} 0.0 0.03 1.0 1.0\n'.format(
cew = 'E' year, origin.time.month, origin.time.day, origin.time.hour, origin.time.minute, origin.time.second,
# get last two integers of origin year origin.latitude, lat_dir, origin.longitude, lon_dir, origin.depth, eventinfo.magnitudes[0].mag, 0))
stime = eventsource['time'] for key, value in arrivals.items():
if stime.year - 2000 >= 0: for phase in ['P', 'S']:
syear = stime.year - 2000 if phase in value and value[phase].get('weight', 0) < 4:
else: onset = value[phase]['mpp']
syear = stime.year - 1900 rt = (onset - origin.time).total_seconds()
ifx = 0 # default value, see VELEST manual, pp. 22-23 fid.write('{:<4}{}{}{:6.2f}\n'.format(key[:4], phase, value[phase].get('weight', 0), rt))
# write header
fid.write('%s%02d%02d %02d%02d %05.2f %7.4f%c %8.4f%c %7.2f %6.2f %02.0f 0.0 0.03 1.0 1.0\n' % (
syear, stime.month, stime.day, stime.hour, stime.minute, stime.second, eventsource['latitude'],
cns, eventsource['longitude'], cew, eventsource['depth'], eventinfo.magnitudes[0]['mag'], ifx))
n = 0
# check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT)
if isinstance(arrivals, dict) == False:
# convert pick object (PyLoT) into dictionary
evt = ope.Event(resource_id=eventinfo['resource_id'])
evt.picks = arrivals
arrivals = picksdict_from_picks(evt)
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
# P onsets
if 'P' in usedarrivals[key]:
if usedarrivals[key]['P']['weight'] < 4:
n += 1
stat = key
if len(stat) > 4: # VELEST handles only 4-string station IDs
stat = stat[1:5]
Ponset = usedarrivals[key]['P']['mpp']
Pweight = usedarrivals[key]['P']['weight']
Prt = Ponset - stime # onset time relative to source time
if n % 6 != 0:
fid.write('%-4sP%d%6.2f' % (stat, Pweight, Prt))
else:
fid.write('%-4sP%d%6.2f\n' % (stat, Pweight, Prt))
# S onsets
if 'S' in usedarrivals[key]:
if usedarrivals[key]['S']['weight'] < 4:
n += 1
stat = key
if len(stat) > 4: # VELEST handles only 4-string station IDs
stat = stat[1:5]
Sonset = usedarrivals[key]['S']['mpp']
Sweight = usedarrivals[key]['S']['weight']
Srt = Ponset - stime # onset time relative to source time
if n % 6 != 0:
fid.write('%-4sS%d%6.2f' % (stat, Sweight, Srt))
else:
fid.write('%-4sS%d%6.2f\n' % (stat, Sweight, Srt))
fid.close()
elif fformat == 'HYPODD': def write_hypodd():
print("Writing phases to %s for hypoDD" % filename) if not eventinfo:
fid = open("%s" % filename, 'w')
# get event information needed for hypoDD-phase file
try:
eventsource = eventinfo.origins[0]
except:
print("No source origin calculated yet, thus no hypoDD-infile creation possible!") print("No source origin calculated yet, thus no hypoDD-infile creation possible!")
return return
stime = eventsource['time'] with open(filename, 'w') as fid:
try: origin = eventinfo.origins[0]
event = eventinfo['pylot_id'] stime = origin.time
hddID = event.split('.')[0][1:5] fid.write('# {} {} {} {} {} {} {:7.4f} +{:6.4f} {:7.4f} {:4.2f} 0.1 0.5 {:4.2f} {}\n'.format(
except:
print("Error 1111111!")
hddID = "00000"
# write header
fid.write('# %d %d %d %d %d %5.2f %7.4f +%6.4f %7.4f %4.2f 0.1 0.5 %4.2f %s\n' % (
stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second, stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second,
eventsource['latitude'], eventsource['longitude'], eventsource['depth'] / 1000, origin.latitude, origin.longitude, origin.depth / 1000, eventinfo.magnitudes[0].mag,
eventinfo.magnitudes[0]['mag'], eventsource['quality']['standard_error'], hddID)) origin.quality.standard_error, "00000"))
# check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT) for key, value in arrivals.items():
if isinstance(arrivals, dict) == False: for phase in ['P', 'S']:
# convert pick object (PyLoT) into dictionary if phase in value and value[phase].get('weight', 0) < 4:
evt = ope.Event(resource_id=eventinfo['resource_id']) onset = value[phase]['mpp']
evt.picks = arrivals rt = (onset - stime).total_seconds()
arrivals = picksdict_from_picks(evt) fid.write('{} {:6.3f} 1 {}\n'.format(key, rt, phase))
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
if 'P' in usedarrivals[key]:
# P onsets
if usedarrivals[key]['P']['weight'] < 4:
Ponset = usedarrivals[key]['P']['mpp']
Prt = Ponset - stime # onset time relative to source time
fid.write('%s %6.3f 1 P\n' % (key, Prt))
if 'S' in usedarrivals[key]:
# S onsets
if usedarrivals[key]['S']['weight'] < 4:
Sonset = usedarrivals[key]['S']['mpp']
Srt = Sonset - stime # onset time relative to source time
fid.write('%-5s %6.3f 1 S\n' % (key, Srt))
fid.close() def write_focmec():
if not eventinfo:
elif fformat == 'FOCMEC':
print("Writing phases to %s for FOCMEC" % filename)
fid = open("%s" % filename, 'w')
# get event information needed for FOCMEC-input file
try:
eventsource = eventinfo.origins[0]
except:
print("No source origin calculated yet, thus no FOCMEC-infile creation possible!") print("No source origin calculated yet, thus no FOCMEC-infile creation possible!")
return return
stime = eventsource['time'] with open(filename, 'w') as fid:
origin = eventinfo.origins[0]
# avoid printing '*' in focmec-input file stime = origin.time
if parameter.get('eventid') == '*' or parameter.get('eventid') is None: fid.write('{} {}{:02d}{:02d}{:02d}{:02d}{:02.0f} {:7.4f} {:6.4f} {:3.1f} {:3.1f}\n'.format(
evID = 'e0000' parameter.get('eventid', 'e0000'), stime.year, stime.month, stime.day, stime.hour, stime.minute,
else: stime.second, origin.latitude, origin.longitude, origin.depth / 1000, eventinfo.magnitudes[0].mag))
evID = parameter.get('eventid') for key, value in arrivals.items():
if 'P' in value and value['P'].get('weight', 0) < 4 and value['P'].get('fm'):
# write header line including event information for pick in eventinfo.picks:
fid.write('%s %d%02d%02d%02d%02d%02.0f %7.4f %6.4f %3.1f %3.1f\n' % (evID, if pick.waveform_id.station_code == key:
stime.year, stime.month, stime.day, for arrival in origin.arrivals:
stime.hour, stime.minute, stime.second, if arrival.pick_id == pick.resource_id and arrival.phase == 'P':
eventsource['latitude'], stat = key[:4]
eventsource['longitude'], az = arrival.azimuth
eventsource['depth'] / 1000, inz = arrival.takeoff_angle
eventinfo.magnitudes[0]['mag'])) fid.write('{:<4} {:6.2f} {:6.2f}{}\n'.format(stat, az, inz, value['P']['fm']))
picks = eventinfo.picks
# check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT)
if isinstance(arrivals, dict) == False:
# convert pick object (PyLoT) into dictionary
evt = ope.Event(resource_id=eventinfo['resource_id'])
evt.picks = arrivals
arrivals = picksdict_from_picks(evt)
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
if 'P' in usedarrivals[key]:
if usedarrivals[key]['P']['weight'] < 4 and usedarrivals[key]['P']['fm'] is not None:
stat = key
for i in range(len(picks)):
station = picks[i].waveform_id.station_code
if station == stat:
# get resource ID
resid_picks = picks[i].get('resource_id')
# find same ID in eventinfo
# there it is the pick_id!!
for j in range(len(eventinfo.origins[0].arrivals)):
resid_eventinfo = eventinfo.origins[0].arrivals[j].get('pick_id')
if resid_eventinfo == resid_picks and eventinfo.origins[0].arrivals[j].phase == 'P':
if len(stat) > 4: # FOCMEC handles only 4-string station IDs
stat = stat[1:5]
az = eventinfo.origins[0].arrivals[j].get('azimuth')
inz = eventinfo.origins[0].arrivals[j].get('takeoff_angle')
fid.write('%-4s %6.2f %6.2f%s \n' % (stat,
az,
inz,
usedarrivals[key]['P']['fm']))
break break
fid.close() def write_hash():
# Define filenames for HASH driver 1 and 2
filename1 = f"{filename}drv1.phase"
filename2 = f"{filename}drv2.phase"
elif fformat == 'HASH': print(f"Writing phases to {filename1} for HASH-driver 1")
# two different input files for print(f"Writing phases to {filename2} for HASH-driver 2")
# HASH-driver 1 and 2 (see HASH manual!)
filename1 = filename + 'drv1' + '.phase' # Open files for writing
filename2 = filename + 'drv2' + '.phase' with open(filename1, 'w') as fid1, open(filename2, 'w') as fid2:
print("Writing phases to %s for HASH for HASH-driver 1" % filename1) # Get event information needed for HASH-input file
fid1 = open("%s" % filename1, 'w')
print("Writing phases to %s for HASH for HASH-driver 2" % filename2)
fid2 = open("%s" % filename2, 'w')
# get event information needed for HASH-input file
try: try:
eventsource = eventinfo.origins[0] eventsource = eventinfo.origins[0]
except: except IndexError:
print("No source origin calculated yet, thus no cnv-file creation possible!") print("No source origin calculated yet, thus no cnv-file creation possible!")
return return
eventsource = eventinfo.origins[0]
event = parameter.get('eventID') event = parameter.get('eventID')
hashID = event.split('.')[0][1:5] hashID = event.split('.')[0][1:5]
latdeg = eventsource['latitude'] latdeg = eventsource['latitude']
latmin = eventsource['latitude'] * 60 / 10000 latmin = (eventsource['latitude'] * 60) / 10000
londeg = eventsource['longitude'] londeg = eventsource['longitude']
lonmin = eventsource['longitude'] * 60 / 10000 lonmin = (eventsource['longitude'] * 60) / 10000
erh = 1 / 2 * (eventsource.origin_uncertainty['min_horizontal_uncertainty'] +
eventsource.origin_uncertainty['max_horizontal_uncertainty']) / 1000 erh = (eventsource.origin_uncertainty['min_horizontal_uncertainty'] +
eventsource.origin_uncertainty['max_horizontal_uncertainty']) / 2000
erz = eventsource.depth_errors['uncertainty'] erz = eventsource.depth_errors['uncertainty']
stime = eventsource['time'] stime = eventsource['time']
if stime.year - 2000 >= 0: syear = stime.year % 100 # Calculate two-digit year
syear = stime.year - 2000
else:
syear = stime.year - 1900
picks = eventinfo.picks picks = eventinfo.picks
# write header line including event information
# for HASH-driver 1 # Write header line including event information for HASH-driver 1
fid1.write('%s%02d%02d%02d%02d%5.2f%2dN%5.2f%3dE%5.2f%6.3f%4.2f%5.2f%5.2f%s\n' % (syear, fid1.write(f"{syear:02d}{stime.month:02d}{stime.day:02d}{stime.hour:02d}{stime.minute:02d}"
stime.month, stime.day, f"{stime.second:05.2f}{latdeg:2d}N{latmin:05.2f}{londeg:3d}E{lonmin:05.2f}"
stime.hour, stime.minute, f"{eventsource['depth']:6.2f}{eventinfo.magnitudes[0]['mag']:4.2f}{erh:5.2f}{erz:5.2f}{hashID}\n")
stime.second,
latdeg, latmin, londeg, # Write header line including event information for HASH-driver 2
lonmin, eventsource['depth'], fid2.write(f"{syear:02d}{stime.month:02d}{stime.day:02d}{stime.hour:02d}{stime.minute:02d}"
eventinfo.magnitudes[0][ f"{stime.second:05.2f}{latdeg}N{latmin:05.2f}{londeg}E{lonmin:6.2f}{eventsource['depth']:5.2f}"
'mag'], erh, erz, f"{eventsource['quality']['used_phase_count']:3d}{erh:5.2f}{erz:5.2f}"
hashID)) f"{eventinfo.magnitudes[0]['mag']:4.2f}{hashID}\n")
# write header line including event information
# for HASH-driver 2 # Write phase lines
fid2.write( for key, arrival in arrivals.items():
'%d%02d%02d%02d%02d%5.2f%dN%5.2f%3dE%6.2f%5.2f %d %5.2f %5.2f %4.2f %s \n' % ( if 'P' in arrival and arrival['P']['weight'] < 4 and arrival['P']['fm'] is not None:
syear, stime.month, stime.day,
stime.hour, stime.minute, stime.second,
latdeg, latmin, londeg, lonmin,
eventsource['depth'],
eventsource['quality']['used_phase_count'],
erh, erz, eventinfo.magnitudes[0]['mag'],
hashID))
# Prefer Manual Picks over automatic ones if possible
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
# write phase lines
for key in arrivals:
if 'P' in arrivals[key]:
if arrivals[key]['P']['weight'] < 4 and arrivals[key]['P']['fm'] is not None:
stat = key stat = key
ccode = arrivals[key]['P']['channel'] ccode = arrival['P']['channel']
ncode = arrivals[key]['P']['network'] ncode = arrival['P']['network']
Pqual = 'I' if arrival['P']['weight'] < 2 else 'E'
if arrivals[key]['P']['weight'] < 2: for pick in picks:
Pqual = 'I' if pick.waveform_id.station_code == stat:
else: resid_picks = pick.get('resource_id')
Pqual = 'E' for origin_arrival in eventinfo.origins[0].arrivals:
if (origin_arrival.get('pick_id') == resid_picks and
for i in range(len(picks)): origin_arrival.phase == 'P'):
station = picks[i].waveform_id.station_code if len(stat) > 4: # HASH handles only 4-character station IDs
if station == stat:
# get resource ID
resid_picks = picks[i].get('resource_id')
# find same ID in eventinfo
# there it is the pick_id!!
for j in range(len(eventinfo.origins[0].arrivals)):
resid_eventinfo = eventinfo.origins[0].arrivals[j].get('pick_id')
if resid_eventinfo == resid_picks and eventinfo.origins[0].arrivals[j].phase == 'P':
if len(stat) > 4: # HASH handles only 4-string station IDs
stat = stat[1:5] stat = stat[1:5]
az = eventinfo.origins[0].arrivals[j].get('azimuth')
inz = eventinfo.origins[0].arrivals[j].get('takeoff_angle') az = origin_arrival.get('azimuth')
dist = eventinfo.origins[0].arrivals[j].get('distance') inz = origin_arrival.get('takeoff_angle')
# write phase line for HASH-driver 1 dist = origin_arrival.get('distance')
fid1.write(
'%-4s%sP%s%d 0 %3.1f %03d %03d 2 1 %s\n' % ( # Write phase line for HASH-driver 1
stat, Pqual, arrivals[key]['P']['fm'], arrivals[key]['P']['weight'], fid1.write(f"{stat:<4}{Pqual}P{arrival['P']['fm']}{arrival['P']['weight']:d}"
dist, inz, az, ccode)) f"{dist:3.1f}{inz:03d}{az:03d}{ccode}\n")
# write phase line for HASH-driver 2
fid2.write('%-4s %s %s %s %s \n' % ( # Write phase line for HASH-driver 2
stat, fid2.write(f"{stat:<4} {ncode} {ccode} {Pqual} {arrival['P']['fm']}\n")
ncode,
ccode,
Pqual,
arrivals[key]['P']['fm']))
break break
fid1.write(' %s' % hashID) fid1.write(f"{'':<36}{hashID}")
fid1.close()
fid2.close() # Prefer Manual Picks over automatic ones if possible
arrivals = chooseArrivals(arrivals) # Function not defined, assumed to exist
if fformat == 'NLLoc':
write_nlloc()
elif fformat == 'HYPO71':
write_hypo71()
elif fformat == 'HYPOSAT':
write_hyposat()
elif fformat == 'VELEST':
write_velest()
elif fformat == 'HYPODD':
write_hypodd()
elif fformat == 'FOCMEC':
write_focmec()
elif fformat == 'HASH':
write_hash()
def chooseArrivals(arrivals): def chooseArrivals(arrivals):

177
pylot/core/io/project.py Normal file
View File

@ -0,0 +1,177 @@
import os
from pylot.core.util.event import Event
class Project(object):
'''
Pickable class containing information of a PyLoT project, like event lists and file locations.
'''
# TODO: remove rootpath
def __init__(self):
self.eventlist = []
self.location = None
self.rootpath = None
self.datapath = None
self.dirty = False
self.parameter = None
self._table = None
def add_eventlist(self, eventlist):
'''
Add events from an eventlist containing paths to event directories.
Will skip existing paths.
'''
if len(eventlist) == 0:
return
for item in eventlist:
event = Event(item)
event.rootpath = self.parameter['rootpath']
event.database = self.parameter['database']
event.datapath = self.parameter['datapath']
if not event.path in self.getPaths():
self.eventlist.append(event)
self.setDirty()
else:
print('Skipping event with path {}. Already part of project.'.format(event.path))
self.eventlist.sort(key=lambda x: x.pylot_id)
self.search_eventfile_info()
def remove_event(self, event):
self.eventlist.remove(event)
def remove_event_by_id(self, eventID):
for event in self.eventlist:
if eventID in str(event.resource_id):
self.remove_event(event)
break
def read_eventfile_info(self, filename, separator=','):
'''
Try to read event information from file (:param:filename) comparing specific event datetimes.
File structure (each row): event, date, time, magnitude, latitude, longitude, depth
separated by :param:separator each.
'''
with open(filename, 'r') as infile:
for line in infile.readlines():
eventID, date, time, mag, lat, lon, depth = line.split(separator)[:7]
# skip first line
try:
day, month, year = date.split('/')
except:
continue
year = int(year)
# hardcoded, if year only consists of 2 digits (e.g. 16 instead of 2016)
if year < 100:
year += 2000
datetime = '{}-{}-{}T{}'.format(year, month, day, time)
try:
datetime = UTCDateTime(datetime)
except Exception as e:
print(e, datetime, filename)
continue
for event in self.eventlist:
if eventID in str(event.resource_id) or eventID in event.origins:
if event.origins:
origin = event.origins[0] # should have only one origin
if origin.time == datetime:
origin.latitude = float(lat)
origin.longitude = float(lon)
origin.depth = float(depth)
else:
continue
elif not event.origins:
origin = Origin(resource_id=event.resource_id,
time=datetime, latitude=float(lat),
longitude=float(lon), depth=float(depth))
event.origins.append(origin)
event.magnitudes.append(Magnitude(resource_id=event.resource_id,
mag=float(mag),
mag_type='M'))
break
def search_eventfile_info(self):
'''
Search all datapaths in rootpath for filenames with given file extension fext
and try to read event info from it
'''
datapaths = []
fext = '.csv'
for event in self.eventlist:
if not event.datapath in datapaths:
datapaths.append(event.datapath)
for datapath in datapaths:
# datapath = os.path.join(self.rootpath, datapath)
if os.path.isdir(datapath):
for filename in os.listdir(datapath):
filename = os.path.join(datapath, filename)
if os.path.isfile(filename) and filename.endswith(fext):
try:
self.read_eventfile_info(filename)
except Exception as e:
print('Failed on reading eventfile info from file {}: {}'.format(filename, e))
else:
print("Directory %s does not exist!" % datapath)
def getPaths(self):
'''
Returns paths (eventlist) of all events saved in the project.
'''
paths = []
for event in self.eventlist:
paths.append(event.path)
return paths
def setDirty(self, value=True):
self.dirty = value
def getEventFromPath(self, path):
'''
Search for an event in the project by event path.
'''
for event in self.eventlist:
if event.path == path:
return event
def save(self, filename=None):
'''
Save PyLoT Project to a file.
Can be loaded by using project.load(filename).
'''
try:
import pickle
except ImportError:
import _pickle as pickle
if filename:
self.location = filename
else:
filename = self.location
table = self._table # MP: see below
try:
outfile = open(filename, 'wb')
self._table = [] # MP: Workaround as long as table cannot be saved as part of project
pickle.dump(self, outfile, protocol=pickle.HIGHEST_PROTOCOL)
self.setDirty(False)
self._table = table # MP: see above
return True
except Exception as e:
print('Could not pickle PyLoT project. Reason: {}'.format(e))
self.setDirty()
self._table = table # MP: see above
return False
@staticmethod
def load(filename):
'''
Load project from filename.
'''
import pickle
infile = open(filename, 'rb')
project = pickle.load(infile)
infile.close()
project.location = filename
print('Loaded %s' % filename)
return project

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)]

123
pylot/core/io/waveformdata.py Normal file
View File

@ -0,0 +1,123 @@
import logging
from dataclasses import dataclass, field
from typing import Union, List
from obspy import Stream, read
from obspy.io.sac import SacIOError
from pylot.core.io.utils import validate_filenames
from pylot.core.util.dataprocessing import Metadata
from pylot.core.util.utils import get_stations, check_for_nan, check4rotated
@dataclass
class WaveformData:
wfdata: Stream = field(default_factory=Stream)
wforiginal: Union[Stream, None] = None
wf_alt: Stream = field(default_factory=Stream)
dirty: bool = False
def load_waveforms(self, fnames: List[str], fnames_alt: List[str] = None, check_rotated=False, metadata=None, tstart=0, tstop=0):
fn_list = validate_filenames(fnames)
if not fn_list:
logging.warning('No valid filenames given for loading waveforms')
else:
self.clear()
self.add_waveforms(fn_list)
if fnames_alt is None:
pass
else:
alt_fn_list = validate_filenames(fnames_alt)
if not alt_fn_list:
logging.warning('No valid alternative filenames given for loading waveforms')
else:
self.add_waveforms(alt_fn_list, alternative=True)
if not fn_list and not alt_fn_list:
logging.error('No filenames or alternative filenames given for loading waveforms')
return False
self.merge()
self.replace_nan()
if not check_rotated or not metadata:
pass
else:
self.rotate_zne()
self.trim_station_traces()
self.wforiginal = self.wfdata.copy()
self.dirty = False
return True
def add_waveforms(self, fnames: List[str], alternative: bool = False):
data_stream = self.wf_alt if alternative else self.wfdata
warnmsg = ''
for fname in set(fnames):
try:
data_stream += read(fname)
except TypeError:
try:
data_stream += read(fname, format='GSE2')
except Exception as e:
try:
data_stream += read(fname, format='SEGY')
except Exception as e:
warnmsg += f'{fname}\n{e}\n'
except SacIOError as se:
warnmsg += f'{fname}\n{se}\n'
if warnmsg:
print(f'WARNING in add_waveforms: unable to read waveform data\n{warnmsg}')
def clear(self):
self.wfdata = Stream()
self.wforiginal = None
self.wf_alt = Stream()
def reset(self):
"""
Resets the waveform data to its original state.
"""
if self.wforiginal:
self.wfdata = self.wforiginal.copy()
else:
self.wfdata = Stream()
self.dirty = False
def merge(self):
"""
check for gaps in Stream and merge if gaps are found
"""
gaps = self.wfdata.get_gaps()
if gaps:
merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
self.wfdata.merge(method=1)
logging.info('Merged the following stations because of gaps:')
for station in merged:
logging.info(station)
def replace_nan(self):
"""
Replace all NaNs in data with 0. (in place)
"""
self.wfdata = check_for_nan(self.wfdata)
def rotate_zne(self, metadata: Metadata = None):
"""
Check all traces in stream for rotation. If a trace is not in ZNE rotation (last symbol of channel code is numeric) and the trace
is in the metadata with azimuth and dip, rotate it to classical ZNE orientation.
Rotating the traces requires them to be of the same length, so, all traces will be trimmed to a common length as a
side effect.
"""
self.wfdata = check4rotated(self.wfdata, metadata)
def trim_station_traces(self):
"""
trim data stream to common time window
"""
for station in get_stations(self.wfdata):
station_traces = self.wfdata.select(station=station)
station_traces.trim(starttime=max([trace.stats.starttime for trace in station_traces]),
endtime=min([trace.stats.endtime for trace in station_traces]))

4
pylot/core/loc/focmec.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString() __version__ = _getVersionString()
@ -25,4 +25,4 @@ def export(picks, fnout, parameter, eventinfo):
:type eventinfo: list object :type eventinfo: list object
''' '''
# write phases to FOCMEC-phase file # write phases to FOCMEC-phase file
writephases(picks, 'FOCMEC', fnout, parameter, eventinfo) write_phases(picks, 'FOCMEC', fnout, parameter, eventinfo)

4
pylot/core/loc/hash.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString() __version__ = _getVersionString()
@ -25,4 +25,4 @@ def export(picks, fnout, parameter, eventinfo):
:type eventinfo: list object :type eventinfo: list object
''' '''
# write phases to HASH-phase file # write phases to HASH-phase file
writephases(picks, 'HASH', fnout, parameter, eventinfo) write_phases(picks, 'HASH', fnout, parameter, eventinfo)

4
pylot/core/loc/hypo71.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString() __version__ = _getVersionString()
@ -22,4 +22,4 @@ def export(picks, fnout, parameter):
:type parameter: object :type parameter: object
''' '''
# write phases to HYPO71-phase file # write phases to HYPO71-phase file
writephases(picks, 'HYPO71', fnout, parameter) write_phases(picks, 'HYPO71', fnout, parameter)

4
pylot/core/loc/hypodd.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString() __version__ = _getVersionString()
@ -25,4 +25,4 @@ def export(picks, fnout, parameter, eventinfo):
:type eventinfo: list object :type eventinfo: list object
''' '''
# write phases to hypoDD-phase file # write phases to hypoDD-phase file
writephases(picks, 'HYPODD', fnout, parameter, eventinfo) write_phases(picks, 'HYPODD', fnout, parameter, eventinfo)

4
pylot/core/loc/hyposat.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString() __version__ = _getVersionString()
@ -22,4 +22,4 @@ def export(picks, fnout, parameter):
:type parameter: object :type parameter: object
''' '''
# write phases to HYPOSAT-phase file # write phases to HYPOSAT-phase file
writephases(picks, 'HYPOSAT', fnout, parameter) write_phases(picks, 'HYPOSAT', fnout, parameter)

4
pylot/core/loc/nll.py
View File

@ -7,7 +7,7 @@ import subprocess
from obspy import read_events from obspy import read_events
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.gui import which from pylot.core.util.gui import which
from pylot.core.util.utils import getPatternLine, runProgram from pylot.core.util.utils import getPatternLine, runProgram
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
@ -34,7 +34,7 @@ def export(picks, fnout, parameter):
:type parameter: object :type parameter: object
''' '''
# write phases to NLLoc-phase file # write phases to NLLoc-phase file
writephases(picks, 'NLLoc', fnout, parameter) write_phases(picks, 'NLLoc', fnout, parameter)
def modify_inputs(ctrfn, root, nllocoutn, phasefn, tttn): def modify_inputs(ctrfn, root, nllocoutn, phasefn, tttn):

4
pylot/core/loc/velest.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pylot.core.io.phases import writephases from pylot.core.io.phases import write_phases
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString() __version__ = _getVersionString()
@ -25,4 +25,4 @@ def export(picks, fnout, eventinfo, parameter=None):
:type parameter: object :type parameter: object
''' '''
# write phases to VELEST-phase file # write phases to VELEST-phase file
writephases(picks, 'VELEST', fnout, parameter, eventinfo) write_phases(picks, 'VELEST', fnout, parameter, eventinfo)

4
pylot/core/pick/autopick.py
View File

@ -22,7 +22,7 @@ from pylot.core.pick.picker import AICPicker, PragPicker
from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \ from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
getSNR, fmpicker, checkPonsets, wadaticheck, get_quality_class, PickingFailedException, MissingTraceException getSNR, fmpicker, checkPonsets, wadaticheck, get_quality_class, PickingFailedException, MissingTraceException
from pylot.core.util.utils import getPatternLine, gen_Pool, \ from pylot.core.util.utils import getPatternLine, gen_Pool, \
get_bool, identifyPhaseID, get_None, correct_iplot get_bool, identifyPhaseID, get_none, correct_iplot
def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None, ncores=0, metadata=None, origin=None): def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None, ncores=0, metadata=None, origin=None):
@ -258,7 +258,7 @@ class AutopickStation(object):
self.pickparams = copy.deepcopy(pickparam) self.pickparams = copy.deepcopy(pickparam)
self.verbose = verbose self.verbose = verbose
self.iplot = correct_iplot(iplot) self.iplot = correct_iplot(iplot)
self.fig_dict = get_None(fig_dict) self.fig_dict = get_none(fig_dict)
self.metadata = metadata self.metadata = metadata
self.origin = origin self.origin = origin

18
pylot/core/pick/utils.py
View File

@ -15,7 +15,7 @@ import numpy as np
from obspy.core import Stream, UTCDateTime from obspy.core import Stream, UTCDateTime
from scipy.signal import argrelmax from scipy.signal import argrelmax
from pylot.core.util.utils import get_bool, get_None, SetChannelComponents from pylot.core.util.utils import get_bool, get_none, SetChannelComponents
def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecolor='k'): def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecolor='k'):
@ -136,7 +136,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol
PickError = symmetrize_error(diffti_te, diffti_tl) PickError = symmetrize_error(diffti_te, diffti_tl)
if iplot > 1: if iplot > 1:
if get_None(fig) is None: if get_none(fig) is None:
fig = plt.figure() # iplot) fig = plt.figure() # iplot)
plt_flag = 1 plt_flag = 1
fig._tight = True fig._tight = True
@ -275,7 +275,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
try: try:
P1 = np.polyfit(xslope1, xraw[islope1], 1) P1 = np.polyfit(xslope1, xraw[islope1], 1)
datafit1 = np.polyval(P1, xslope1) datafit1 = np.polyval(P1, xslope1)
except Exception as e: except ValueError as e:
print("fmpicker: Problems with data fit! {}".format(e)) print("fmpicker: Problems with data fit! {}".format(e))
print("Skip first motion determination!") print("Skip first motion determination!")
return FM return FM
@ -321,7 +321,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
try: try:
P2 = np.polyfit(xslope2, xfilt[islope2], 1) P2 = np.polyfit(xslope2, xfilt[islope2], 1)
datafit2 = np.polyval(P2, xslope2) datafit2 = np.polyval(P2, xslope2)
except Exception as e: except ValueError as e:
emsg = 'fmpicker: polyfit failed: {}'.format(e) emsg = 'fmpicker: polyfit failed: {}'.format(e)
print(emsg) print(emsg)
return FM return FM
@ -344,7 +344,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
print("fmpicker: Found polarity %s" % FM) print("fmpicker: Found polarity %s" % FM)
if iplot > 1: if iplot > 1:
if get_None(fig) is None: if get_none(fig) is None:
fig = plt.figure() # iplot) fig = plt.figure() # iplot)
plt_flag = 1 plt_flag = 1
fig._tight = True fig._tight = True
@ -868,7 +868,7 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
returnflag = 0 returnflag = 0
if iplot > 1: if iplot > 1:
if get_None(fig) is None: if get_none(fig) is None:
fig = plt.figure() # iplot) fig = plt.figure() # iplot)
plt_flag = 1 plt_flag = 1
fig._tight = True fig._tight = True
@ -1213,14 +1213,14 @@ def checkZ4S(X, pick, pickparams, iplot, fig=None, linecolor='k'):
t = np.linspace(diff_dict[key], trace.stats.endtime - trace.stats.starttime + diff_dict[key], t = np.linspace(diff_dict[key], trace.stats.endtime - trace.stats.starttime + diff_dict[key],
trace.stats.npts) trace.stats.npts)
if i == 0: if i == 0:
if get_None(fig) is None: if get_none(fig) is None:
fig = plt.figure() # self.iplot) ### WHY? MP MP fig = plt.figure() # self.iplot) ### WHY? MP MP
plt_flag = 1 plt_flag = 1
ax1 = fig.add_subplot(3, 1, i + 1) ax1 = fig.add_subplot(3, 1, i + 1)
ax = ax1 ax = ax1
ax.set_title('CheckZ4S, Station %s' % zdat[0].stats.station) ax.set_title('CheckZ4S, Station %s' % zdat[0].stats.station)
else: else:
if get_None(fig) is None: if get_none(fig) is None:
fig = plt.figure() # self.iplot) ### WHY? MP MP fig = plt.figure() # self.iplot) ### WHY? MP MP
plt_flag = 1 plt_flag = 1
ax = fig.add_subplot(3, 1, i + 1, sharex=ax1) ax = fig.add_subplot(3, 1, i + 1, sharex=ax1)
@ -1494,7 +1494,7 @@ def getQualityFromUncertainty(uncertainty, Errors):
# set initial quality to 4 (worst) and change only if one condition is hit # set initial quality to 4 (worst) and change only if one condition is hit
quality = 4 quality = 4
if get_None(uncertainty) is None: if get_none(uncertainty) is None:
return quality return quality
if uncertainty <= Errors[0]: if uncertainty <= Errors[0]:

14
pylot/core/util/array_map.py
View File

@ -124,8 +124,8 @@ class Array_map(QtWidgets.QWidget):
self.cmaps_box = QtWidgets.QComboBox() self.cmaps_box = QtWidgets.QComboBox()
self.cmaps_box.setMaxVisibleItems(20) self.cmaps_box.setMaxVisibleItems(20)
[self.cmaps_box.addItem(map_name) for map_name in sorted(plt.colormaps())] [self.cmaps_box.addItem(map_name) for map_name in sorted(plt.colormaps())]
# try to set to hsv as default # try to set to viridis as default
self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('hsv')) self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('viridis'))
self.top_row.addWidget(QtWidgets.QLabel('Select a phase: ')) self.top_row.addWidget(QtWidgets.QLabel('Select a phase: '))
self.top_row.addWidget(self.comboBox_phase) self.top_row.addWidget(self.comboBox_phase)
@ -474,17 +474,19 @@ class Array_map(QtWidgets.QWidget):
transform=ccrs.PlateCarree(), label='deleted')) transform=ccrs.PlateCarree(), label='deleted'))
def openPickDlg(self, ind): def openPickDlg(self, ind):
data = self._parent.get_data().getWFData() wfdata = self._parent.get_data().get_wf_data()
wfdata_comp = self._parent.get_data().get_wf_dataComp()
for index in ind: for index in ind:
network, station = self._station_onpick_ids[index].split('.')[:2] network, station = self._station_onpick_ids[index].split('.')[:2]
pyl_mw = self._parent pyl_mw = self._parent
try: try:
data = data.select(station=station) wfdata = wfdata.select(station=station)
if not data: wfdata_comp = wfdata_comp.select(station=station)
if not wfdata:
self._warn('No data for station {}'.format(station)) self._warn('No data for station {}'.format(station))
return return
pickDlg = PickDlg(self._parent, parameter=self.parameter, pickDlg = PickDlg(self._parent, parameter=self.parameter,
data=data, network=network, station=station, data=wfdata.copy(), data_compare=wfdata_comp.copy(), network=network, station=station,
picks=self._parent.get_current_event().getPick(station), picks=self._parent.get_current_event().getPick(station),
autopicks=self._parent.get_current_event().getAutopick(station), autopicks=self._parent.get_current_event().getAutopick(station),
filteroptions=self._parent.filteroptions, metadata=self.metadata, filteroptions=self._parent.filteroptions, metadata=self.metadata,

174
pylot/core/util/dataprocessing.py
View File

@ -338,19 +338,6 @@ class Metadata(object):
return inv, exc return inv, exc
def time_from_header(header):
"""
Function takes in the second line from a .gse file and takes out the date and time from that line.
:param header: second line from .gse file
:type header: string
:return: a list of integers of form [year, month, day, hour, minute, second, microsecond]
"""
timeline = header.split(' ')
time = timeline[1].split('/') + timeline[2].split(':')
time = time[:-1] + time[-1].split('.')
return [int(t) for t in time]
def check_time(datetime): def check_time(datetime):
""" """
Function takes in date and time as list and validates it's values by trying to make an UTCDateTime object from it Function takes in date and time as list and validates it's values by trying to make an UTCDateTime object from it
@ -388,167 +375,6 @@ def check_time(datetime):
return False return False
# TODO: change root to datapath
def get_file_list(root_dir):
"""
Function uses a directorie to get all the *.gse files from it.
:param root_dir: a directorie leading to the .gse files
:type root_dir: string
:return: returns a list of filenames (without path to them)
"""
file_list = glob.glob1(root_dir, '*.gse')
return file_list
def checks_station_second(datetime, file):
"""
Function uses the given list to check if the parameter 'second' is set to 60 by mistake
and sets the time correctly if so. Can only correct time if no date change would be necessary.
:param datetime: [year, month, day, hour, minute, second, microsecond]
:return: returns the input with the correct value for second
"""
if datetime[5] == 60:
if datetime[4] == 59:
if datetime[3] == 23:
err_msg = 'Date should be next day. ' \
'File not changed: {0}'.format(file)
raise ValueError(err_msg)
else:
datetime[3] += 1
datetime[4] = 0
datetime[5] = 0
else:
datetime[4] += 1
datetime[5] = 0
return datetime
def make_time_line(line, datetime):
"""
Function takes in the original line from a .gse file and a list of date and
time values to make a new line with corrected date and time.
:param line: second line from .gse file.
:type line: string
:param datetime: list of integers [year, month, day, hour, minute, second, microsecond]
:type datetime: list
:return: returns a string to write it into a file.
"""
ins_form = '{0:02d}:{1:02d}:{2:02d}.{3:03d}'
insertion = ins_form.format(int(datetime[3]),
int(datetime[4]),
int(datetime[5]),
int(datetime[6] * 1e-3))
newline = line[:16] + insertion + line[28:]
return newline
def evt_head_check(root_dir, out_dir=None):
"""
A function to make sure that an arbitrary number of .gse files have correct values in their header.
:param root_dir: a directory leading to the .gse files.
:type root_dir: string
:param out_dir: a directory to store the new files somwhere els.
:return: returns nothing
"""
if not out_dir:
print('WARNING files are going to be overwritten!')
inp = str(input('Continue? [y/N]'))
if not inp == 'y':
sys.exit()
filelist = get_file_list(root_dir)
nfiles = 0
for file in filelist:
infile = open(os.path.join(root_dir, file), 'r')
lines = infile.readlines()
infile.close()
datetime = time_from_header(lines[1])
if check_time(datetime):
continue
else:
nfiles += 1
datetime = checks_station_second(datetime, file)
print('writing ' + file)
# write File
lines[1] = make_time_line(lines[1], datetime)
if not out_dir:
out = open(os.path.join(root_dir, file), 'w')
out.writelines(lines)
out.close()
else:
out = open(os.path.join(out_dir, file), 'w')
out.writelines(lines)
out.close()
print(nfiles)
def read_metadata(path_to_inventory):
"""
take path_to_inventory and return either the corresponding list of files
found or the Parser object for a network dataless seed volume to prevent
read overhead for large dataless seed volumes
:param path_to_inventory:
:return: tuple containing a either list of files or `obspy.io.xseed.Parser`
object and the inventory type found
:rtype: tuple
"""
dlfile = list()
invfile = list()
respfile = list()
# possible file extensions specified here:
inv = dict(dless=dlfile, xml=invfile, resp=respfile, dseed=dlfile[:])
if os.path.isfile(path_to_inventory):
ext = os.path.splitext(path_to_inventory)[1].split('.')[1]
inv[ext] += [path_to_inventory]
else:
for ext in inv.keys():
inv[ext] += glob.glob1(path_to_inventory, '*.{0}'.format(ext))
invtype = key_for_set_value(inv)
if invtype is None:
print("Neither dataless-SEED file, inventory-xml file nor "
"RESP-file found!")
print("!!WRONG CALCULATION OF SOURCE PARAMETERS!!")
robj = None,
elif invtype == 'dless': # prevent multiple read of large dlsv
print("Reading metadata information from dataless-SEED file ...")
if len(inv[invtype]) == 1:
fullpath_inv = os.path.join(path_to_inventory, inv[invtype][0])
robj = Parser(fullpath_inv)
else:
robj = inv[invtype]
else:
print("Reading metadata information from inventory-xml file ...")
robj = read_inventory(inv[invtype])
return invtype, robj
# idea to optimize read_metadata
# def read_metadata_new(path_to_inventory):
# metadata_objects = []
# # read multiple files from directory
# if os.path.isdir(path_to_inventory):
# fnames = os.listdir(path_to_inventory)
# # read single file
# elif os.path.isfile(path_to_inventory):
# fnames = [path_to_inventory]
# else:
# print("Neither dataless-SEED file, inventory-xml file nor "
# "RESP-file found!")
# print("!!WRONG CALCULATION OF SOURCE PARAMETERS!!")
# fnames = []
#
# for fname in fnames:
# path_to_inventory_filename = os.path.join(path_to_inventory, fname)
# try:
# ftype, robj = read_metadata_file(path_to_inventory_filename)
# metadata_objects.append((ftype, robj))
# except Exception as e:
# print('Could not read metadata file {} '
# 'because of the following Exception: {}'.format(path_to_inventory_filename, e))
# return metadata_objects
def restitute_trace(input_tuple): def restitute_trace(input_tuple):
def no_metadata(tr, seed_id): def no_metadata(tr, seed_id):
print('no metadata file found ' print('no metadata file found '

2
pylot/core/util/generate_array_maps.py
View File

@ -20,7 +20,7 @@ import matplotlib
matplotlib.use('Qt5Agg') matplotlib.use('Qt5Agg')
sys.path.append(os.path.join('/'.join(sys.argv[0].split('/')[:-1]), '../../..')) sys.path.append(os.path.join('/'.join(sys.argv[0].split('/')[:-1]), '../../..'))
from PyLoT import Project from pylot.core.io.project import Project
from pylot.core.util.dataprocessing import Metadata from pylot.core.util.dataprocessing import Metadata
from pylot.core.util.array_map import Array_map from pylot.core.util.array_map import Array_map

4
pylot/core/util/structure.py
View File

@ -6,7 +6,7 @@ Created on Wed Jan 26 17:47:25 2015
@author: sebastianw @author: sebastianw
""" """
from pylot.core.io.data import SeiscompDataStructure, PilotDataStructure, ObspyDMTdataStructure from pylot.core.io.data import SeiscompDataStructure, PilotDataStructure
DATASTRUCTURE = {'PILOT': PilotDataStructure, 'SeisComP': SeiscompDataStructure, DATASTRUCTURE = {'PILOT': PilotDataStructure, 'SeisComP': SeiscompDataStructure,
'obspyDMT': ObspyDMTdataStructure, None: PilotDataStructure} 'obspyDMT': PilotDataStructure, None: PilotDataStructure}

388
pylot/core/util/utils.py
View File

@ -20,6 +20,10 @@ from pylot.core.io.inputs import PylotParameter, FilterOptions
from pylot.core.util.obspyDMT_interface import check_obspydmt_eventfolder from pylot.core.util.obspyDMT_interface import check_obspydmt_eventfolder
from pylot.styles import style_settings from pylot.styles import style_settings
Rgba: Type[tuple] = Tuple[int, int, int, int]
Mplrgba: Type[tuple] = Tuple[float, float, float, float]
Mplrgbastr: Type[tuple] = Tuple[str, str, str, str]
def _pickle_method(m): def _pickle_method(m):
if m.im_self is None: if m.im_self is None:
@ -83,25 +87,6 @@ def fit_curve(x, y):
return splev, splrep(x, y) return splev, splrep(x, y)
def getindexbounds(f, eta):
"""
Get indices of values closest below and above maximum value in an array
:param f: array
:type f: `~numpy.ndarray`
:param eta: look for value in array that is closes to max_value * eta
:type eta: float
:return: tuple containing index of max value, index of value closest below max value,
index of value closest above max value
:rtype: (int, int, int)
"""
mi = f.argmax() # get indices of max values
m = max(f) # get maximum value
b = m * eta #
l = find_nearest(f[:mi], b) # find closest value below max value
u = find_nearest(f[mi:], b) + mi # find closest value above max value
return mi, l, u
def gen_Pool(ncores=0): def gen_Pool(ncores=0):
""" """
Generate mulitprocessing pool object utilizing ncores amount of cores Generate mulitprocessing pool object utilizing ncores amount of cores
@ -121,7 +106,7 @@ def gen_Pool(ncores=0):
print('gen_Pool: Generated multiprocessing Pool with {} cores\n'.format(ncores)) print('gen_Pool: Generated multiprocessing Pool with {} cores\n'.format(ncores))
pool = multiprocessing.Pool(ncores, maxtasksperchild=100) pool = multiprocessing.Pool(ncores)
return pool return pool
@ -167,11 +152,11 @@ def clims(lim1, lim2):
""" """
takes two pairs of limits and returns one pair of common limts takes two pairs of limits and returns one pair of common limts
:param lim1: limit 1 :param lim1: limit 1
:type lim1: int :type lim1: List[int]
:param lim2: limit 2 :param lim2: limit 2
:type lim2: int :type lim2: List[int]
:return: new upper and lower limit common to both given limits :return: new upper and lower limit common to both given limits
:rtype: [int, int] :rtype: List[int]
>>> clims([0, 4], [1, 3]) >>> clims([0, 4], [1, 3])
[0, 4] [0, 4]
@ -303,7 +288,7 @@ def fnConstructor(s):
if type(s) is str: if type(s) is str:
s = s.split(':')[-1] s = s.split(':')[-1]
else: else:
s = getHash(UTCDateTime()) s = get_hash(UTCDateTime())
badchars = re.compile(r'[^A-Za-z0-9_. ]+|^\.|\.$|^ | $|^$') badchars = re.compile(r'[^A-Za-z0-9_. ]+|^\.|\.$|^ | $|^$')
badsuffix = re.compile(r'(aux|com[1-9]|con|lpt[1-9]|prn)(\.|$)') badsuffix = re.compile(r'(aux|com[1-9]|con|lpt[1-9]|prn)(\.|$)')
@ -315,15 +300,32 @@ def fnConstructor(s):
return fn return fn
def get_None(value): def get_none(value):
""" """
Convert "None" to None Convert "None" to None
:param value: :param value:
:type value: str, bool :type value: str, NoneType
:return: :return:
:rtype: bool :rtype: type(value) or NoneType
>>> st = read()
>>> print(get_none(st))
3 Trace(s) in Stream:
BW.RJOB..EHZ | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples
BW.RJOB..EHN | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples
BW.RJOB..EHE | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples
>>> get_none('Stream')
'Stream'
>>> get_none(0)
0
>>> get_none(0.)
0.0
>>> print(get_none('None'))
None
>>> print(get_none(None))
None
""" """
if value == 'None': if value is None or (type(value) is str and value == 'None'):
return None return None
else: else:
return value return value
@ -331,20 +333,47 @@ def get_None(value):
def get_bool(value): def get_bool(value):
""" """
Convert string representations of bools to their true boolean value Convert string representations of bools to their true boolean value. Return value if it cannot be identified as bool.
:param value: :param value:
:type value: str, bool :type value: str, bool, int, float
:return: true boolean value :return: true boolean value
:rtype: bool :rtype: bool
>>> get_bool(True)
True
>>> get_bool(False)
False
>>> get_bool(0)
False
>>> get_bool(0.)
False
>>> get_bool(0.1)
True
>>> get_bool(2)
True
>>> get_bool(-1)
False
>>> get_bool(-0.3)
False
>>> get_bool(None)
None
>>> get_bool('Stream')
'Stream'
""" """
if type(value) is bool: if type(value) == bool:
return value return value
elif value in ['True', 'true']: elif value in ['True', 'true']:
return True return True
elif value in ['False', 'false']: elif value in ['False', 'false']:
return False return False
elif isinstance(value, float) or isinstance(value, int):
if value > 0. or value > 0:
return True
else: else:
return bool(value) return False
else:
return value
def four_digits(year): def four_digits(year):
""" """
@ -355,8 +384,8 @@ def four_digits(year):
:return: four digit year correspondent :return: four digit year correspondent
:rtype: int :rtype: int
>>> four_digits(20) >>> four_digits(75)
1920 1975
>>> four_digits(16) >>> four_digits(16)
2016 2016
>>> four_digits(00) >>> four_digits(00)
@ -438,36 +467,53 @@ def backtransformFilterString(st):
return st return st
def getHash(time): def get_hash(time):
""" """
takes a time object and returns the corresponding SHA1 hash of the formatted date string takes a time object and returns the corresponding SHA1 hash of the formatted date string
:param time: time object for which a hash should be calculated :param time: time object for which a hash should be calculated
:type time: `~obspy.core.utcdatetime.UTCDateTime` :type time: `~obspy.core.utcdatetime.UTCDateTime`
:return: SHA1 hash :return: SHA1 hash
:rtype: str :rtype: str
>>> time = UTCDateTime(0)
>>> get_hash(time)
'7627cce3b1b58dd21b005dac008b34d18317dd15'
>>> get_hash(0)
Traceback (most recent call last):
...
AssertionError: 'time' is not an ObsPy UTCDateTime object
""" """
assert isinstance(time, UTCDateTime), '\'time\' is not an ObsPy UTCDateTime object'
hg = hashlib.sha1() hg = hashlib.sha1()
hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f')) hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f').encode('utf-8'))
return hg.hexdigest() return hg.hexdigest()
def getLogin(): def get_login():
""" """
returns the actual user's login ID returns the actual user's name
:return: login ID :return: login name
:rtype: str :rtype: str
""" """
import getpass import getpass
return getpass.getuser() return getpass.getuser()
def getOwner(fn): def get_owner(fn):
""" """
takes a filename and return the login ID of the actual owner of the file takes a filename and return the login ID of the actual owner of the file
:param fn: filename of the file tested :param fn: filename of the file tested
:type fn: str :type fn: str
:return: login ID of the file's owner :return: login ID of the file's owner
:rtype: str :rtype: str
>>> import tempfile
>>> with tempfile.NamedTemporaryFile() as tmpfile:
... tmpfile.write(b'') and True
... tmpfile.flush()
... get_owner(tmpfile.name) == os.path.expanduser('~').split('/')[-1]
0
True
""" """
system_name = platform.system() system_name = platform.system()
if system_name in ["Linux", "Darwin"]: if system_name in ["Linux", "Darwin"]:
@ -513,6 +559,11 @@ def is_executable(fn):
:param fn: path to the file to be tested :param fn: path to the file to be tested
:return: True or False :return: True or False
:rtype: bool :rtype: bool
>>> is_executable('/bin/ls')
True
>>> is_executable('/var/log/system.log')
False
""" """
return os.path.isfile(fn) and os.access(fn, os.X_OK) return os.path.isfile(fn) and os.access(fn, os.X_OK)
@ -539,24 +590,36 @@ def isSorted(iterable):
>>> isSorted([2,3,1,4]) >>> isSorted([2,3,1,4])
False False
""" """
assert isIterable(iterable), 'object is not iterable; object: {' \ assert is_iterable(iterable), "object is not iterable; object: {}".format(iterable)
'}'.format(iterable)
if type(iterable) is str: if type(iterable) is str:
iterable = [s for s in iterable] iterable = [s for s in iterable]
return sorted(iterable) == iterable return sorted(iterable) == iterable
def isIterable(obj): def is_iterable(obj):
""" """
takes a python object and returns True is the object is iterable and takes a python object and returns True is the object is iterable and
False otherwise False otherwise
:param obj: a python object :param obj: a python object
:type obj: object :type obj: obj
:return: True of False :return: True of False
:rtype: bool :rtype: bool
>>> is_iterable(1)
False
>>> is_iterable(True)
False
>>> is_iterable(0.)
False
>>> is_iterable((0,1,3,4))
True
>>> is_iterable([1])
True
>>> is_iterable('a')
True
""" """
try: try:
iterator = iter(obj) iter(obj)
except TypeError as te: except TypeError as te:
return False return False
return True return True
@ -565,13 +628,19 @@ def isIterable(obj):
def key_for_set_value(d): def key_for_set_value(d):
""" """
takes a dictionary and returns the first key for which's value the takes a dictionary and returns the first key for which's value the
boolean is True boolean representation is True
:param d: dictionary containing values :param d: dictionary containing values
:type d: dict :type d: dict
:return: key to the first non-False value found; None if no value's :return: key to the first non-False value found; None if no value's
boolean equals True boolean equals True
:rtype: :rtype: bool or NoneType
>>> key_for_set_value({'one': 0, 'two': 1})
'two'
>>> print(key_for_set_value({1: 0, 2: False}))
None
""" """
assert type(d) is dict, "Function only defined for inputs of type 'dict'."
r = None r = None
for k, v in d.items(): for k, v in d.items():
if v: if v:
@ -579,32 +648,53 @@ def key_for_set_value(d):
return r return r
def prepTimeAxis(stime, trace, verbosity=0): def prep_time_axis(offset, trace, verbosity=0):
""" """
takes a starttime and a trace object and returns a valid time axis for takes an offset and a trace object and returns a valid time axis for
plotting plotting
:param stime: start time of the actual seismogram as UTCDateTime :param offset: offset of the actual seismogram on plotting axis
:type stime: `~obspy.core.utcdatetime.UTCDateTime` :type offset: float or int
:param trace: seismic trace object :param trace: seismic trace object
:type trace: `~obspy.core.trace.Trace` :type trace: `~obspy.core.trace.Trace`
:param verbosity: if != 0, debug output will be written to console :param verbosity: if != 0, debug output will be written to console
:type verbosity: int :type verbosity: int
:return: valid numpy array with time stamps for plotting :return: valid numpy array with time stamps for plotting
:rtype: `~numpy.ndarray` :rtype: `~numpy.ndarray`
>>> tr = read()[0]
>>> prep_time_axis(0., tr)
array([0.00000000e+00, 1.00033344e-02, 2.00066689e-02, ...,
2.99799933e+01, 2.99899967e+01, 3.00000000e+01])
>>> prep_time_axis(22.5, tr)
array([22.5 , 22.51000333, 22.52000667, ..., 52.47999333,
52.48999667, 52.5 ])
>>> prep_time_axis(tr.stats.starttime, tr)
Traceback (most recent call last):
...
AssertionError: 'offset' is not of type 'float' or 'int'; type: <class 'obspy.core.utcdatetime.UTCDateTime'>
>>> tr.stats.npts -= 1
>>> prep_time_axis(0, tr)
array([0.00000000e+00, 1.00033356e-02, 2.00066711e-02, ...,
2.99699933e+01, 2.99799967e+01, 2.99900000e+01])
>>> tr.stats.npts += 2
>>> prep_time_axis(0, tr)
array([0.00000000e+00, 1.00033333e-02, 2.00066667e-02, ...,
2.99899933e+01, 2.99999967e+01, 3.00100000e+01])
""" """
assert isinstance(offset, (float, int)), "'offset' is not of type 'float' or 'int'; type: {}".format(type(offset))
nsamp = trace.stats.npts nsamp = trace.stats.npts
srate = trace.stats.sampling_rate srate = trace.stats.sampling_rate
tincr = trace.stats.delta tincr = trace.stats.delta
etime = stime + nsamp / srate etime = offset + nsamp / srate
time_ax = np.linspace(stime, etime, nsamp) time_ax = np.linspace(offset, etime, nsamp)
if len(time_ax) < nsamp: if len(time_ax) < nsamp:
if verbosity: if verbosity:
print('elongate time axes by one datum') print('elongate time axes by one datum')
time_ax = np.arange(stime, etime + tincr, tincr) time_ax = np.arange(offset, etime + tincr, tincr)
elif len(time_ax) > nsamp: elif len(time_ax) > nsamp:
if verbosity: if verbosity:
print('shorten time axes by one datum') print('shorten time axes by one datum')
time_ax = np.arange(stime, etime - tincr, tincr) time_ax = np.arange(offset, etime - tincr, tincr)
if len(time_ax) != nsamp: if len(time_ax) != nsamp:
print('Station {0}, {1} samples of data \n ' print('Station {0}, {1} samples of data \n '
'{2} length of time vector \n' '{2} length of time vector \n'
@ -620,13 +710,13 @@ def find_horizontals(data):
:param data: waveform data :param data: waveform data
:type data: `obspy.core.stream.Stream` :type data: `obspy.core.stream.Stream`
:return: components list :return: components list
:rtype: list :rtype: List(str)
..example:: ..example::
>>> st = read() >>> st = read()
>>> find_horizontals(st) >>> find_horizontals(st)
[u'N', u'E'] ['N', 'E']
""" """
rval = [] rval = []
for tr in data: for tr in data:
@ -637,7 +727,7 @@ def find_horizontals(data):
return rval return rval
def pick_color(picktype, phase, quality=0): def pick_color(picktype: Literal['manual', 'automatic'], phase: Literal['P', 'S'], quality: int = 0) -> Rgba:
""" """
Create pick color by modifying the base color by the quality. Create pick color by modifying the base color by the quality.
@ -650,7 +740,7 @@ def pick_color(picktype, phase, quality=0):
:param quality: quality of pick. Decides the new intensity of the modifier color :param quality: quality of pick. Decides the new intensity of the modifier color
:type quality: int :type quality: int
:return: tuple containing modified rgba color values :return: tuple containing modified rgba color values
:rtype: (int, int, int, int) :rtype: Rgba
""" """
min_quality = 3 min_quality = 3
bpc = base_phase_colors(picktype, phase) # returns dict like {'modifier': 'g', 'rgba': (0, 0, 255, 255)} bpc = base_phase_colors(picktype, phase) # returns dict like {'modifier': 'g', 'rgba': (0, 0, 255, 255)}
@ -706,17 +796,17 @@ def pick_linestyle_plt(picktype, key):
return linestyles[picktype][key] return linestyles[picktype][key]
def modify_rgba(rgba, modifier, intensity): def modify_rgba(rgba: Rgba, modifier: Literal['r', 'g', 'b'], intensity: float) -> Rgba:
""" """
Modify rgba color by adding the given intensity to the modifier color Modify rgba color by adding the given intensity to the modifier color
:param rgba: tuple containing rgba values :param rgba: tuple containing rgba values
:type rgba: (int, int, int, int) :type rgba: Rgba
:param modifier: which color should be modified, eg. 'r', 'g', 'b' :param modifier: which color should be modified; options: 'r', 'g', 'b'
:type modifier: str :type modifier: Literal['r', 'g', 'b']
:param intensity: intensity to be added to selected color :param intensity: intensity to be added to selected color
:type intensity: float :type intensity: float
:return: tuple containing rgba values :return: tuple containing rgba values
:rtype: (int, int, int, int) :rtype: Rgba
""" """
rgba = list(rgba) rgba = list(rgba)
index = {'r': 0, index = {'r': 0,
@ -750,18 +840,20 @@ def transform_colors_mpl_str(colors, no_alpha=False):
Transforms rgba color values to a matplotlib string of color values with a range of [0, 1] Transforms rgba color values to a matplotlib string of color values with a range of [0, 1]
:param colors: tuple of rgba color values ranging from [0, 255] :param colors: tuple of rgba color values ranging from [0, 255]
:type colors: (float, float, float, float) :type colors: (float, float, float, float)
:param no_alpha: Wether to return a alpha value in the matplotlib color string :param no_alpha: Whether to return an alpha value in the matplotlib color string
:type no_alpha: bool :type no_alpha: bool
:return: String containing r, g, b values and alpha value if no_alpha is False (default) :return: String containing r, g, b values and alpha value if no_alpha is False (default)
:rtype: str :rtype: str
>>> transform_colors_mpl_str((255., 255., 255., 255.), True)
'(1.0, 1.0, 1.0)'
>>> transform_colors_mpl_str((255., 255., 255., 255.))
'(1.0, 1.0, 1.0, 1.0)'
""" """
colors = list(colors)
colors_mpl = tuple([color / 255. for color in colors])
if no_alpha: if no_alpha:
colors_mpl = '({}, {}, {})'.format(*colors_mpl) return '({}, {}, {})'.format(*transform_colors_mpl(colors))
else: else:
colors_mpl = '({}, {}, {}, {})'.format(*colors_mpl) return '({}, {}, {}, {})'.format(*transform_colors_mpl(colors))
return colors_mpl
def transform_colors_mpl(colors): def transform_colors_mpl(colors):
@ -771,27 +863,16 @@ def transform_colors_mpl(colors):
:type colors: (float, float, float, float) :type colors: (float, float, float, float)
:return: tuple of rgba color values ranging from [0, 1] :return: tuple of rgba color values ranging from [0, 1]
:rtype: (float, float, float, float) :rtype: (float, float, float, float)
>>> transform_colors_mpl((127.5, 0., 63.75, 255.))
(0.5, 0.0, 0.25, 1.0)
>>> transform_colors_mpl(())
""" """
colors = list(colors) colors = list(colors)
colors_mpl = tuple([color / 255. for color in colors]) colors_mpl = tuple([color / 255. for color in colors])
return colors_mpl return colors_mpl
def remove_underscores(data):
"""
takes a `obspy.core.stream.Stream` object and removes all underscores
from station names
:param data: stream of seismic data
:type data: `~obspy.core.stream.Stream`
:return: data stream
:rtype: `~obspy.core.stream.Stream`
"""
# for tr in data:
# # remove underscores
# tr.stats.station = tr.stats.station.strip('_')
return data
def trim_station_components(data, trim_start=True, trim_end=True): def trim_station_components(data, trim_start=True, trim_end=True):
""" """
cut a stream so only the part common to all three traces is kept to avoid dealing with offsets cut a stream so only the part common to all three traces is kept to avoid dealing with offsets
@ -820,6 +901,19 @@ def trim_station_components(data, trim_start=True, trim_end=True):
return data return data
def merge_stream(stream):
gaps = stream.get_gaps()
if gaps:
# list of merged stations (seed_ids)
merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
stream.merge(method=1)
print('Merged the following stations because of gaps:')
for merged_station in merged:
print(merged_station)
return stream, gaps
def check4gapsAndRemove(data): def check4gapsAndRemove(data):
""" """
check for gaps in Stream and remove them check for gaps in Stream and remove them
@ -840,12 +934,12 @@ def check4gapsAndRemove(data):
return data return data
def check_for_gaps_and_merge(data): def check4gapsAndMerge(data):
""" """
check for gaps in Stream and merge if gaps are found check for gaps in Stream and merge if gaps are found
:param data: stream of seismic data :param data: stream of seismic data
:type data: `~obspy.core.stream.Stream` :type data: `~obspy.core.stream.Stream`
:return: data stream, gaps returned from obspy get_gaps :return: data stream
:rtype: `~obspy.core.stream.Stream` :rtype: `~obspy.core.stream.Stream`
""" """
gaps = data.get_gaps() gaps = data.get_gaps()
@ -856,7 +950,7 @@ def check_for_gaps_and_merge(data):
for merged_station in merged: for merged_station in merged:
print(merged_station) print(merged_station)
return data, gaps return data
def check4doubled(data): def check4doubled(data):
@ -928,11 +1022,11 @@ def get_possible_pylot_eventfile_extensions(event, fext):
def get_stations(data): def get_stations(data):
""" """
Get list of all station names in data stream Get list of all station names in data-stream
:param data: stream containing seismic traces :param data: stream containing seismic traces
:type data: `~obspy.core.stream.Stream` :type data: `~obspy.core.stream.Stream`
:return: list of all station names in data, no duplicates :return: list of all station names in data, no duplicates
:rtype: list of str :rtype: List(str)
""" """
stations = [] stations = []
for tr in data: for tr in data:
@ -959,66 +1053,87 @@ def check4rotated(data, metadata=None, verbosity=1):
:rtype: `~obspy.core.stream.Stream` :rtype: `~obspy.core.stream.Stream`
""" """
def rotate_components(wfstream, metadata=None): def rotation_required(trace_ids):
"""
Derive if any rotation is required from the orientation code of the input.
:param trace_ids: string identifier of waveform data trace
:type trace_ids: List(str)
:return: boolean representing if rotation is necessary for any of the traces
:rtype: bool
"""
orientations = [trace_id[-1] for trace_id in trace_ids]
return any([orientation.isnumeric() for orientation in orientations])
def rotate_components(wfs_in, metadata=None):
""" """
Rotate components if orientation code is numeric (= non traditional orientation). 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. 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 Returns unrotated traces of no metadata is provided
:param wfstream: stream containing seismic traces of a station :param wfs_in: stream containing seismic traces of a station
:type wfstream: `~obspy.core.stream.Stream` :type wfs_in: `~obspy.core.stream.Stream`
:param metadata: tuple containing metadata type string and metadata parser object :param metadata: tuple containing metadata type string and metadata parser object
:type metadata: (str, `~obspy.io.xseed.parser.Parser`) :type metadata: (str, `~obspy.io.xseed.parser.Parser`)
:return: stream object with traditionally oriented traces (ZNE) :return: stream object with traditionally oriented traces (ZNE)
:rtype: `~obspy.core.stream.Stream` :rtype: `~obspy.core.stream.Stream`
""" """
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 # check if any traces in this station need to be rotated
trace_ids = [trace.id for trace in wfstream] trace_ids = [trace.id for trace in wfs_in]
orientations = [trace_id[-1] for trace_id in trace_ids] if not rotation_required(trace_ids):
rotation_required = [orientation.isnumeric() for orientation in orientations] logging.debug(f"Stream does not need any rotation: Traces are {trace_ids=}")
if any(rotation_required): return wfs_in
t_start = full_range(wfstream)
# check metadata quality
t_start = full_range(wfs_in)
try: try:
azimuts = [] azimuths = []
dips = [] dips = []
for tr_id in trace_ids: for tr_id in trace_ids:
azimuts.append(metadata.get_coordinates(tr_id, t_start)['azimuth']) azimuths.append(metadata.get_coordinates(tr_id, t_start)['azimuth'])
dips.append(metadata.get_coordinates(tr_id, t_start)['dip']) dips.append(metadata.get_coordinates(tr_id, t_start)['dip'])
except (KeyError, TypeError) as e: except (KeyError, TypeError) as err:
print('Failed to rotate trace {}, no azimuth or dip available in metadata'.format(tr_id)) logging.error(f"{type(err)=} occurred: {err=} Rotating not possible, not all azimuth and dip information "
return wfstream f"available in metadata. Stream remains unchanged.")
if len(wfstream) < 3: return wfs_in
print('Failed to rotate Stream {}, not enough components available.'.format(wfstream)) except Exception as err:
return wfstream print(f"Unexpected {err=}, {type(err)=}")
raise
# to rotate all traces must have same length, so trim them # to rotate all traces must have same length, so trim them
wfstream = trim_station_components(wfstream, trim_start=True, trim_end=True) wfs_out = trim_station_components(wfs_in, trim_start=True, trim_end=True)
try: try:
z, n, e = rotate2zne(wfstream[0], azimuts[0], dips[0], z, n, e = rotate2zne(wfs_out[0], azimuths[0], dips[0],
wfstream[1], azimuts[1], dips[1], wfs_out[1], azimuths[1], dips[1],
wfstream[2], azimuts[2], dips[2]) wfs_out[2], azimuths[2], dips[2])
print('check4rotated: rotated trace {} to ZNE'.format(trace_ids)) print('check4rotated: rotated trace {} to ZNE'.format(trace_ids))
# replace old data with rotated data, change the channel code to ZNE # replace old data with rotated data, change the channel code to ZNE
z_index = dips.index(min( 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) dips)) # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90
wfstream[z_index].data = z # being vertical)
wfstream[z_index].stats.channel = wfstream[z_index].stats.channel[0:-1] + 'Z' 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] del trace_ids[z_index]
for trace_id in trace_ids: for trace_id in trace_ids:
coordinates = metadata.get_coordinates(trace_id, t_start) coordinates = metadata.get_coordinates(trace_id, t_start)
dip, az = coordinates['dip'], coordinates['azimuth'] dip, az = coordinates['dip'], coordinates['azimuth']
trace = wfstream.select(id=trace_id)[0] trace = wfs_out.select(id=trace_id)[0]
if az > 315 or az <= 45 or az > 135 and az <= 225: if az > 315 or az <= 45 or 135 < az <= 225:
trace.data = n trace.data = n
trace.stats.channel = trace.stats.channel[0:-1] + 'N' trace.stats.channel = trace.stats.channel[0:-1] + 'N'
elif az > 45 and az <= 135 or az > 225 and az <= 315: elif 45 < az <= 135 or 225 < az <= 315:
trace.data = e trace.data = e
trace.stats.channel = trace.stats.channel[0:-1] + 'E' trace.stats.channel = trace.stats.channel[0:-1] + 'E'
except (ValueError) as e: except ValueError as err:
print(e) print(f"{err=} Rotation failed. Stream remains unchanged.")
return wfstream return wfs_in
return wfstream return wfs_out
if metadata is None: if metadata is None:
if verbosity: if verbosity:
@ -1032,38 +1147,6 @@ def check4rotated(data, metadata=None, verbosity=1):
return data return data
def scaleWFData(data, factor=None, components='all'):
"""
produce scaled waveforms from given waveform data and a scaling factor,
waveform may be selected by their components name
:param data: waveform data to be scaled
:type data: `~obspy.core.stream.Stream` object
:param factor: scaling factor
:type factor: float
:param components: components labels for the traces in data to be scaled by
the scaling factor (optional, default: 'all')
:type components: tuple
:return: scaled waveform data
:rtype: `~obspy.core.stream.Stream` object
"""
if components != 'all':
for comp in components:
if factor is None:
max_val = np.max(np.abs(data.select(component=comp)[0].data))
data.select(component=comp)[0].data /= 2 * max_val
else:
data.select(component=comp)[0].data /= 2 * factor
else:
for tr in data:
if factor is None:
max_val = float(np.max(np.abs(tr.data)))
tr.data /= 2 * max_val
else:
tr.data /= 2 * factor
return data
def runProgram(cmd, parameter=None): def runProgram(cmd, parameter=None):
""" """
run an external program specified by cmd with parameters input returning the run an external program specified by cmd with parameters input returning the
@ -1135,7 +1218,6 @@ def identifyPhase(phase):
return False return False
@lru_cache
def identifyPhaseID(phase): def identifyPhaseID(phase):
""" """
Returns phase id (capital P or S) Returns phase id (capital P or S)

509
pylot/core/util/widgets.py
View File

@ -7,6 +7,7 @@ Created on Wed Mar 19 11:27:35 2014
import copy import copy
import datetime import datetime
import getpass import getpass
import glob
import multiprocessing import multiprocessing
import os import os
import subprocess import subprocess
@ -16,6 +17,7 @@ import traceback
import matplotlib import matplotlib
import numpy as np import numpy as np
from pylot.core.io.phases import getQualitiesfromxml
matplotlib.use('QT5Agg') matplotlib.use('QT5Agg')
@ -49,11 +51,11 @@ from pylot.core.pick.utils import getSNR, earllatepicker, getnoisewin, \
from pylot.core.pick.compare import Comparison from pylot.core.pick.compare import Comparison
from pylot.core.pick.autopick import fmpicker from pylot.core.pick.autopick import fmpicker
from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS
from pylot.core.util.utils import prepTimeAxis, full_range, demeanTrace, isSorted, findComboBoxIndex, clims, \ from pylot.core.util.utils import prep_time_axis, full_range, demeanTrace, isSorted, findComboBoxIndex, clims, \
pick_linestyle_plt, pick_color_plt, \ pick_linestyle_plt, pick_color_plt, \
check4rotated, check4doubled, check_for_gaps_and_merge, check_for_nan, identifyPhase, \ check4rotated, check4doubled, check_for_gaps_and_merge, check_for_nan, identifyPhase, \
loopIdentifyPhase, trim_station_components, transformFilteroptions2String, \ loopIdentifyPhase, trim_station_components, transformFilteroptions2String, \
identifyPhaseID, get_bool, get_None, pick_color, getAutoFilteroptions, SetChannelComponents, \ identifyPhaseID, get_bool, get_none, pick_color, getAutoFilteroptions, SetChannelComponents, \
station_id_remove_channel, get_pylot_eventfile_with_extension, get_possible_pylot_eventfile_extensions station_id_remove_channel, get_pylot_eventfile_with_extension, get_possible_pylot_eventfile_extensions
from autoPyLoT import autoPyLoT from autoPyLoT import autoPyLoT
from pylot.core.util.thread import Thread from pylot.core.util.thread import Thread
@ -138,18 +140,6 @@ class LogWidget(QtWidgets.QWidget):
self.stderr.append(60 * '#' + '\n\n') 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, def plot_pdf(_axes, x, y, annotation, bbox_props, xlabel=None, ylabel=None,
title=None): title=None):
# try method or data # try method or data
@ -793,7 +783,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
def connect_signals(self): def connect_signals(self):
self.qcombo_processed.activated.connect(self.parent().newWF) 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): def init_labels(self):
self.label_layout.addWidget(self.status_label) self.label_layout.addWidget(self.status_label)
@ -804,13 +794,13 @@ class WaveformWidgetPG(QtWidgets.QWidget):
# use widgets as placeholder, so that child widgets keep position when others are hidden # use widgets as placeholder, so that child widgets keep position when others are hidden
mid_layout = QHBoxLayout() mid_layout = QHBoxLayout()
right_layout = QHBoxLayout() right_layout = QHBoxLayout()
mid_layout.addWidget(self.syn_checkbox) mid_layout.addWidget(self.comp_checkbox)
right_layout.addWidget(self.qcombo_processed) right_layout.addWidget(self.qcombo_processed)
mid_widget.setLayout(mid_layout) mid_widget.setLayout(mid_layout)
right_widget.setLayout(right_layout) right_widget.setLayout(right_layout)
self.label_layout.addWidget(mid_widget) self.label_layout.addWidget(mid_widget)
self.label_layout.addWidget(right_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(0, 4)
self.label_layout.setStretch(1, 0) self.label_layout.setStretch(1, 0)
self.label_layout.setStretch(2, 0) self.label_layout.setStretch(2, 0)
@ -825,7 +815,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
label = QtWidgets.QLabel() label = QtWidgets.QLabel()
self.perm_labels.append(label) self.perm_labels.append(label)
self.qcombo_processed = QtWidgets.QComboBox() 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('processed', 'green')
self.addQCboxItem('raw', 'black') self.addQCboxItem('raw', 'black')
# self.perm_qcbox_right.setAlignment(2) # self.perm_qcbox_right.setAlignment(2)
@ -834,9 +824,11 @@ class WaveformWidgetPG(QtWidgets.QWidget):
def getPlotDict(self): def getPlotDict(self):
return self.plotdict return self.plotdict
def activateObspyDMToptions(self, activate): def activateObspyDMToptions(self, activate: bool) -> None:
self.syn_checkbox.setVisible(activate) self.qcombo_processed.setEnabled(activate)
self.qcombo_processed.setVisible(activate)
def activateCompareOptions(self, activate: bool) -> None:
self.comp_checkbox.setEnabled(activate)
def setPermText(self, number, text=None, color='black'): def setPermText(self, number, text=None, color='black'):
if not 0 <= number < len(self.perm_labels): if not 0 <= number < len(self.perm_labels):
@ -936,10 +928,10 @@ class WaveformWidgetPG(QtWidgets.QWidget):
msg = 'plotting %s channel of station %s' % (channel, station) msg = 'plotting %s channel of station %s' % (channel, station)
print(msg) print(msg)
stime = trace.stats.starttime - self.wfstart stime = trace.stats.starttime - self.wfstart
time_ax = prepTimeAxis(stime, trace) time_ax = prep_time_axis(stime, trace)
if st_syn: if st_syn:
stime_syn = trace_syn.stats.starttime - self.wfstart stime_syn = trace_syn.stats.starttime - self.wfstart
time_ax_syn = prepTimeAxis(stime_syn, trace_syn) time_ax_syn = prep_time_axis(stime_syn, trace_syn)
if method == 'fast': if method == 'fast':
trace.data, time_ax = self.minMax(trace, time_ax) trace.data, time_ax = self.minMax(trace, time_ax)
@ -959,7 +951,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
[time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else []) [time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else [])
trace.data = np.array( trace.data = np.array(
[datum * gain + n for index, datum in enumerate(trace.data) if not index % nth_sample]) [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 []) if not index % nth_sample] if st_syn else [])
plots.append((times, trace.data, plots.append((times, trace.data,
times_syn, trace_syn.data)) times_syn, trace_syn.data))
@ -1005,15 +997,6 @@ class WaveformWidgetPG(QtWidgets.QWidget):
time_ax = np.linspace(time_ax[0], time_ax[-1], num=len(data)) time_ax = np.linspace(time_ax[0], time_ax[-1], num=len(data))
return data, time_ax 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): def setXLims(self, lims):
vb = self.plotWidget.getPlotItem().getViewBox() vb = self.plotWidget.getPlotItem().getViewBox()
vb.setXRange(float(lims[0]), float(lims[1]), padding=0) vb.setXRange(float(lims[0]), float(lims[1]), padding=0)
@ -1148,12 +1131,12 @@ class PylotCanvas(FigureCanvas):
ax.set_xlim(self.cur_xlim) ax.set_xlim(self.cur_xlim)
ax.set_ylim(self.cur_ylim) ax.set_ylim(self.cur_ylim)
self.refreshPickDlgText() self.refreshPickDlgText()
ax.figure.canvas.draw() ax.figure.canvas.draw_idle()
def panRelease(self, gui_event): def panRelease(self, gui_event):
self.press = None self.press = None
self.press_rel = None self.press_rel = None
self.figure.canvas.draw() self.figure.canvas.draw_idle()
def panZoom(self, gui_event, threshold=2., factor=1.1): def panZoom(self, gui_event, threshold=2., factor=1.1):
if not gui_event.x and not gui_event.y: if not gui_event.x and not gui_event.y:
@ -1167,8 +1150,6 @@ class PylotCanvas(FigureCanvas):
break break
if not ax_check: return if not ax_check: return
# self.updateCurrentLimits() #maybe put this down to else:
# calculate delta (relative values in axis) # calculate delta (relative values in axis)
old_x, old_y = self.press_rel old_x, old_y = self.press_rel
xdiff = gui_event.x - old_x xdiff = gui_event.x - old_x
@ -1371,110 +1352,145 @@ class PylotCanvas(FigureCanvas):
plot_positions[channel] = plot_pos plot_positions[channel] = plot_pos
return plot_positions 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, noiselevel=None, scaleddata=False, mapping=True,
component='*', nth_sample=1, iniPick=None, verbosity=0, component='*', nth_sample=1, iniPick=None, verbosity=0,
plot_additional=False, additional_channel=None, scaleToChannel=None, plot_additional=False, additional_channel=None, scaleToChannel=None,
snr=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 = self.axes[0]
ax.cla() ax.cla()
return ax
self.clearPlotDict() def get_wf_range(self, wfdata):
wfstart, wfend = full_range(wfdata) return full_range(wfdata)
nmax = 0
def get_comp_class(self):
settings = QSettings() settings = QSettings()
compclass = SetChannelComponents.from_qsettings(settings) return SetChannelComponents.from_qsettings(settings)
if not component == '*': def get_plot_streams(self, wfdata, wfdata_compare, component, compclass):
alter_comp = compclass.getCompPosition(component) def get_wf_dict(data=Stream(), linecolor='k', offset=0., **plot_kwargs):
# alter_comp = str(alter_comp[0]) return dict(data=data, linecolor=linecolor, offset=offset, plot_kwargs=plot_kwargs)
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()
linecolor = (0., 0., 0., 1.) if not self.style else self.style['linecolor']['rgba_mpl'] 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): for n, seed_id in enumerate(nslc):
network, station, location, channel = seed_id.split('.') network, station, location, channel = seed_id.split('.')
st = st_select.select(id=seed_id) for wf_name, wf_dict in plot_streams.items():
trace = st[0].copy() st_select = wf_dict.get('data')
if not st_select:
continue
trace = st_select.select(id=seed_id)[0].copy()
if mapping: if mapping:
n = plot_positions[trace.stats.channel] n = plot_positions[trace.stats.channel]
if n > nmax: if n > nmax:
nmax = n nmax = n
if verbosity: if verbosity:
msg = 'plotting %s channel of station %s' % (channel, station) print(f'plotting {channel} channel of station {station}')
print(msg) time_ax = prep_time_axis(trace.stats.starttime - wfstart, trace)
stime = trace.stats.starttime - wfstart self.plot_trace(ax, trace, time_ax, wf_dict, n, scaleToChannel, noiselevel, scaleddata, nth_sample)
time_ax = prepTimeAxis(stime, trace) 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 time_ax is not None:
if scaleToChannel: 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) st_scale = wfdata.select(channel=scaleToChannel)
if st_scale: if st_scale:
tr = st_scale[0] tr = st_scale[0]
trace.detrend('constant') trace.detrend('constant')
trace.normalize(np.max(np.abs(tr.data)) * 2) 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] def plot_noise_level(self, ax, time_ax, noiselevel, channel, n, linecolor):
data = [datum + n for index, datum in enumerate(trace.data) if not index % nth_sample]
ax.axhline(n, color="0.5", lw=0.5)
ax.plot(times, data, color=linecolor, linewidth=0.7)
if noiselevel is not None:
for level in [-noiselevel[channel], noiselevel[channel]]: for level in [-noiselevel[channel], noiselevel[channel]]:
ax.plot([time_ax[0], time_ax[-1]], ax.plot([time_ax[0], time_ax[-1]], [n + level, n + level], color=linecolor, linestyle='dashed')
[n + level, n + level],
color=linecolor, def plot_additional_trace(self, ax, wfdata, additional_channel, scaleToChannel, scaleddata, nth_sample, wfstart):
linestyle='dashed')
self.setPlotDict(n, seed_id)
if plot_additional and additional_channel:
compare_stream = wfdata.select(channel=additional_channel) compare_stream = wfdata.select(channel=additional_channel)
if compare_stream: if compare_stream:
trace = compare_stream[0] trace = compare_stream[0]
if scaleToChannel: if scaleToChannel:
st_scale = wfdata.select(channel=scaleToChannel) self.scale_trace(trace, scaleToChannel)
if st_scale:
tr = st_scale[0]
trace.detrend('constant')
trace.normalize(np.max(np.abs(tr.data)) * 2)
scaleddata = True scaleddata = True
if not scaleddata: if not scaleddata:
trace.detrend('constant') trace.detrend('constant')
trace.normalize(np.max(np.abs(trace.data)) * 2) trace.normalize(np.max(np.abs(trace.data)) * 2)
time_ax = prepTimeAxis(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] times = [time for index, time in enumerate(time_ax) if not index % nth_sample]
p_data = compare_stream[0].data p_data = trace.data
# #normalize
# p_max = max(abs(p_data))
# p_data /= p_max
for index in range(3): for index in range(3):
ax.plot(times, p_data, color='red', alpha=0.5, linewidth=0.7) ax.plot(times, p_data, color='red', alpha=0.5, linewidth=0.7)
p_data += 1 p_data += 1
def finalize_plot(self, ax, wfstart, wfend, nmax, zoomx, zoomy, iniPick, title, snr):
if iniPick: if iniPick:
ax.vlines(iniPick, ax.get_ylim()[0], ax.get_ylim()[1], ax.vlines(iniPick, ax.get_ylim()[0], ax.get_ylim()[1], colors='m', linestyles='dashed', linewidth=2)
colors='m', linestyles='dashed', xlabel = f'seconds since {wfstart}'
linewidth=2)
xlabel = 'seconds since {0}'.format(wfstart)
ylabel = '' ylabel = ''
self.updateWidget(xlabel, ylabel, title) self.updateWidget(xlabel, ylabel, title)
self.setXLims(ax, [0, wfend - wfstart]) self.setXLims(ax, [0, wfend - wfstart])
@ -1484,15 +1500,14 @@ class PylotCanvas(FigureCanvas):
if zoomy is not None: if zoomy is not None:
self.setYLims(ax, zoomy) self.setYLims(ax, zoomy)
if snr is not None: if snr is not None:
if snr < 2: self.plot_snr_warning(ax, snr)
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.draw() 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 @staticmethod
def getXLims(ax): def getXLims(ax):
return ax.get_xlim() return ax.get_xlim()
@ -1574,6 +1589,8 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog):
self.events = events self.events = events
self.filepaths = [] self.filepaths = []
self.file_extensions = [] self.file_extensions = []
self.check_all_state = True
self.merge_strategy = None
self.default_text = default_text self.default_text = default_text
self.label = label self.label = label
self.setButtons() self.setButtons()
@ -1581,16 +1598,17 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog):
self.connectSignals() self.connectSignals()
self.showPaths() self.showPaths()
self.refreshSelectionBox() self.refreshSelectionBox()
self.refresh_timer = QTimer(self) # self.refresh_timer = QTimer(self)
self.refresh_timer.timeout.connect(self.showPaths) # self.refresh_timer.timeout.connect(self.showPaths)
self.refresh_timer.start(10000) # self.refresh_timer.start(10000)
self.resize(800, 450) self.resize(800, 450)
def setupUi(self): def setupUi(self):
ncol = 4
self.main_layout = QtWidgets.QVBoxLayout() self.main_layout = QtWidgets.QVBoxLayout()
self.header_layout = QtWidgets.QHBoxLayout() self.header_layout = QtWidgets.QHBoxLayout()
self.footer_layout = QtWidgets.QHBoxLayout()
# #
self.setLayout(self.main_layout) self.setLayout(self.main_layout)
@ -1604,11 +1622,24 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog):
self.searchButton = QtWidgets.QPushButton('Search') self.searchButton = QtWidgets.QPushButton('Search')
self.searchButton.setVisible(False) self.searchButton.setVisible(False)
# check/uncheck button for table
self.checkAllButton = QtWidgets.QPushButton('Check/Uncheck all')
# radiobutton for merge selection
self.mergeRadioButtonGroup = QtWidgets.QButtonGroup()
self.merge_button = QtWidgets.QRadioButton('Merge')
self.overwrite_button = QtWidgets.QRadioButton('Overwrite')
self.mergeRadioButtonGroup.addButton(self.merge_button)
self.mergeRadioButtonGroup.addButton(self.overwrite_button)
self.merge_button.setChecked(True)
self.merge_strategy = self.merge_button.text()
# table
self.tableWidget = QtWidgets.QTableWidget() self.tableWidget = QtWidgets.QTableWidget()
tableWidget = self.tableWidget tableWidget = self.tableWidget
tableWidget.setColumnCount(3) tableWidget.setColumnCount(ncol)
tableWidget.setRowCount(len(self.events)) tableWidget.setRowCount(len(self.events))
tableWidget.setHorizontalHeaderLabels(('Event ID', 'Filename', 'Last modified')) tableWidget.setHorizontalHeaderLabels(('', 'Event ID', 'Filename', 'Last modified'))
tableWidget.setEditTriggers(tableWidget.NoEditTriggers) tableWidget.setEditTriggers(tableWidget.NoEditTriggers)
tableWidget.setSortingEnabled(True) tableWidget.setSortingEnabled(True)
header = tableWidget.horizontalHeader() header = tableWidget.horizontalHeader()
@ -1621,9 +1652,17 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog):
self.header_layout.addWidget(self.comboBox) self.header_layout.addWidget(self.comboBox)
self.header_layout.addWidget(self.searchButton) self.header_layout.addWidget(self.searchButton)
self.footer_layout.addWidget(self.checkAllButton)
self.footer_layout.addWidget(self.statusText)
self.footer_layout.addWidget(self.merge_button)
self.footer_layout.addWidget(self.overwrite_button)
self.footer_layout.setStretch(0, 0)
self.footer_layout.setStretch(1, 1)
self.main_layout.addLayout(self.header_layout) self.main_layout.addLayout(self.header_layout)
self.main_layout.addWidget(self.tableWidget) self.main_layout.addWidget(self.tableWidget)
self.main_layout.addWidget(self.statusText) self.main_layout.addLayout(self.footer_layout)
self.main_layout.addWidget(self._buttonbox) self.main_layout.addWidget(self._buttonbox)
def showPaths(self): def showPaths(self):
@ -1632,23 +1671,23 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog):
self.tableWidget.clearContents() self.tableWidget.clearContents()
for index, event in enumerate(self.events): for index, event in enumerate(self.events):
filename = get_pylot_eventfile_with_extension(event, fext) filename = get_pylot_eventfile_with_extension(event, fext)
self.tableWidget.setItem(index, 0, QtWidgets.QTableWidgetItem(f'{event.pylot_id}')) pf_selected_item = QtWidgets.QTableWidgetItem()
check_state = QtCore.Qt.Checked if filename else QtCore.Qt.Unchecked
pf_selected_item.setCheckState(check_state)
self.tableWidget.setItem(index, 0, pf_selected_item)
self.tableWidget.setItem(index, 1, QtWidgets.QTableWidgetItem(f'{event.pylot_id}'))
if filename: if filename:
self.filepaths.append(filename) self.filepaths.append(filename)
ts = int(os.path.getmtime(filename)) ts = int(os.path.getmtime(filename))
# create QTableWidgetItems of filepath and last modification time # create QTableWidgetItems of filepath and last modification time
fname_item = QtWidgets.QTableWidgetItem(f'{os.path.split(filename)[-1]}') fname_item = QtWidgets.QTableWidgetItem(f'{os.path.split(filename)[-1]}')
fname_item.setData(3, filename)
ts_item = QtWidgets.QTableWidgetItem(f'{datetime.datetime.fromtimestamp(ts)}') ts_item = QtWidgets.QTableWidgetItem(f'{datetime.datetime.fromtimestamp(ts)}')
self.tableWidget.setItem(index, 1, fname_item) self.tableWidget.setItem(index, 2, fname_item)
self.tableWidget.setItem(index, 2, ts_item) self.tableWidget.setItem(index, 3, ts_item)
# TODO: Idea -> only refresh if table contents changed. Use selection to load only a subset of files self.update_status()
if len(self.filepaths) > 0:
status_text = f'Found {len(self.filepaths)} eventfiles. Do you want to load them?'
else:
status_text = 'Did not find any files for specified file mask.'
self.statusText.setText(status_text)
def refreshSelectionBox(self): def refreshSelectionBox(self):
fext = self.comboBox.currentText() fext = self.comboBox.currentText()
@ -1668,12 +1707,52 @@ class SearchFileByExtensionDialog(QtWidgets.QDialog):
self._buttonbox = QDialogButtonBox(QDialogButtonBox.Ok | self._buttonbox = QDialogButtonBox(QDialogButtonBox.Ok |
QDialogButtonBox.Cancel) QDialogButtonBox.Cancel)
def toggleCheckAll(self):
self.check_all_state = not self.check_all_state
self.checkAll(self.check_all_state)
def checkAll(self, state):
state = QtCore.Qt.Checked if state else QtCore.Qt.Unchecked
for row_ind in range(self.tableWidget.rowCount()):
item = self.tableWidget.item(row_ind, 0)
item.setCheckState(state)
def getChecked(self):
filepaths = []
for row_ind in range(self.tableWidget.rowCount()):
item_check = self.tableWidget.item(row_ind, 0)
if item_check.checkState() == QtCore.Qt.Checked:
item_fname = self.tableWidget.item(row_ind, 2)
if item_fname:
filepath = item_fname.data(3)
filepaths.append(filepath)
return filepaths
def update_status(self, row=None, col=None):
if col is not None and col != 0:
return
filepaths = self.getChecked()
if len(filepaths) > 0:
status_text = f"Found {len(filepaths)} eventfile{'s' if len(filepaths) > 1 else ''}. Do you want to load them?"
else:
status_text = 'Did not find/select any files for specified file mask.'
self.statusText.setText(status_text)
def update_merge_strategy(self):
for button in (self.merge_button, self.overwrite_button):
if button.isChecked():
self.merge_strategy = button.text()
def connectSignals(self): def connectSignals(self):
self._buttonbox.accepted.connect(self.accept) self._buttonbox.accepted.connect(self.accept)
self._buttonbox.rejected.connect(self.reject) self._buttonbox.rejected.connect(self.reject)
self.comboBox.editTextChanged.connect(self.showPaths) self.comboBox.editTextChanged.connect(self.showPaths)
self.searchButton.clicked.connect(self.showPaths) self.searchButton.clicked.connect(self.showPaths)
self.checkAllButton.clicked.connect(self.toggleCheckAll)
self.checkAllButton.clicked.connect(self.update_status)
self.tableWidget.cellClicked.connect(self.update_status)
self.merge_button.clicked.connect(self.update_merge_strategy)
self.overwrite_button.clicked.connect(self.update_merge_strategy)
class SingleTextLineDialog(QtWidgets.QDialog): class SingleTextLineDialog(QtWidgets.QDialog):
@ -1780,8 +1859,8 @@ class PhaseDefaults(QtWidgets.QDialog):
class PickDlg(QDialog): class PickDlg(QDialog):
update_picks = QtCore.Signal(dict) update_picks = QtCore.Signal(dict)
def __init__(self, parent=None, data=None, station=None, network=None, location=None, picks=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, autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None, show_comp_data=False,
event=None, filteroptions=None, model=None, wftype=None): event=None, filteroptions=None, model=None, wftype=None):
super(PickDlg, self).__init__(parent, Qt.Window) super(PickDlg, self).__init__(parent, Qt.Window)
self.orig_parent = parent self.orig_parent = parent
@ -1790,6 +1869,7 @@ class PickDlg(QDialog):
# initialize attributes # initialize attributes
self.parameter = parameter self.parameter = parameter
self._embedded = embedded self._embedded = embedded
self.showCompData = show_comp_data
self.station = station self.station = station
self.network = network self.network = network
self.location = location self.location = location
@ -1828,11 +1908,32 @@ class PickDlg(QDialog):
else: else:
self.filteroptions = FILTERDEFAULTS self.filteroptions = FILTERDEFAULTS
self.pick_block = False 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 ') self.nextStation = QtWidgets.QCheckBox('Continue with next station ')
# comparison channel # comparison channel
self.compareChannel = QtWidgets.QComboBox() self.referenceChannel = QtWidgets.QComboBox()
self.compareChannel.activated.connect(self.resetPlot) 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 # scale channel
self.scaleChannel = QtWidgets.QComboBox() self.scaleChannel = QtWidgets.QComboBox()
@ -1845,19 +1946,7 @@ class PickDlg(QDialog):
self.cur_xlim = None self.cur_xlim = None
self.cur_ylim = None self.cur_ylim = None
# set attribute holding data self.stime, self.etime = full_range(self.get_wf_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())
# initialize plotting widget # initialize plotting widget
self.multicompfig = PylotCanvas(parent=self, multicursor=True) self.multicompfig = PylotCanvas(parent=self, multicursor=True)
@ -1868,12 +1957,12 @@ class PickDlg(QDialog):
self.setupUi() self.setupUi()
# fill compare and scale channels # fill compare and scale channels
self.compareChannel.addItem('-', None) self.referenceChannel.addItem('-', None)
self.scaleChannel.addItem('individual', None) self.scaleChannel.addItem('individual', None)
for trace in self.getWFData(): for trace in self.get_wf_data():
channel = trace.stats.channel 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']: if not channel[-1] in ['Z', 'N', 'E', '1', '2', '3']:
print('Skipping unknown channel for scaling: {}'.format(channel)) print('Skipping unknown channel for scaling: {}'.format(channel))
continue continue
@ -1890,7 +1979,7 @@ class PickDlg(QDialog):
if self.wftype is not None: if self.wftype is not None:
title += ' | ({})'.format(self.wftype) 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) title=title)
self.multicompfig.setZoomBorders2content() self.multicompfig.setZoomBorders2content()
@ -2066,8 +2155,11 @@ class PickDlg(QDialog):
_dialtoolbar.addWidget(est_label) _dialtoolbar.addWidget(est_label)
_dialtoolbar.addWidget(self.plot_arrivals_button) _dialtoolbar.addWidget(self.plot_arrivals_button)
_dialtoolbar.addSeparator() _dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtWidgets.QLabel('Compare to channel: ')) _dialtoolbar.addWidget(QtWidgets.QLabel('Plot reference channel: '))
_dialtoolbar.addWidget(self.compareChannel) _dialtoolbar.addWidget(self.referenceChannel)
_dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtWidgets.QLabel('Compare: '))
_dialtoolbar.addWidget(self.compareCB)
_dialtoolbar.addSeparator() _dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtWidgets.QLabel('Scaling: ')) _dialtoolbar.addWidget(QtWidgets.QLabel('Scaling: '))
_dialtoolbar.addWidget(self.scaleChannel) _dialtoolbar.addWidget(self.scaleChannel)
@ -2152,10 +2244,12 @@ class PickDlg(QDialog):
station_id = trace.get_id() station_id = trace.get_id()
starttime = trace.stats.starttime starttime = trace.stats.starttime
station_coords = self.metadata.get_coordinates(station_id, starttime) station_coords = self.metadata.get_coordinates(station_id, starttime)
if not station_coords:
print('get_arrivals: No station coordinates found. Return!')
return
origins = self.pylot_event.origins origins = self.pylot_event.origins
if phases == ['None', 'None']: if phases == ['None', 'None']:
print("get_arrivals: Creation info (manual or auto) not available!") print("get_arrivals: Creation info (manual or auto) not available! Return!")
print("Return!")
return return
if origins: if origins:
source_origin = origins[0] source_origin = origins[0]
@ -2267,7 +2361,7 @@ class PickDlg(QDialog):
# create action and add to menu # create action and add to menu
# phase name transferred using lambda function # phase name transferred using lambda function
slot = lambda phase=phase, phaseID=phaseID: phaseSelect[phaseID](phase) slot = lambda ph=phase, phID=phaseID: phaseSelect[phID](ph)
picksAction = createAction(parent=self, text=phase, picksAction = createAction(parent=self, text=phase,
slot=slot, slot=slot,
shortcut=shortcut) shortcut=shortcut)
@ -2402,7 +2496,7 @@ class PickDlg(QDialog):
def activatePicking(self): def activatePicking(self):
self.leave_rename_phase() self.leave_rename_phase()
self.renamePhaseAction.setEnabled(False) self.renamePhaseAction.setEnabled(False)
self.compareChannel.setEnabled(False) self.referenceChannel.setEnabled(False)
self.scaleChannel.setEnabled(False) self.scaleChannel.setEnabled(False)
phase = self.currentPhase phase = self.currentPhase
phaseID = self.getPhaseID(phase) phaseID = self.getPhaseID(phase)
@ -2420,7 +2514,7 @@ class PickDlg(QDialog):
if self.autoFilterAction.isChecked(): if self.autoFilterAction.isChecked():
for filteraction in [self.filterActionP, self.filterActionS]: for filteraction in [self.filterActionP, self.filterActionS]:
filteraction.setChecked(False) filteraction.setChecked(False)
self.filterWFData() self.filter_wf_data()
self.draw() self.draw()
else: else:
self.draw() self.draw()
@ -2434,7 +2528,7 @@ class PickDlg(QDialog):
self.disconnectPressEvent() self.disconnectPressEvent()
self.multicompfig.connectEvents() self.multicompfig.connectEvents()
self.renamePhaseAction.setEnabled(True) self.renamePhaseAction.setEnabled(True)
self.compareChannel.setEnabled(True) self.referenceChannel.setEnabled(True)
self.scaleChannel.setEnabled(True) self.scaleChannel.setEnabled(True)
self.connect_pick_delete() self.connect_pick_delete()
self.draw() self.draw()
@ -2504,9 +2598,15 @@ class PickDlg(QDialog):
def getGlobalLimits(self, ax, axis): def getGlobalLimits(self, ax, axis):
return self.multicompfig.getGlobalLimits(ax, axis) return self.multicompfig.getGlobalLimits(ax, axis)
def getWFData(self): def get_wf_data(self):
return self.data return self.data
def get_wf_dataComp(self):
if self.showCompData:
return self.data_compare
else:
return Stream()
def selectWFData(self, channel): def selectWFData(self, channel):
component = channel[-1].upper() component = channel[-1].upper()
wfdata = Stream() wfdata = Stream()
@ -2516,17 +2616,17 @@ class PickDlg(QDialog):
return tr return tr
if component == 'E' or component == 'N': if component == 'E' or component == 'N':
for trace in self.getWFData(): for trace in self.get_wf_data():
trace = selectTrace(trace, 'NE') trace = selectTrace(trace, 'NE')
if trace: if trace:
wfdata.append(trace) wfdata.append(trace)
elif component == '1' or component == '2': elif component == '1' or component == '2':
for trace in self.getWFData(): for trace in self.get_wf_data():
trace = selectTrace(trace, '12') trace = selectTrace(trace, '12')
if trace: if trace:
wfdata.append(trace) wfdata.append(trace)
else: else:
wfdata = self.getWFData().select(component=component) wfdata = self.get_wf_data().select(component=component)
return wfdata return wfdata
def getPicks(self, picktype='manual'): def getPicks(self, picktype='manual'):
@ -2628,11 +2728,16 @@ class PickDlg(QDialog):
stime = self.getStartTime() stime = self.getStartTime()
# copy data for plotting # copy wfdata for plotting
data = self.getWFData().copy() wfdata = self.get_wf_data().copy()
data = self.getPickPhases(data, phase) wfdata_comp = self.get_wf_dataComp().copy()
data.normalize() wfdata = self.getPickPhases(wfdata, phase)
if not data: 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', QtWidgets.QMessageBox.warning(self, 'No channel to plot',
'No channel to plot for phase: {}. ' 'No channel to plot for phase: {}. '
'Make sure to select the correct channels for P and S ' 'Make sure to select the correct channels for P and S '
@ -2640,14 +2745,16 @@ class PickDlg(QDialog):
self.leave_picking_mode() self.leave_picking_mode()
return 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() filterphase = self.currentFilterPhase()
if filterphase: if filterphase:
filteroptions = self.getFilterOptions(filterphase).parseFilterOptions() filteroptions = self.getFilterOptions(filterphase).parseFilterOptions()
try: try:
data.detrend('linear') for wfd in [wfdata, wfdata_comp]:
data.filter(**filteroptions) if wfd:
# wfdata.filter(**filteroptions)# MP MP removed filtering of original data wfd.detrend('linear')
wfd.filter(**filteroptions)
# wfdata.filter(**filteroptions)# MP MP removed filtering of original wfdata
except ValueError as e: except ValueError as e:
self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Information, self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Information,
'Denied', 'Denied',
@ -2657,8 +2764,8 @@ class PickDlg(QDialog):
snr = [] snr = []
noiselevels = {} noiselevels = {}
# determine SNR and noiselevel # determine SNR and noiselevel
for trace in data.traces: for trace in wfdata.traces:
st = data.select(channel=trace.stats.channel) st = wfdata.select(channel=trace.stats.channel)
stime_diff = trace.stats.starttime - stime stime_diff = trace.stats.starttime - stime
result = getSNR(st, (noise_win, gap_win, signal_win), ini_pick - stime_diff) result = getSNR(st, (noise_win, gap_win, signal_win), ini_pick - stime_diff)
snr.append(result[0]) snr.append(result[0])
@ -2669,12 +2776,14 @@ class PickDlg(QDialog):
noiselevel = nfac noiselevel = nfac
noiselevels[trace.stats.channel] = noiselevel noiselevels[trace.stats.channel] = noiselevel
# prepare plotting of data # prepare plotting of wfdata
for trace in data: for wfd in [wfdata, wfdata_comp]:
t = prepTimeAxis(trace.stats.starttime - stime, trace) if wfd:
for trace in wfd:
t = prep_time_axis(trace.stats.starttime - stime, trace)
inoise = getnoisewin(t, ini_pick, noise_win, gap_win) inoise = getnoisewin(t, ini_pick, noise_win, gap_win)
trace = demeanTrace(trace, inoise) 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 channel = trace.stats.channel
noiselevel = noiselevels[channel] noiselevel = noiselevels[channel]
noiseScaleFactor = self.calcNoiseScaleFactor(noiselevel, zoomfactor=5.) noiseScaleFactor = self.calcNoiseScaleFactor(noiselevel, zoomfactor=5.)
@ -2685,7 +2794,7 @@ class PickDlg(QDialog):
x_res = getResolutionWindow(mean_snr, parameter.get('extent')) x_res = getResolutionWindow(mean_snr, parameter.get('extent'))
xlims = [ini_pick - x_res, ini_pick + x_res] 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 = self.getStation() + ' picking mode'
title += ' | SNR: {}'.format(mean_snr) title += ' | SNR: {}'.format(mean_snr)
@ -2693,9 +2802,10 @@ class PickDlg(QDialog):
filtops_str = transformFilteroptions2String(filteroptions) filtops_str = transformFilteroptions2String(filteroptions)
title += ' | Filteroptions: {}'.format(filtops_str) title += ' | Filteroptions: {}'.format(filtops_str)
plot_additional = bool(self.compareChannel.currentText()) plot_additional = bool(self.referenceChannel.currentText())
additional_channel = self.compareChannel.currentText() additional_channel = self.referenceChannel.currentText()
self.multicompfig.plotWFData(wfdata=data, self.multicompfig.plotWFData(wfdata=wfdata,
wfdata_compare=wfdata_comp,
title=title, title=title,
zoomx=xlims, zoomx=xlims,
zoomy=ylims, zoomy=ylims,
@ -2729,7 +2839,7 @@ class PickDlg(QDialog):
filteroptions = None filteroptions = None
# copy and filter data for earliest and latest possible picks # 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: if filteroptions:
try: try:
wfdata.detrend('linear') wfdata.detrend('linear')
@ -2776,7 +2886,7 @@ class PickDlg(QDialog):
minFMSNR = parameter.get('minFMSNR') minFMSNR = parameter.get('minFMSNR')
quality = get_quality_class(spe, parameter.get('timeerrorsP')) quality = get_quality_class(spe, parameter.get('timeerrorsP'))
if quality <= minFMweight and snr >= minFMSNR: 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) pick - stime_diff)
# save pick times for actual phase # save pick times for actual phase
@ -3068,7 +3178,7 @@ class PickDlg(QDialog):
def togglePickBlocker(self): def togglePickBlocker(self):
return not self.pick_block return not self.pick_block
def filterWFData(self, phase=None): def filter_wf_data(self, phase=None):
if not phase: if not phase:
phase = self.currentPhase phase = self.currentPhase
if self.getPhaseID(phase) == 'P': if self.getPhaseID(phase) == 'P':
@ -3086,7 +3196,8 @@ class PickDlg(QDialog):
self.cur_xlim = self.multicompfig.axes[0].get_xlim() self.cur_xlim = self.multicompfig.axes[0].get_xlim()
self.cur_ylim = self.multicompfig.axes[0].get_ylim() self.cur_ylim = self.multicompfig.axes[0].get_ylim()
# self.multicompfig.updateCurrentLimits() # self.multicompfig.updateCurrentLimits()
data = self.getWFData().copy() wfdata = self.get_wf_data().copy()
wfdata_comp = self.get_wf_dataComp().copy()
title = self.getStation() title = self.getStation()
if filter: if filter:
filtoptions = None filtoptions = None
@ -3094,19 +3205,22 @@ class PickDlg(QDialog):
filtoptions = self.getFilterOptions(self.getPhaseID(phase), gui_filter=True).parseFilterOptions() filtoptions = self.getFilterOptions(self.getPhaseID(phase), gui_filter=True).parseFilterOptions()
if filtoptions is not None: if filtoptions is not None:
data.detrend('linear') for wfd in [wfdata, wfdata_comp]:
data.taper(0.02, type='cosine') if wfd:
data.filter(**filtoptions) wfd.detrend('linear')
wfd.taper(0.02, type='cosine')
wfd.filter(**filtoptions)
filtops_str = transformFilteroptions2String(filtoptions) filtops_str = transformFilteroptions2String(filtoptions)
title += ' | Filteroptions: {}'.format(filtops_str) title += ' | Filteroptions: {}'.format(filtops_str)
if self.wftype is not None: if self.wftype is not None:
title += ' | ({})'.format(self.wftype) title += ' | ({})'.format(self.wftype)
plot_additional = bool(self.compareChannel.currentText()) plot_additional = bool(self.referenceChannel.currentText())
additional_channel = self.compareChannel.currentText() additional_channel = self.referenceChannel.currentText()
scale_channel = self.scaleChannel.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(), zoomx=self.getXLims(),
zoomy=self.getYLims(), zoomy=self.getYLims(),
plot_additional=plot_additional, plot_additional=plot_additional,
@ -3120,14 +3234,14 @@ class PickDlg(QDialog):
def filterP(self): def filterP(self):
self.filterActionS.setChecked(False) self.filterActionS.setChecked(False)
if self.filterActionP.isChecked(): if self.filterActionP.isChecked():
self.filterWFData('P') self.filter_wf_data('P')
else: else:
self.refreshPlot() self.refreshPlot()
def filterS(self): def filterS(self):
self.filterActionP.setChecked(False) self.filterActionP.setChecked(False)
if self.filterActionS.isChecked(): if self.filterActionS.isChecked():
self.filterWFData('S') self.filter_wf_data('S')
else: else:
self.refreshPlot() self.refreshPlot()
@ -3179,11 +3293,14 @@ class PickDlg(QDialog):
self.resetZoom() self.resetZoom()
self.refreshPlot() self.refreshPlot()
def switchCompData(self):
self.showCompData = self.compareCB.isChecked()
def refreshPlot(self): def refreshPlot(self):
if self.autoFilterAction.isChecked(): if self.autoFilterAction.isChecked():
self.filterActionP.setChecked(False) self.filterActionP.setChecked(False)
self.filterActionS.setChecked(False) self.filterActionS.setChecked(False)
# data = self.getWFData().copy() # data = self.get_wf_data().copy()
# title = self.getStation() # title = self.getStation()
filter = False filter = False
phase = None phase = None
@ -3683,8 +3800,8 @@ class TuneAutopicker(QWidget):
fnames = self.station_ids[self.get_current_station_id()] fnames = self.station_ids[self.get_current_station_id()]
if not fnames == self.fnames: if not fnames == self.fnames:
self.fnames = fnames self.fnames = fnames
self.data.setWFData(fnames) self.data.set_wf_data(fnames)
wfdat = self.data.getWFData() # all available streams wfdat = self.data.get_wf_data() # all available streams
# remove possible underscores in station names # remove possible underscores in station names
# wfdat = remove_underscores(wfdat) # wfdat = remove_underscores(wfdat)
# rotate misaligned stations to ZNE # rotate misaligned stations to ZNE
@ -3789,12 +3906,14 @@ class TuneAutopicker(QWidget):
network = None network = None
location = None location = None
wfdata = self.data.getWFData() wfdata = self.data.get_wf_data()
wfdata_comp = self.data.get_wf_dataComp()
metadata = self.parent().metadata metadata = self.parent().metadata
event = self.get_current_event() event = self.get_current_event()
filteroptions = self.parent().filteroptions filteroptions = self.parent().filteroptions
wftype = self.wftype if self.obspy_dmt else '' wftype = self.wftype if self.obspy_dmt else ''
self.pickDlg = PickDlg(self.parent(), data=wfdata.select(station=station).copy(), self.pickDlg = PickDlg(self.parent(), data=wfdata.select(station=station).copy(),
data_comp=wfdata_comp.select(station=station).copy(),
station=station, network=network, station=station, network=network,
location=location, parameter=self.parameter, location=location, parameter=self.parameter,
picks=self.get_current_event_picks(station), picks=self.get_current_event_picks(station),
@ -3843,7 +3962,7 @@ class TuneAutopicker(QWidget):
for plotitem in self._manual_pick_plots: for plotitem in self._manual_pick_plots:
self.clear_plotitem(plotitem) self.clear_plotitem(plotitem)
self._manual_pick_plots = [] self._manual_pick_plots = []
st = self.data.getWFData() st = self.data.get_wf_data()
tr = st.select(station=self.get_current_station())[0] tr = st.select(station=self.get_current_station())[0]
starttime = tr.stats.starttime starttime = tr.stats.starttime
# create two lists with figure names and subindices (for subplots) to get the correct axes # create two lists with figure names and subindices (for subplots) to get the correct axes
@ -3879,7 +3998,7 @@ class TuneAutopicker(QWidget):
self.plot_manual_pick_to_ax(ax=ax, picks=picks, phase='S', self.plot_manual_pick_to_ax(ax=ax, picks=picks, phase='S',
starttime=starttime, quality=qualitySpick) starttime=starttime, quality=qualitySpick)
for canvas in self.parent().canvas_dict.values(): for canvas in self.parent().canvas_dict.values():
canvas.draw() canvas.draw_idle()
def plot_manual_pick_to_ax(self, ax, picks, phase, starttime, quality): def plot_manual_pick_to_ax(self, ax, picks, phase, starttime, quality):
mpp = picks[phase]['mpp'] - starttime mpp = picks[phase]['mpp'] - starttime
@ -4794,8 +4913,8 @@ class InputsTab(PropTab):
self.tstopBox = QSpinBox() self.tstopBox = QSpinBox()
for spinbox in [self.tstartBox, self.tstopBox]: for spinbox in [self.tstartBox, self.tstopBox]:
spinbox.setRange(-99999, 99999) spinbox.setRange(-99999, 99999)
self.tstartBox.setValue(float(settings.value('tstart')) if get_None(settings.value('tstart')) else 0) self.tstartBox.setValue(float(settings.value('tstart')) if get_none(settings.value('tstart')) else 0)
self.tstopBox.setValue(float(settings.value('tstop')) if get_None(settings.value('tstop')) else 1e6) self.tstopBox.setValue(float(settings.value('tstop')) if get_none(settings.value('tstop')) else 1e6)
self.cuttimesLayout.addWidget(self.tstartBox, 10) self.cuttimesLayout.addWidget(self.tstartBox, 10)
self.cuttimesLayout.addWidget(QLabel('[s] and'), 0) self.cuttimesLayout.addWidget(QLabel('[s] and'), 0)
self.cuttimesLayout.addWidget(self.tstopBox, 10) self.cuttimesLayout.addWidget(self.tstopBox, 10)
@ -5791,7 +5910,7 @@ class ChooseWaveFormWindow(QWidget):
#self.currentSpectro = self.traces[ #self.currentSpectro = self.traces[
# self.chooseBoxTraces.currentText()[3:]][self.chooseBoxComponent.currentText()].spectrogram(show=False, title=t) # self.chooseBoxTraces.currentText()[3:]][self.chooseBoxComponent.currentText()].spectrogram(show=False, title=t)
#self.currentSpectro.show() #self.currentSpectro.show()
applyFFT() self.applyFFT()
def applyFFT(self, trace): def applyFFT(self, trace):
tra = self.traces[self.chooseBoxTraces.currentText()[3:]]['Z'] tra = self.traces[self.chooseBoxTraces.currentText()[3:]]['Z']

15
requirements.txt
View File

@ -1,12 +1,9 @@
# This file may be used to create an environment using: # This file may be used to create an environment using:
# $ conda create --name <env> --file <this file> # $ conda create --name <env> --file <this file>
# platform: win-64 # platform: win-64
cartopy=0.20.2 cartopy>=0.20.2
matplotlib-base=3.3.4 numpy<2
numpy=1.22.3 obspy>=1.3.0
obspy=1.3.0 pyqtgraph>=0.12.4
pyqtgraph=0.12.4 pyside2>=5.13.2
pyside2=5.13.2 scipy>=1.8.0
python=3.8.12
qt=5.12.9
scipy=1.8.0