52 changed files with 800 additions and 69633 deletions
-1
View File
@@ -2,4 +2,3 @@
*~ *~
.idea .idea
pylot/RELEASE-VERSION pylot/RELEASE-VERSION
/tests/test_autopicker/dmt_database_test/
-39
View File
@@ -1,39 +0,0 @@
Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <Darius_A@web.de>
Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <darius.arnold@rub.de>
Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <darius.arnold@ruhr-uni-bochum.de>
Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <mail@dariusarnold.de>
Dennis Wlecklik <dennisw@minos02.geophysik.ruhr-uni-bochum.de>
Jeldrik Gaal <jeldrikgaal@gmail.com>
Kaan Coekerim <kaan.coekerim@ruhr-uni-bochum.de>
Kaan Coekerim <kaan.coekerim@ruhr-uni-bochum.de> <kaan.coekerim@rub.de>
Ludger Kueperkoch <kueperkoch@igem-energie.de> <kueperkoch@bestec-for-nature.com>
Ludger Kueperkoch <kueperkoch@igem-energie.de> <ludger@quake2.(none)>
Ludger Kueperkoch <kueperkoch@igem-energie.de> <ludger@sauron.bestec-for-nature>
Marc S. Boxberg <marc.boxberg@rub.de>
Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel.paffrath@rub.de>
Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@minos01.geophysik.ruhr-uni-bochum.de>
Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@minos02.geophysik.ruhr-uni-bochum.de>
Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@minos25.geophysik.ruhr-uni-bochum.de>
Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@email.com>
Sally Zimmermann <sally.zimmermann@ruhr-uni-bochum.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastianw@minos01.geophysik.ruhr-uni-bochum.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastianw@minos02.geophysik.ruhr-uni-bochum.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastianw@minos22.geophysik.ruhr-uni-bochum.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastian.wehling-benatelli@scisys.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastian.wehling@rub.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastian.wehling@rub.de>
Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <DarkBeQst@users.noreply.github.com>
Thomas Moeller <thomas.moeller@rub.de>
Ann-Christin Koch <ann-christin.koch@ruhr-uni-bochum.de> <Ann-Christin.Koch@ruhr-uni-bochum.de>
Sebastian Priebe <sebastian.priebe@rub.de>
Regular → Executable
+62 -67
View File
@@ -83,7 +83,7 @@ from pylot.core.util.event import Event
from pylot.core.io.location import create_creation_info, create_event from pylot.core.io.location import create_creation_info, create_event
from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \
PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \ PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \
ComparisonWidget, TuneAutopicker, PylotParameterWidget, AutoPickDlg, CanvasWidget, AutoPickWidget, \ ComparisonWidget, TuneAutopicker, PylotParaBox, AutoPickDlg, CanvasWidget, AutoPickWidget, \
CompareEventsWidget, ProgressBarWidget, AddMetadataWidget, SingleTextLineDialog, LogWidget, PickQualitiesFromXml, \ CompareEventsWidget, ProgressBarWidget, AddMetadataWidget, SingleTextLineDialog, LogWidget, PickQualitiesFromXml, \
SpectrogramTab, SearchFileByExtensionDialog SpectrogramTab, SearchFileByExtensionDialog
from pylot.core.util.array_map import Array_map from pylot.core.util.array_map import Array_map
@@ -136,7 +136,7 @@ class MainWindow(QMainWindow):
self.project.parameter = self._inputs self.project.parameter = self._inputs
self.tap = None self.tap = None
self.apw = None self.apw = None
self.parameterWidget = None self.paraBox = None
self.array_map = None self.array_map = None
self._metadata = Metadata(verbosity=0) self._metadata = Metadata(verbosity=0)
self._eventChanged = [False, False] self._eventChanged = [False, False]
@@ -188,6 +188,7 @@ class MainWindow(QMainWindow):
self.table_headers = ['', 'Event', 'Time', 'Lat', 'Lon', 'Depth', 'Ml', 'Mw', '[N] MP', '[N] AP', 'Tuning Set', self.table_headers = ['', 'Event', 'Time', 'Lat', 'Lon', 'Depth', 'Ml', 'Mw', '[N] MP', '[N] AP', 'Tuning Set',
'Test Set', 'Notes'] 'Test Set', 'Notes']
# TODO: refactor rootpath to datapath
while True: while True:
try: try:
if settings.value("user/FullName", None) is None: if settings.value("user/FullName", None) is None:
@@ -685,9 +686,10 @@ class MainWindow(QMainWindow):
# add scroll area used in case number of traces gets too high # add scroll area used in case number of traces gets too high
self.wf_scroll_area = QtWidgets.QScrollArea(self) self.wf_scroll_area = QtWidgets.QScrollArea(self)
self.wf_scroll_area.setVisible(False) self.wf_scroll_area.setVisible(False)
self.no_data_label = QLabel('No Data. If data were already loaded, try to select the event again in the eventbox.') self.no_data_label = QLabel('No Data')
self.no_data_label.setStyleSheet('color: red') self.no_data_label.setStyleSheet('color: red')
self.no_data_label.setAlignment(Qt.AlignCenter) self.no_data_label.setAlignment(Qt.AlignCenter)
# create central matplotlib figure canvas widget # create central matplotlib figure canvas widget
self.init_wfWidget() self.init_wfWidget()
@@ -716,14 +718,14 @@ class MainWindow(QMainWindow):
self.tabs.addTab(wf_tab, 'Waveform Plot') self.tabs.addTab(wf_tab, 'Waveform Plot')
self.tabs.addTab(array_tab, 'Array Map') self.tabs.addTab(array_tab, 'Array Map')
self.tabs.addTab(events_tab, 'Eventlist') self.tabs.addTab(events_tab, 'Eventlist')
#self.tabs.addTab(spectro_tab, 'Spectro') self.tabs.addTab(spectro_tab, 'Spectro')
self.wf_layout.addWidget(self.no_data_label) self.wf_layout.addWidget(self.no_data_label)
self.wf_layout.addWidget(self.wf_scroll_area) self.wf_layout.addWidget(self.wf_scroll_area)
self.wf_scroll_area.setWidgetResizable(True) self.wf_scroll_area.setWidgetResizable(True)
self.init_array_tab() self.init_array_tab()
self.init_event_table() self.init_event_table()
#self.init_spectro_tab() self.init_spectro_tab()
self.tabs.setCurrentIndex(0) self.tabs.setCurrentIndex(0)
self.eventLabel = QLabel() self.eventLabel = QLabel()
@@ -1030,24 +1032,8 @@ class MainWindow(QMainWindow):
fname = self.filename_from_action(action) fname = self.filename_from_action(action)
if not fname: if not fname:
return return
if not event: if not event:
event = self.get_current_event() event = self.get_current_event()
if event.picks:
qmb = QMessageBox(self, icon=QMessageBox.Question,
text='Do you want to overwrite the data?',)
overwrite_button = qmb.addButton('Overwrite', QMessageBox.YesRole)
merge_button = qmb.addButton('Merge', QMessageBox.NoRole)
qmb.exec_()
if qmb.clickedButton() == overwrite_button:
merge_strategy = 'Overwrite'
elif qmb.clickedButton() == merge_button:
merge_strategy = 'Merge'
else:
return
data = Data(self, event) data = Data(self, event)
try: try:
data_new = Data(self, evtdata=str(fname)) data_new = Data(self, evtdata=str(fname))
@@ -1210,7 +1196,7 @@ class MainWindow(QMainWindow):
with open(eventlist_file, 'r') as infile: with open(eventlist_file, 'r') as infile:
eventlist_subset = [os.path.join(basepath, filename.split('\n')[0]) for filename in eventlist_subset = [os.path.join(basepath, filename.split('\n')[0]) for filename in
infile.readlines()] infile.readlines()]
msg = 'Found file "eventlist.txt" in datapath. WILL ONLY USE SELECTED EVENTS out of {} events ' \ msg = 'Found file "eventlist.txt" in database path. WILL ONLY USE SELECTED EVENTS out of {} events ' \
'contained in this subset' 'contained in this subset'
print(msg.format(len(eventlist_subset))) print(msg.format(len(eventlist_subset)))
eventlist = [eventname for eventname in eventlist if eventname in eventlist_subset] eventlist = [eventname for eventname in eventlist if eventname in eventlist_subset]
@@ -1235,34 +1221,49 @@ class MainWindow(QMainWindow):
# get path from first event in list and split them # get path from first event in list and split them
path = eventlist[0] path = eventlist[0]
try: try:
datapath = os.path.split(path)[0] system_name = platform.system()
dirs = { if system_name in ["Linux", "Darwin"]:
'datapath': datapath, dirs = {
} 'database': path.split('/')[-2],
'datapath': os.path.split(path)[0], # path.split('/')[-3],
'rootpath': '/' + os.path.join(*path.split('/')[:-3])
}
elif system_name == "Windows":
rootpath = path.split('/')[:-3]
rootpath[0] += '/'
dirs = {
# TODO: Arrange path to meet Win standards
'database': path.split('/')[-2],
'datapath': path.split('/')[-3],
'rootpath': os.path.join(*rootpath)
}
except Exception as e: except Exception as e:
dirs = { dirs = {
'database': '',
'datapath': '', 'datapath': '',
'rootpath': ''
} }
print('Warning: Could not automatically init folder structure. ({})'.format(e)) print('Warning: Could not automatically init folder structure. ({})'.format(e))
settings = QSettings() settings = QSettings()
settings.setValue("data/dataRoot", dirs['datapath']) settings.setValue("data/dataRoot", dirs['datapath']) # d irs['rootpath'])
settings.sync() settings.sync()
if not self.project.eventlist: if not self.project.eventlist:
# init parameter object # init parameter object
self.setParameter(show=False) self.setParameter(show=False)
# hide all parameter (show all needed parameter later) # hide all parameter (show all needed parameter later)
self.parameterWidget.hide_parameter() self.paraBox.hide_parameter()
for directory in dirs.keys(): for directory in dirs.keys():
# set parameter # set parameter
box = self.parameterWidget.boxes[directory] box = self.paraBox.boxes[directory]
self.parameterWidget.setValue(box, dirs[directory]) self.paraBox.setValue(box, dirs[directory])
# show needed parameter in box # show needed parameter in box
self.parameterWidget.show_parameter(directory) self.paraBox.show_parameter(directory)
dirs_box = self.parameterWidget.get_groupbox_dialog('Directories') dirs_box = self.paraBox.get_groupbox_dialog('Directories')
if not dirs_box.exec_(): if not dirs_box.exec_():
return return
self.project.rootpath = dirs['rootpath']
self.project.datapath = dirs['datapath'] self.project.datapath = dirs['datapath']
else: else:
if hasattr(self.project, 'datapath'): if hasattr(self.project, 'datapath'):
@@ -1271,6 +1272,7 @@ class MainWindow(QMainWindow):
'Datapath missmatch to current project!') 'Datapath missmatch to current project!')
return return
else: else:
self.project.rootpath = dirs['rootpath']
self.project.datapath = dirs['datapath'] self.project.datapath = dirs['datapath']
self.project.add_eventlist(eventlist) self.project.add_eventlist(eventlist)
@@ -1358,10 +1360,11 @@ class MainWindow(QMainWindow):
return True return True
def modify_project_path(self, new_rootpath): def modify_project_path(self, new_rootpath):
self.project.datapath = new_rootpath # TODO: change root to datapath
self.project.rootpath = new_rootpath
for event in self.project.eventlist: for event in self.project.eventlist:
event.datapath = new_rootpath event.rootpath = new_rootpath
event.path = os.path.join(event.datapath, event.pylot_id) event.path = os.path.join(event.rootpath, event.datapath, event.database, event.pylot_id)
event.path = event.path.replace('\\', '/') event.path = event.path.replace('\\', '/')
event.path = event.path.replace('//', '/') event.path = event.path.replace('//', '/')
@@ -1537,8 +1540,8 @@ class MainWindow(QMainWindow):
return True return True
return not bool(os.listdir(wf_path)) return not bool(os.listdir(wf_path))
def filename_from_action(self, action=None): def filename_from_action(self, action):
if not action or action.data() is None: if action.data() is None:
filt = "Supported file formats" \ filt = "Supported file formats" \
" (*.mat *.qml *.xml *.kor *.evt)" " (*.mat *.qml *.xml *.kor *.evt)"
caption = "Open an event file" caption = "Open an event file"
@@ -1550,12 +1553,7 @@ class MainWindow(QMainWindow):
fname = str(action.data().toString()) fname = str(action.data().toString())
return fname return fname
def getEventFileName(self, type='manual'): def saveData(self, event=None, directory=None, outformats=['.xml', '.cnv', '.obs', '_focmec.in', '.pha']):
if self.get_fnames(type) is None:
self.set_fname(self.get_data().getEventFileName(), type)
return self.get_fnames(type)
def saveData(self, event=None, directory=None, outformats=None):
''' '''
Save event data to directory with specified output formats. Save event data to directory with specified output formats.
:param event: PyLoT Event, if not set current event will be used :param event: PyLoT Event, if not set current event will be used
@@ -1563,8 +1561,6 @@ class MainWindow(QMainWindow):
:param outformats: str/list of output formats :param outformats: str/list of output formats
:return: :return:
''' '''
if outformats is None:
outformats = ['.xml', '.cnv', '.obs', '_focmec.in', '.pha']
if not event: if not event:
event = self.get_current_event() event = self.get_current_event()
if not type(outformats) == list: if not type(outformats) == list:
@@ -1700,7 +1696,7 @@ class MainWindow(QMainWindow):
# WIP JG # WIP JG
def eventlistXml(self): def eventlistXml(self):
path = self._inputs['datapath'] path = self._inputs['rootpath'] + '/' + self._inputs['datapath'] + '/' + self._inputs['database']
outpath = self.project.location[:self.project.location.rfind('/')] outpath = self.project.location[:self.project.location.rfind('/')]
geteventlistfromxml(path, outpath) geteventlistfromxml(path, outpath)
return return
@@ -1976,6 +1972,7 @@ class MainWindow(QMainWindow):
self.dataPlot.activateObspyDMToptions(self.obspy_dmt) self.dataPlot.activateObspyDMToptions(self.obspy_dmt)
if self.obspy_dmt: if self.obspy_dmt:
self.prepareObspyDMT_data(eventpath) self.prepareObspyDMT_data(eventpath)
self.dataPlot.activateCompareOptions(True)
def loadWaveformData(self): def loadWaveformData(self):
''' '''
@@ -2152,11 +2149,10 @@ class MainWindow(QMainWindow):
self.wf_scroll_area.setVisible(len(plots) > 0) self.wf_scroll_area.setVisible(len(plots) > 0)
self.no_data_label.setVisible(not len(plots) > 0) self.no_data_label.setVisible(not len(plots) > 0)
for times, data, times_syn, data_syn in plots: for times, data, times_syn, data_syn in plots:
self.dataPlot.plotWidget.getPlotItem().plot(np.array(times), np.array(data), self.dataPlot.plotWidget.getPlotItem().plot(times, data,
pen=self.dataPlot.pen_linecolor, pen=self.dataPlot.pen_linecolor)
skipFiniteCheck=True)
if len(data_syn) > 0: if len(data_syn) > 0:
self.dataPlot.plotWidget.getPlotItem().plot(np.array(times_syn), np.array(data_syn), self.dataPlot.plotWidget.getPlotItem().plot(times_syn, data_syn,
pen=self.dataPlot.pen_linecolor_syn) pen=self.dataPlot.pen_linecolor_syn)
self.dataPlot.reinitMoveProxy() self.dataPlot.reinitMoveProxy()
self.highlight_stations() self.highlight_stations()
@@ -2196,8 +2192,7 @@ class MainWindow(QMainWindow):
if event.pylot_autopicks: if event.pylot_autopicks:
self.drawPicks(picktype='auto') self.drawPicks(picktype='auto')
if event.pylot_picks or event.pylot_autopicks: if event.pylot_picks or event.pylot_autopicks:
if not self._inputs.get('extent') == 'global': self.locateEventAction.setEnabled(True)
self.locateEventAction.setEnabled(True)
self.qualities_action.setEnabled(True) self.qualities_action.setEnabled(True)
self.eventlist_xml_action.setEnabled(True) self.eventlist_xml_action.setEnabled(True)
@@ -2427,7 +2422,7 @@ class MainWindow(QMainWindow):
filterS = filteroptions['S'] filterS = filteroptions['S']
minP, maxP = filterP.getFreq() minP, maxP = filterP.getFreq()
minS, maxS = filterS.getFreq() minS, maxS = filterS.getFreq()
self.parameterWidget.params_to_gui() self.paraBox.params_to_gui()
def getFilterOptions(self): def getFilterOptions(self):
return self.filteroptions return self.filteroptions
@@ -2632,6 +2627,7 @@ class MainWindow(QMainWindow):
picks=self.getPicksOnStation(station, 'manual'), picks=self.getPicksOnStation(station, 'manual'),
autopicks=self.getPicksOnStation(station, 'auto'), autopicks=self.getPicksOnStation(station, 'auto'),
metadata=self.metadata, event=event, metadata=self.metadata, event=event,
model=self.inputs.get('taup_model'),
filteroptions=self.filteroptions, wftype=wftype, filteroptions=self.filteroptions, wftype=wftype,
show_comp_data=self.dataPlot.comp_checkbox.isChecked()) show_comp_data=self.dataPlot.comp_checkbox.isChecked())
if self.filterActionP.isChecked(): if self.filterActionP.isChecked():
@@ -3096,7 +3092,7 @@ class MainWindow(QMainWindow):
if self.pg: if self.pg:
if spe: if spe:
if not self.plot_method == 'fast' and picks['epp'] and picks['lpp']: if picks['epp'] and picks['lpp']:
pen = make_pen(picktype, phaseID, 'epp', quality) pen = make_pen(picktype, phaseID, 'epp', quality)
self.drawnPicks[picktype][station].append(pw.plot([epp, epp], ylims, self.drawnPicks[picktype][station].append(pw.plot([epp, epp], ylims,
alpha=.25, pen=pen, name='EPP')) alpha=.25, pen=pen, name='EPP'))
@@ -3175,8 +3171,8 @@ class MainWindow(QMainWindow):
ttt = parameter['ttpatter'] ttt = parameter['ttpatter']
outfile = parameter['outpatter'] outfile = parameter['outpatter']
eventname = self.get_current_event_name() eventname = self.get_current_event_name()
obsdir = os.path.join(self._inputs['datapath'], eventname) obsdir = os.path.join(self._inputs['rootpath'], self._inputs['datapath'], self._inputs['database'], eventname)
self.saveData(event=self.get_current_event(), directory=obsdir, outformats=['.obs']) self.saveData(event=self.get_current_event(), directory=obsdir, outformats='.obs')
filename = 'PyLoT_' + eventname filename = 'PyLoT_' + eventname
locpath = os.path.join(locroot, 'loc', filename) locpath = os.path.join(locroot, 'loc', filename)
phasefile = os.path.join(obsdir, filename + '.obs') phasefile = os.path.join(obsdir, filename + '.obs')
@@ -3733,7 +3729,6 @@ class MainWindow(QMainWindow):
if self.project.parameter: if self.project.parameter:
# do this step to update default parameter on older PyLoT projects # do this step to update default parameter on older PyLoT projects
self.project.parameter.reinit_default_parameters() self.project.parameter.reinit_default_parameters()
PylotParameter.check_deprecated_parameters(self.project.parameter)
self._inputs = self.project.parameter self._inputs = self.project.parameter
self.updateFilteroptions() self.updateFilteroptions()
@@ -3851,13 +3846,13 @@ class MainWindow(QMainWindow):
def setParameter(self, checked=0, show=True): def setParameter(self, checked=0, show=True):
if checked: pass # dummy argument to receive trigger signal (checked) if called by QAction if checked: pass # dummy argument to receive trigger signal (checked) if called by QAction
if not self.parameterWidget: if not self.paraBox:
self.parameterWidget = PylotParameterWidget(self._inputs, parent=self, windowflag=Qt.Window) self.paraBox = PylotParaBox(self._inputs, parent=self, windowflag=Qt.Window)
self.parameterWidget.accepted.connect(self._setDirty) self.paraBox.accepted.connect(self._setDirty)
self.parameterWidget.accepted.connect(self.filterOptionsFromParameter) self.paraBox.accepted.connect(self.filterOptionsFromParameter)
if show: if show:
self.parameterWidget.params_to_gui() self.paraBox.params_to_gui()
self.parameterWidget.show() self.paraBox.show()
def deleteAllAutopicks(self): def deleteAllAutopicks(self):
qmb = QMessageBox(self, icon=QMessageBox.Question, qmb = QMessageBox(self, icon=QMessageBox.Question,
@@ -3904,18 +3899,16 @@ class Project(object):
Pickable class containing information of a PyLoT project, like event lists and file locations. Pickable class containing information of a PyLoT project, like event lists and file locations.
''' '''
# TODO: remove rootpath
def __init__(self): def __init__(self):
self.eventlist = [] self.eventlist = []
self.location = None self.location = None
self.rootpath = None
self.datapath = None self.datapath = None
self.dirty = False self.dirty = False
self.parameter = None self.parameter = None
self._table = None self._table = None
@property
def rootpath(self):
return self.datapath
def add_eventlist(self, eventlist): def add_eventlist(self, eventlist):
''' '''
Add events from an eventlist containing paths to event directories. Add events from an eventlist containing paths to event directories.
@@ -3925,6 +3918,8 @@ class Project(object):
return return
for item in eventlist: for item in eventlist:
event = Event(item) event = Event(item)
event.rootpath = self.parameter['rootpath']
event.database = self.parameter['database']
event.datapath = self.parameter['datapath'] event.datapath = self.parameter['datapath']
if not event.path in self.getPaths(): if not event.path in self.getPaths():
self.eventlist.append(event) self.eventlist.append(event)
+30 -14
View File
@@ -11,7 +11,7 @@ PILOT has originally been developed in Mathworks' MatLab. In order to distribute
problems, it has been decided to redevelop the software package in Python. The great work of the ObsPy group allows easy problems, it has been decided to redevelop the software package in Python. The great work of the ObsPy group allows easy
handling of a bunch of seismic data and PyLoT will benefit a lot compared to the former MatLab version. handling of a bunch of seismic data and PyLoT will benefit a lot compared to the former MatLab version.
The development of PyLoT is part of the joint research project MAGS2, AlpArray and AdriaArray. The development of PyLoT is part of the joint research project MAGS2 and AlpArray.
## Installation ## Installation
@@ -27,41 +27,58 @@ Afterwards run (from the PyLoT main directory where the files *requirements.txt*
conda env create -f pylot.yml conda env create -f pylot.yml
or or
conda create -c conda-forge --name pylot_311 python=3.11 --file requirements.txt conda create --name pylot_38 --file requirements.txt
to create a new Anaconda environment called *pylot_311*. to create a new Anaconda environment called "pylot_38".
Afterwards activate the environment by typing Afterwards activate the environment by typing
conda activate pylot_311 conda activate pylot_38
#### Prerequisites: #### Prerequisites:
In order to run PyLoT you need to install: In order to run PyLoT you need to install:
- Python 3 - Python 3
- cartopy
- joblib
- obspy - obspy
- pyaml
- pyqtgraph
- pyside2 - pyside2
- pyqtgraph
- cartopy
(the following are already dependencies of the above packages): (the following are already dependencies of the above packages):
- scipy - scipy
- numpy - numpy
- matplotlib - matplotlib <= 3.3.x
#### Some handwork: #### Some handwork:
Some extra information on error estimates (just needed for reading old PILOT data) and the Richter magnitude scaling PyLoT needs a properties folder on your system to work. It should be situated in your home directory
(on Windows usually C:/Users/*username*):
mkdir ~/.pylot
In the next step you have to copy some files to this directory:
*for local distance seismicity*
cp path-to-pylot/inputs/pylot_local.in ~/.pylot/pylot.in
*for regional distance seismicity*
cp path-to-pylot/inputs/pylot_regional.in ~/.pylot/pylot.in
*for global distance seismicity*
cp path-to-pylot/inputs/pylot_global.in ~/.pylot/pylot.in
and some extra information on error estimates (just needed for reading old PILOT data) and the Richter magnitude scaling
relation relation
cp path-to-pylot/inputs/PILOT_TimeErrors.in path-to-pylot/inputs/richter_scaling.data ~/.pylot/ cp path-to-pylot/inputs/PILOT_TimeErrors.in path-to-pylot/inputs/richter_scaling.data ~/.pylot/
You may need to do some modifications to these files. Especially folder names should be reviewed. You may need to do some modifications to these files. Especially folder names should be reviewed.
PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10/11. PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10.
## Release notes ## Release notes
@@ -70,7 +87,6 @@ PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10/11.
- event organisation in project files and waveform visualisation - event organisation in project files and waveform visualisation
- consistent manual phase picking through predefined SNR dependant zoom level - consistent manual phase picking through predefined SNR dependant zoom level
- consistent automatic phase picking routines using Higher Order Statistics, AIC and Autoregression - consistent automatic phase picking routines using Higher Order Statistics, AIC and Autoregression
- pick correlation correction for teleseismic waveforms
- interactive tuning of auto-pick parameters - interactive tuning of auto-pick parameters
- uniform uncertainty estimation from waveform's properties for automatic and manual picks - uniform uncertainty estimation from waveform's properties for automatic and manual picks
- pdf representation and comparison of picks taking the uncertainty intrinsically into account - pdf representation and comparison of picks taking the uncertainty intrinsically into account
@@ -79,7 +95,7 @@ PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10/11.
#### Known issues: #### Known issues:
Current release is still in development progress and has several issues. We are currently lacking manpower, but hope to assess many of the issues in the near future. We hope to solve these with the next release.
## Staff ## Staff
@@ -92,4 +108,4 @@ Others: A. Bruestle, T. Meier, W. Friederich
[ObsPy]: http://github.com/obspy/obspy/wiki [ObsPy]: http://github.com/obspy/obspy/wiki
September 2024 April 2022
+13 -12
View File
@@ -136,9 +136,11 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
if parameter.hasParam('datastructure'): if parameter.hasParam('datastructure'):
# getting information on data structure # getting information on data structure
datastructure = DATASTRUCTURE[parameter.get('datastructure')]() datastructure = DATASTRUCTURE[parameter.get('datastructure')]()
dsfields = {'dpath': parameter.get('datapath'),} dsfields = {'root': parameter.get('rootpath'),
'dpath': parameter.get('datapath'),
'dbase': parameter.get('database')}
exf = ['dpath'] exf = ['root', 'dpath', 'dbase']
if parameter['eventID'] != '*' and fnames == 'None': if parameter['eventID'] != '*' and fnames == 'None':
dsfields['eventID'] = parameter['eventID'] dsfields['eventID'] = parameter['eventID']
@@ -184,15 +186,15 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
if not input_dict: if not input_dict:
# started in production mode # started in production mode
datapath = datastructure.expandDataPath() datapath = datastructure.expandDataPath()
if fnames in [None, 'None'] and parameter['eventID'] == '*': if fnames == 'None' and parameter['eventID'] == '*':
# multiple event processing # multiple event processing
# read each event in database # read each event in database
events = [event for event in glob.glob(os.path.join(datapath, '*')) if events = [event for event in glob.glob(os.path.join(datapath, '*')) if
(os.path.isdir(event) and not event.endswith('EVENTS-INFO'))] (os.path.isdir(event) and not event.endswith('EVENTS-INFO'))]
elif fnames in [None, 'None'] and parameter['eventID'] != '*' and not type(parameter['eventID']) == list: elif fnames == 'None' and parameter['eventID'] != '*' and not type(parameter['eventID']) == list:
# single event processing # single event processing
events = glob.glob(os.path.join(datapath, parameter['eventID'])) events = glob.glob(os.path.join(datapath, parameter['eventID']))
elif fnames in [None, 'None'] and type(parameter['eventID']) == list: elif fnames == 'None' and type(parameter['eventID']) == list:
# multiple event processing # multiple event processing
events = [] events = []
for eventID in parameter['eventID']: for eventID in parameter['eventID']:
@@ -204,10 +206,12 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
locflag = 2 locflag = 2
else: else:
# started in tune or interactive mode # started in tune or interactive mode
datapath = parameter['datapath'] datapath = os.path.join(parameter['rootpath'],
parameter['datapath'])
events = [] events = []
for eventID in eventid: for eventID in eventid:
events.append(os.path.join(datapath, events.append(os.path.join(datapath,
parameter['database'],
eventID)) eventID))
if not events: if not events:
@@ -234,15 +238,12 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
data.get_evt_data().path = eventpath data.get_evt_data().path = eventpath
print('Reading event data from filename {}...'.format(filename)) print('Reading event data from filename {}...'.format(filename))
except Exception as e: except Exception as e:
if type(e) == FileNotFoundError: print('Could not read event from file {}: {}'.format(filename, e))
print('Creating new event file.')
else:
print('Could not read event from file {}: {}'.format(filename, e))
data = Data() data = Data()
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 in [None, 'None']: if fnames == 'None':
data.setWFData(glob.glob(os.path.join(event_datapath, '*'))) data.setWFData(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
+2
View File
@@ -203,6 +203,8 @@ The meaning of the header entries is:
PyLoT GUI starts with an empty project. To add events, use the add event data button. Select one or multiple folders PyLoT GUI starts with an empty project. To add events, use the add event data button. Select one or multiple folders
containing events. containing events.
[//]: <> (TODO: explain _Directories: Root path, Data path, Database path_)
### Saving projects ### Saving projects
Save the current project from the menu with File->Save project or File->Save project as. PyLoT uses ``.plp`` files to Save the current project from the menu with File->Save project or File->Save project as. PyLoT uses ``.plp`` files to
+3 -2
View File
@@ -4,8 +4,10 @@
%Parameters are optimized for %extent data sets! %Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings# #main settings#
#rootpath# %project path
#datapath# %data path #datapath# %data path
#eventID# %event ID for single event processing (* for all events found in datapath) #database# %name of data base
#eventID# %event ID for single event processing (* for all events found in database)
#invdir# %full path to inventory or dataless-seed file #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output True #apverbose# %choose 'True' or 'False' for terminal output
@@ -41,7 +43,6 @@ global #extent# %extent of a
1150.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking 1150.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
iasp91 #taup_model# %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6 iasp91 #taup_model# %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
P,Pdiff #taup_phases# %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
0.05 0.5 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 0.05 0.5 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
0.001 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] 0.001 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
0.05 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] 0.05 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
+4 -3
View File
@@ -4,8 +4,10 @@
%Parameters are optimized for %extent data sets! %Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings# #main settings#
/DATA/Insheim/EVENT_DATA/LOCAL/2018.02_Insheim #datapath# %data path /DATA/Insheim #rootpath# %project path
e0006.038.18 #eventID# %event ID for single event processing (* for all events found in datapath) EVENT_DATA/LOCAL #datapath# %data path
2018.02_Insheim #database# %name of data base
e0006.038.18 #eventID# %event ID for single event processing (* for all events found in database)
/DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file /DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output True #apverbose# %choose 'True' or 'False' for terminal output
@@ -41,7 +43,6 @@ local #extent# %extent of a
10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking 10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
False #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF False #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
iasp91 #taup_model# %define TauPy model for traveltime estimation iasp91 #taup_model# %define TauPy model for traveltime estimation
P #taup_phases# %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
2.0 20.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 2.0 20.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
2.0 30.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] 2.0 30.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
2.0 10.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] 2.0 10.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
+3 -2
View File
@@ -4,8 +4,10 @@
%Parameters are optimized for %extent data sets! %Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings# #main settings#
#rootpath# %project path
#datapath# %data path #datapath# %data path
#eventID# %event ID for single event processing (* for all events found in datapath) #database# %name of data base
#eventID# %event ID for single event processing (* for all events found in database)
#invdir# %full path to inventory or dataless-seed file #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output True #apverbose# %choose 'True' or 'False' for terminal output
@@ -41,7 +43,6 @@ local #extent# %extent of a
10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking 10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
iasp91 #taup_model# %define TauPy model for traveltime estimation iasp91 #taup_model# %define TauPy model for traveltime estimation
P #taup_phases# %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
2.0 10.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 2.0 10.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
2.0 12.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] 2.0 12.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
2.0 8.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] 2.0 8.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
+10 -8
View File
@@ -1,12 +1,14 @@
name: pylot_311 name: pylot_38
channels: channels:
- conda-forge - conda-forge
- defaults - defaults
dependencies: dependencies:
- cartopy=0.23.0=py311hcf9f919_1 - cartopy=0.20.2
- joblib=1.4.2=pyhd8ed1ab_0 - matplotlib-base=3.3.4
- obspy=1.4.1=py311he736701_3 - numpy=1.22.3
- pyaml=24.7.0=pyhd8ed1ab_0 - obspy=1.3.0
- pyqtgraph=0.13.7=pyhd8ed1ab_0 - pyqtgraph=0.12.4
- pyside2=5.15.8=py311h3d699ce_4 - pyside2>=5.13.2
- pytest=8.3.2=pyhd8ed1ab_0 - python=3.8.12
- qt>=5.12.9
- scipy=1.8.0
+1 -1
View File
@@ -9,7 +9,7 @@ PyLoT - the Python picking and Localization Tool
This python library contains a graphical user interfaces for picking This python library contains a graphical user interfaces for picking
seismic phases. This software needs ObsPy (http://github.com/obspy/obspy/wiki) seismic phases. This software needs ObsPy (http://github.com/obspy/obspy/wiki)
and the Qt libraries to be installed first. and the Qt4 libraries to be installed first.
PILOT has been developed in Mathworks' MatLab. In order to distribute PILOT has been developed in Mathworks' MatLab. In order to distribute
PILOT without facing portability problems, it has been decided to re- PILOT without facing portability problems, it has been decided to re-
+170 -242
View File
@@ -4,13 +4,16 @@
import copy import copy
import logging import logging
import os import os
import fnmatch
from dataclasses import dataclass, field
from typing import List
from PySide2.QtWidgets import QMessageBox from PySide2.QtWidgets import QMessageBox
from obspy import read_events from obspy import read, read_events, Stream, Catalog, UTCDateTime
from obspy.core import read, Stream, UTCDateTime
from obspy.core.event import Event as ObsPyEvent from obspy.core.event import Event as ObsPyEvent
from obspy.io.sac import SacIOError from obspy.io.sac import SacIOError
import pylot.core.loc.focmec as focmec import pylot.core.loc.focmec as focmec
import pylot.core.loc.hypodd as hypodd import pylot.core.loc.hypodd as hypodd
import pylot.core.loc.velest as velest import pylot.core.loc.velest as velest
@@ -139,29 +142,6 @@ class Data(object):
def setNew(self): def setNew(self):
self._new = True self._new = True
def getCutTimes(self):
"""
Returns earliest start and latest end of all waveform data
:return: minimum start time and maximum end time as a tuple
:rtype: (UTCDateTime, UTCDateTime)
"""
if self.cuttimes is None:
self.updateCutTimes()
return self.cuttimes
def updateCutTimes(self):
"""
Update cuttimes to contain earliest start and latest end time
of all waveform data
:rtype: None
"""
self.cuttimes = full_range(self.getWFData())
def getEventFileName(self):
ID = self.getID()
# handle forbidden filenames especially on windows systems
return fnConstructor(str(ID))
def checkEvent(self, event, fcheck, forceOverwrite=False): def checkEvent(self, event, fcheck, forceOverwrite=False):
""" """
Check information in supplied event and own event and replace with own Check information in supplied event and own event and replace with own
@@ -252,184 +232,6 @@ class Data(object):
if picktype in str(pick.method_id.id): if picktype in str(pick.method_id.id):
picks.append(pick) 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:
parameter = PylotParameter()
logging.warning('Using default input parameter')
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:
parameter = PylotParameter()
logging.warning('Using default input parameter')
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): def getID(self):
""" """
Get unique resource id Get unique resource id
@@ -474,21 +276,6 @@ class Data(object):
fnames = check_fname_exists(fnames) fnames = check_fname_exists(fnames)
fnames_alt = check_fname_exists(fnames_alt) fnames_alt = check_fname_exists(fnames_alt)
# 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: if fnames is not None:
self.appendWFData(fnames) self.appendWFData(fnames)
if fnames_alt is not None: if fnames_alt is not None:
@@ -496,9 +283,6 @@ class Data(object):
else: else:
return False return False
# various pre-processing steps:
# remove possible underscores in station names
# self.wfdata = remove_underscores(self.wfdata)
# check for gaps and merge # check for gaps and merge
self.wfdata, _ = check_for_gaps_and_merge(self.wfdata) self.wfdata, _ = check_for_gaps_and_merge(self.wfdata)
# check for nans # check for nans
@@ -620,11 +404,6 @@ class Data(object):
picks = picks_from_picksdict(picks) picks = picks_from_picksdict(picks)
break break
self.get_evt_data().picks = picks 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): def applyEvent(event):
""" """
@@ -656,6 +435,171 @@ class Data(object):
applydata[typ](data) applydata[typ](data)
self._new = False self._new = False
@dataclass
class SeismicEventData:
event_id: str = ""
catalog: Catalog = field(default_factory=Catalog)
def find_event_files(self, directory: str, extensions: List[str]) -> List[str]:
"""
Browse the directory to find event files with specified extensions.
Parameters:
directory (str): The directory path to search for event files.
extensions (List[str]): List of file extensions to search for.
Returns:
List[str]: List of file paths that match the given extensions.
Example:
>>> sed = SeismicEventData()
>>> sed.find_event_files('test_directory', ['.xml', '.quakeml']) # doctest: +SKIP
['test_directory/event1.xml', 'test_directory/event2.quakeml']
"""
matches = []
for root, _, files in os.walk(directory):
for ext in extensions:
for filename in fnmatch.filter(files, f'*{ext}'):
matches.append(os.path.join(root, filename))
return matches
def read_event_from_directory(self, directory: str, extensions: List[str], format: str) -> None:
"""
Read a seismic event from the first found file in the directory with specified format.
Parameters:
directory (str): The directory path to search for event files.
extensions (List[str]): List of file extensions to search for.
format (str): The format to read the event file.
Example:
>>> sed = SeismicEventData()
>>> sed.read_event_from_directory('test_directory', ['.xml', '.quakeml'], 'QUAKEML') # doctest: +SKIP
"""
event_files = self.find_event_files(directory, extensions)
if event_files:
self.read_event(event_files[0], format)
else:
raise FileNotFoundError(f"No event files found in directory {directory} with extensions {extensions}.")
def read_event(self, file_path: str, format: str) -> None:
"""
Read a seismic event from a file with specified format.
Parameters:
file_path (str): The path to the event file.
format (str): The format to read the event file.
Example:
>>> sed = SeismicEventData()
>>> sed.read_event('test_directory/event1.xml', 'QUAKEML') # doctest: +SKIP
"""
if os.path.exists(file_path):
self.catalog = read_events(file_path, format=format)
self.event_id = self.catalog[0].resource_id.id.split('/')[-1] if self.catalog else ""
else:
raise FileNotFoundError(f"File {file_path} does not exist.")
def write_event(self, file_path: str, format: str) -> None:
"""
Write the seismic event to a file with specified format.
Parameters:
file_path (str): The path to the output file.
format (str): The format to write the event file.
Example:
>>> sed = SeismicEventData(event_id='12345')
>>> sed.write_event('output_directory/event1.xml', 'QUAKEML') # doctest: +SKIP
"""
self.catalog.write(file_path, format=format)
@dataclass
class WaveformData:
stream: Stream = field(default_factory=Stream)
def find_waveform_files(self, directory: str, extensions: List[str]) -> List[str]:
"""
Browse the directory to find waveform files with specified extensions.
Parameters:
directory (str): The directory path to search for waveform files.
extensions (List[str]): List of file extensions to search for.
Returns:
List[str]: List of file paths that match the given extensions.
Example:
>>> wd = WaveformData()
>>> wd.find_waveform_files('test_directory', ['.mseed']) # doctest: +SKIP
['test_directory/waveform1.mseed']
"""
matches = []
for root, _, files in os.walk(directory):
for ext in extensions:
for filename in fnmatch.filter(files, f'*{ext}'):
matches.append(os.path.join(root, filename))
return matches
def read_waveform_from_directory(self, directory: str, extensions: List[str], format: str) -> None:
"""
Read waveform data from the first found file in the directory with specified format.
Parameters:
directory (str): The directory path to search for waveform files.
extensions (List[str]): List of file extensions to search for.
format (str): The format to read the waveform file.
Example:
>>> wd = WaveformData()
>>> wd.read_waveform_from_directory('test_directory', ['.mseed'], 'MSEED') # doctest: +SKIP
"""
waveform_files = self.find_waveform_files(directory, extensions)
if waveform_files:
self.read_waveform(waveform_files[0], format)
else:
raise FileNotFoundError(f"No waveform files found in directory {directory} with extensions {extensions}.")
def read_waveform(self, file_path: str, format: str) -> None:
"""
Read waveform data from a file with specified format.
Parameters:
file_path (str): The path to the waveform file.
format (str): The format to read the waveform file.
Example:
>>> wd = WaveformData()
>>> wd.read_waveform('test_directory/waveform1.mseed', 'MSEED') # doctest: +SKIP
"""
if os.path.exists(file_path):
self.stream = read(file_path, format=format)
else:
raise FileNotFoundError(f"File {file_path} does not exist.")
def write_waveform(self, file_path: str, format: str) -> None:
"""
Write the waveform data to a file with specified format.
Parameters:
file_path (str): The path to the output file.
format (str): The format to write the waveform file.
Example:
>>> wd = WaveformData()
>>> wd.write_waveform('output_directory/waveform1.mseed', 'MSEED') # doctest: +SKIP
"""
self.stream.write(file_path, format=format)
# Example usage:
# seismic_event = SeismicEventData()
# seismic_event.read_event_from_directory("path_to_directory", extensions=[".xml", ".quakeml"], format="QUAKEML")
# seismic_event.write_event("output_event_file.xml", format="QUAKEML")
# waveform_data = WaveformData()
# waveform_data.read_waveform_from_directory("path_to_directory", extensions=[".mseed"], format="MSEED")
# waveform_data.write_waveform("output_waveform_file.mseed", format="MSEED")
class GenericDataStructure(object): class GenericDataStructure(object):
""" """
@@ -839,22 +783,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:
+17 -5
View File
@@ -6,14 +6,24 @@ import numpy as np
Default parameters used for picking Default parameters used for picking
""" """
defaults = {'datapath': {'type': str, defaults = {'rootpath': {'type': str,
'tooltip': 'path to eventfolders', 'tooltip': 'project path',
'value': '',
'namestring': 'Root path'},
'datapath': {'type': str,
'tooltip': 'data path',
'value': '', 'value': '',
'namestring': 'Data path'}, 'namestring': 'Data path'},
'database': {'type': str,
'tooltip': 'name of data base',
'value': '',
'namestring': 'Database path'},
'eventID': {'type': str, 'eventID': {'type': str,
'tooltip': 'event ID for single event processing (* for all events found in datapath)', 'tooltip': 'event ID for single event processing (* for all events found in database)',
'value': '*', 'value': '',
'namestring': 'Event ID'}, 'namestring': 'Event ID'},
'extent': {'type': str, 'extent': {'type': str,
@@ -501,7 +511,7 @@ defaults = {'datapath': {'type': str,
'taup_model': {'type': str, 'taup_model': {'type': str,
'tooltip': 'Define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6', 'tooltip': 'Define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6',
'value': 'iasp91', 'value': None,
'namestring': 'TauPy model'}, 'namestring': 'TauPy model'},
'taup_phases': {'type': str, 'taup_phases': {'type': str,
@@ -512,7 +522,9 @@ defaults = {'datapath': {'type': str,
settings_main = { settings_main = {
'dirs': [ 'dirs': [
'rootpath',
'datapath', 'datapath',
'database',
'eventID', 'eventID',
'invdir', 'invdir',
'datastructure', 'datastructure',
+3 -27
View File
@@ -1,7 +1,5 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import logging
import os
from pylot.core.io import default_parameters from pylot.core.io import default_parameters
from pylot.core.util.errors import ParameterError from pylot.core.util.errors import ParameterError
@@ -90,10 +88,10 @@ class PylotParameter(object):
return bool(self.__parameter) return bool(self.__parameter)
def __getitem__(self, key): def __getitem__(self, key):
if key in self.__parameter: try:
return self.__parameter[key] return self.__parameter[key]
else: except:
logging.warning(f'{key} not found in PylotParameter') return None
def __setitem__(self, key, value): def __setitem__(self, key, value):
try: try:
@@ -420,28 +418,6 @@ class PylotParameter(object):
line = value + name + ttip line = value + name + ttip
fid.write(line) fid.write(line)
@staticmethod
def check_deprecated_parameters(parameters):
if parameters.hasParam('database') and parameters.hasParam('rootpath'):
parameters['datapath'] = os.path.join(parameters['rootpath'], parameters['datapath'],
parameters['database'])
logging.warning(
f'Parameters database and rootpath are deprecated. '
f'Tried to merge them to now path: {parameters["datapath"]}.'
)
remove_keys = []
for key in parameters:
if not key in default_parameters.defaults.keys():
remove_keys.append(key)
logging.warning(f'Removing deprecated parameter: {key}')
for key in remove_keys:
del parameters[key]
parameters._settings_main = default_parameters.settings_main
parameters._settings_special_pick = default_parameters.settings_special_pick
class FilterOptions(object): class FilterOptions(object):
''' '''
+206 -663
View File
@@ -21,25 +21,6 @@ from pylot.core.util.utils import get_owner, full_range, four_digits, transformF
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
@@ -193,39 +174,12 @@ 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, parameter=None, nwst_id: bool = False):
""" """
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
:param evt: Event object contain all available information :param evt: Event object contain all available information
:type evt: `~obspy.core.event.Event` :type evt: `~obspy.core.event.Event`
:param nwst_id: determines if network and station id are used or only station id
:type nwst_id: bool
:return: pick dictionary (auto and manual) :return: pick dictionary (auto and manual)
""" """
picksdict = { picksdict = {
@@ -235,10 +189,7 @@ def picksdict_from_picks(evt, parameter=None, nwst_id: bool = False):
for pick in evt.picks: for pick in evt.picks:
errors = None errors = None
phase = {} phase = {}
if not nwst_id: station = pick.waveform_id.station_code
station_or_nwst = pick.waveform_id.station_code
else:
station_or_nwst = f'{pick.waveform_id.station_code}.{pick.waveform_id.network_code}'
if pick.waveform_id.channel_code is None: if pick.waveform_id.channel_code is None:
channel = '' channel = ''
else: else:
@@ -259,7 +210,7 @@ def picksdict_from_picks(evt, parameter=None, nwst_id: bool = False):
if pick_method == 'None': if pick_method == 'None':
pick_method = 'manual' pick_method = 'manual'
try: try:
onsets = picksdict[pick_method][station_or_nwst] onsets = picksdict[pick_method][station]
except KeyError as e: except KeyError as e:
# print(e) # print(e)
onsets = {} onsets = {}
@@ -306,7 +257,7 @@ def picksdict_from_picks(evt, parameter=None, nwst_id: bool = False):
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()
picksdict[pick_method][station_or_nwst] = onsets.copy() picksdict[pick_method][station] = onsets.copy()
return picksdict return picksdict
@@ -378,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
:type parameter: object
:param parameter: all input information :param eventinfo: Event information needed for specific formats like VELEST, FOCMEC, and HASH
:type parameter: object :type eventinfo: obspy.core.event.Event
:param eventinfo: optional, needed for VELEST-cnv file
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('datapath'), 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:
fid.write(
' {}\n'.format(parameter.get('eventID')))
for key, value in arrivals.items():
if value['P'].get('weight', 0) < 4:
stat = key[:4]
Ponset = value['P']['mpp']
Sonset = value.get('S', {}).get('mpp')
pweight = value['P'].get('weight', 0)
sweight = value.get('S', {}).get('weight', 0)
fm = value['P'].get('fm', '-')
Ao = value.get('S', {}).get('Ao', '')
year = Ponset.year - 2000 if Ponset.year >= 2000 else Ponset.year - 1900
ss_ms = Ponset.second + Ponset.microsecond / 1000000.0
if Sonset:
Sss_ms = Sonset.second + Sonset.microsecond / 1000000.0
fid.write('{}P{}{}{} {}{}{}{}{} {:5.2f} {}{}S {} {}\n'.format(
stat, 'I' if pweight < 2 else 'E', fm, pweight, year, Ponset.month, Ponset.day,
Ponset.hour, Ponset.minute, ss_ms, Sss_ms, 'I' if sweight < 2 else 'E', sweight, Ao))
else:
fid.write('{}P{}{}{} {}{}{}{}{} {:5.2f} {}\n'.format(
stat, 'I' if pweight < 2 else 'E', fm, pweight, year, Ponset.month, Ponset.day,
Ponset.hour, Ponset.minute, ss_ms, Ao))
elif fformat == 'HYPOSAT': def write_hyposat():
print("Writing phases to %s for HYPOSAT" % filename) with open(filename, 'w') as fid:
fid = open("%s" % filename, 'w') fid.write('{}, event {} \n'.format(parameter.get('database'), parameter.get('eventID')))
# write header for key, value in arrivals.items():
fid.write('%s, event %s \n' % (parameter.get('datapath'), parameter.get('eventID'))) for phase in ['P', 'S']:
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere if phase in value and value[phase].get('weight', 0) < 4:
for key in arrivals: onset = value[phase]['mpp']
# P onsets ss_ms = onset.second + onset.microsecond / 1000000.0
if 'P' in arrivals[key] and arrivals[key]['P']['mpp'] is not None: std = value[phase].get('spe', parameter.get('timeerrorsP')[value[phase].get('weight', 0)])
if arrivals[key]['P']['weight'] < 4: fid.write(
Ponset = arrivals[key]['P']['mpp'] '{:<5} {}1 {:4} {:02} {:02} {:02} {:02} {:05.02f} {:5.3f} -999. 0.00 -999. 0.00\n'.format(
pyear = Ponset.year key, phase, onset.year, onset.month, onset.day, onset.hour, onset.minute, ss_ms, std))
pmonth = Ponset.month
pday = Ponset.day
phh = Ponset.hour
pmm = Ponset.minute
pss = Ponset.second
pms = Ponset.microsecond
Pss = pss + pms / 1000000.0
# use symmetrized picking error as std
# (read the HYPOSAT manual)
pstd = arrivals[key]['P']['spe']
if pstd is None:
errorsP = parameter.get('timeerrorsP')
if arrivals[key]['P']['weight'] == 0:
pstd = errorsP[0]
elif arrivals[key]['P']['weight'] == 1:
pstd = errorsP[1]
elif arrivals[key]['P']['weight'] == 2:
pstd = errorsP[2]
elif arrivals[key]['P']['weight'] == 3:
psrd = errorsP[3]
else:
pstd = errorsP[4]
fid.write('%-5s P1 %4.0f %02d %02d %02d %02d %05.02f %5.3f -999. 0.00 -999. 0.00\n'
% (key, pyear, pmonth, pday, phh, pmm, Pss, pstd))
# S onsets
if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None:
if arrivals[key]['S']['weight'] < 4:
Sonset = arrivals[key]['S']['mpp']
syear = Sonset.year
smonth = Sonset.month
sday = Sonset.day
shh = Sonset.hour
smm = Sonset.minute
sss = Sonset.second
sms = Sonset.microsecond
Sss = sss + sms / 1000000.0
sstd = arrivals[key]['S']['spe']
if pstd is None:
errorsS = parameter.get('timeerrorsS')
if arrivals[key]['S']['weight'] == 0:
pstd = errorsS[0]
elif arrivals[key]['S']['weight'] == 1:
pstd = errorsS[1]
elif arrivals[key]['S']['weight'] == 2:
pstd = errorsS[2]
elif arrivals[key]['S']['weight'] == 3:
psrd = errorsS[3]
else:
pstd = errorsP[4]
fid.write('%-5s S1 %4.0f %02d %02d %02d %02d %05.02f %5.3f -999. 0.00 -999. 0.00\n'
% (key, syear, smonth, sday, shh, smm, Sss, sstd))
fid.close()
elif fformat == 'VELEST': def write_velest():
print("Writing phases to %s for VELEST" % filename) if not eventinfo:
fid = open("%s" % filename, 'w')
# get informations needed in cnv-file
# check, whether latitude is N or S and longitude is E or W
try:
eventsource = eventinfo.origins[0]
except:
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: stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second,
print("Error 1111111!") origin.latitude, origin.longitude, origin.depth / 1000, eventinfo.magnitudes[0].mag,
hddID = "00000" origin.quality.standard_error, "00000"))
# write header for key, value in arrivals.items():
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' % ( for phase in ['P', 'S']:
stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second, if phase in value and value[phase].get('weight', 0) < 4:
eventsource['latitude'], eventsource['longitude'], eventsource['depth'] / 1000, onset = value[phase]['mpp']
eventinfo.magnitudes[0]['mag'], eventsource['quality']['standard_error'], hddID)) rt = (onset - stime).total_seconds()
# check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT) fid.write('{} {:6.3f} 1 {}\n'.format(key, rt, phase))
if isinstance(arrivals, dict) == False:
# convert pick object (PyLoT) into dictionary
evt = ope.Event(resource_id=eventinfo['resource_id'])
evt.picks = arrivals
arrivals = picksdict_from_picks(evt)
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
if 'P' in usedarrivals[key]:
# P onsets
if usedarrivals[key]['P']['weight'] < 4:
Ponset = usedarrivals[key]['P']['mpp']
Prt = Ponset - stime # onset time relative to source time
fid.write('%s %6.3f 1 P\n' % (key, Prt))
if 'S' in usedarrivals[key]:
# S onsets
if usedarrivals[key]['S']['weight'] < 4:
Sonset = usedarrivals[key]['S']['mpp']
Srt = Sonset - stime # onset time relative to source time
fid.write('%-5s %6.3f 1 S\n' % (key, Srt))
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"
print(f"Writing phases to {filename1} for HASH-driver 1")
print(f"Writing phases to {filename2} for HASH-driver 2")
# Open files for writing
with open(filename1, 'w') as fid1, open(filename2, 'w') as fid2:
# Get event information needed for HASH-input file
try:
eventsource = eventinfo.origins[0]
except IndexError:
print("No source origin calculated yet, thus no cnv-file creation possible!")
return
event = parameter.get('eventID')
hashID = event.split('.')[0][1:5]
latdeg = eventsource['latitude']
latmin = (eventsource['latitude'] * 60) / 10000
londeg = eventsource['longitude']
lonmin = (eventsource['longitude'] * 60) / 10000
erh = (eventsource.origin_uncertainty['min_horizontal_uncertainty'] +
eventsource.origin_uncertainty['max_horizontal_uncertainty']) / 2000
erz = eventsource.depth_errors['uncertainty']
stime = eventsource['time']
syear = stime.year % 100 # Calculate two-digit year
picks = eventinfo.picks
# Write header line including event information for HASH-driver 1
fid1.write(f"{syear:02d}{stime.month:02d}{stime.day:02d}{stime.hour:02d}{stime.minute:02d}"
f"{stime.second:05.2f}{latdeg:2d}N{latmin:05.2f}{londeg:3d}E{lonmin:05.2f}"
f"{eventsource['depth']:6.2f}{eventinfo.magnitudes[0]['mag']:4.2f}{erh:5.2f}{erz:5.2f}{hashID}\n")
# Write header line including event information for HASH-driver 2
fid2.write(f"{syear:02d}{stime.month:02d}{stime.day:02d}{stime.hour:02d}{stime.minute:02d}"
f"{stime.second:05.2f}{latdeg}N{latmin:05.2f}{londeg}E{lonmin:6.2f}{eventsource['depth']:5.2f}"
f"{eventsource['quality']['used_phase_count']:3d}{erh:5.2f}{erz:5.2f}"
f"{eventinfo.magnitudes[0]['mag']:4.2f}{hashID}\n")
# Write phase lines
for key, arrival in arrivals.items():
if 'P' in arrival and arrival['P']['weight'] < 4 and arrival['P']['fm'] is not None:
stat = key
ccode = arrival['P']['channel']
ncode = arrival['P']['network']
Pqual = 'I' if arrival['P']['weight'] < 2 else 'E'
for pick in picks:
if pick.waveform_id.station_code == stat:
resid_picks = pick.get('resource_id')
for origin_arrival in eventinfo.origins[0].arrivals:
if (origin_arrival.get('pick_id') == resid_picks and
origin_arrival.phase == 'P'):
if len(stat) > 4: # HASH handles only 4-character station IDs
stat = stat[1:5]
az = origin_arrival.get('azimuth')
inz = origin_arrival.get('takeoff_angle')
dist = origin_arrival.get('distance')
# Write phase line for HASH-driver 1
fid1.write(f"{stat:<4}{Pqual}P{arrival['P']['fm']}{arrival['P']['weight']:d}"
f"{dist:3.1f}{inz:03d}{az:03d}{ccode}\n")
# Write phase line for HASH-driver 2
fid2.write(f"{stat:<4} {ncode} {ccode} {Pqual} {arrival['P']['fm']}\n")
break
fid1.write(f"{'':<36}{hashID}")
# Prefer Manual Picks over automatic ones if possible
arrivals = chooseArrivals(arrivals) # Function not defined, assumed to exist
if fformat == 'NLLoc':
write_nlloc()
elif fformat == 'HYPO71':
write_hypo71()
elif fformat == 'HYPOSAT':
write_hyposat()
elif fformat == 'VELEST':
write_velest()
elif fformat == 'HYPODD':
write_hypodd()
elif fformat == 'FOCMEC':
write_focmec()
elif fformat == 'HASH': elif fformat == 'HASH':
# two different input files for write_hash()
# HASH-driver 1 and 2 (see HASH manual!)
filename1 = filename + 'drv1' + '.phase'
filename2 = filename + 'drv2' + '.phase'
print("Writing phases to %s for HASH for HASH-driver 1" % filename1)
fid1 = open("%s" % filename1, 'w')
print("Writing phases to %s for HASH for HASH-driver 2" % filename2)
fid2 = open("%s" % filename2, 'w')
# get event information needed for HASH-input file
try:
eventsource = eventinfo.origins[0]
except:
print("No source origin calculated yet, thus no cnv-file creation possible!")
return
eventsource = eventinfo.origins[0]
event = parameter.get('eventID')
hashID = event.split('.')[0][1:5]
latdeg = eventsource['latitude']
latmin = eventsource['latitude'] * 60 / 10000
londeg = eventsource['longitude']
lonmin = eventsource['longitude'] * 60 / 10000
erh = 1 / 2 * (eventsource.origin_uncertainty['min_horizontal_uncertainty'] +
eventsource.origin_uncertainty['max_horizontal_uncertainty']) / 1000
erz = eventsource.depth_errors['uncertainty']
stime = eventsource['time']
if stime.year - 2000 >= 0:
syear = stime.year - 2000
else:
syear = stime.year - 1900
picks = eventinfo.picks
# write header line including event information
# for HASH-driver 1
fid1.write('%s%02d%02d%02d%02d%5.2f%2dN%5.2f%3dE%5.2f%6.3f%4.2f%5.2f%5.2f%s\n' % (syear,
stime.month, stime.day,
stime.hour, stime.minute,
stime.second,
latdeg, latmin, londeg,
lonmin, eventsource['depth'],
eventinfo.magnitudes[0][
'mag'], erh, erz,
hashID))
# write header line including event information
# for HASH-driver 2
fid2.write(
'%d%02d%02d%02d%02d%5.2f%dN%5.2f%3dE%6.2f%5.2f %d %5.2f %5.2f %4.2f %s \n' % (
syear, stime.month, stime.day,
stime.hour, stime.minute, stime.second,
latdeg, latmin, londeg, lonmin,
eventsource['depth'],
eventsource['quality']['used_phase_count'],
erh, erz, eventinfo.magnitudes[0]['mag'],
hashID))
# Prefer Manual Picks over automatic ones if possible
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
# write phase lines
for key in arrivals:
if 'P' in arrivals[key]:
if arrivals[key]['P']['weight'] < 4 and arrivals[key]['P']['fm'] is not None:
stat = key
ccode = arrivals[key]['P']['channel']
ncode = arrivals[key]['P']['network']
if arrivals[key]['P']['weight'] < 2:
Pqual = 'I'
else:
Pqual = 'E'
for i in range(len(picks)):
station = picks[i].waveform_id.station_code
if station == stat:
# get resource ID
resid_picks = picks[i].get('resource_id')
# find same ID in eventinfo
# there it is the pick_id!!
for j in range(len(eventinfo.origins[0].arrivals)):
resid_eventinfo = eventinfo.origins[0].arrivals[j].get('pick_id')
if resid_eventinfo == resid_picks and eventinfo.origins[0].arrivals[j].phase == 'P':
if len(stat) > 4: # HASH handles only 4-string station IDs
stat = stat[1:5]
az = eventinfo.origins[0].arrivals[j].get('azimuth')
inz = eventinfo.origins[0].arrivals[j].get('takeoff_angle')
dist = eventinfo.origins[0].arrivals[j].get('distance')
# write phase line for HASH-driver 1
fid1.write(
'%-4s%sP%s%d 0 %3.1f %03d %03d 2 1 %s\n' % (
stat, Pqual, arrivals[key]['P']['fm'], arrivals[key]['P']['weight'],
dist, inz, az, ccode))
# write phase line for HASH-driver 2
fid2.write('%-4s %s %s %s %s \n' % (
stat,
ncode,
ccode,
Pqual,
arrivals[key]['P']['fm']))
break
fid1.write(' %s' % hashID)
fid1.close()
fid2.close()
def chooseArrivals(arrivals): def chooseArrivals(arrivals):
@@ -1127,7 +670,7 @@ def getQualitiesfromxml(path, errorsP, errorsS, plotflag=1, figure=None, verbosi
mstation = pick.waveform_id.station_code mstation = pick.waveform_id.station_code
mstation_ext = mstation + '_' mstation_ext = mstation + '_'
for mpick in arrivals_copy: for mpick in arrivals_copy:
phase = identifyPhase(loopIdentifyPhase(pick.phase_hint)) # MP MP catch if this fails? phase = identifyPhase(loopIdentifyPhase(pick.phase_hint)) # MP MP catch if this fails?
if ((mpick.waveform_id.station_code == mstation) or if ((mpick.waveform_id.station_code == mstation) or
(mpick.waveform_id.station_code == mstation_ext)) and \ (mpick.waveform_id.station_code == mstation_ext)) and \
(mpick.method_id.id.split('/')[1] == 'auto') and \ (mpick.method_id.id.split('/')[1] == 'auto') and \
+2 -2
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)
+2 -2
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)
+2 -2
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)
+2 -2
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)
+2 -2
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)
+2 -2
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):
+2 -2
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)
+17 -50
View File
@@ -262,10 +262,6 @@ class AutopickStation(object):
self.metadata = metadata self.metadata = metadata
self.origin = origin self.origin = origin
# initialize TauPy pick estimates
self.estFirstP = None
self.estFirstS = None
# initialize picking results # initialize picking results
self.p_results = PickingResults() self.p_results = PickingResults()
self.s_results = PickingResults() self.s_results = PickingResults()
@@ -447,15 +443,15 @@ class AutopickStation(object):
for arr in arrivals: for arr in arrivals:
phases[identifyPhaseID(arr.phase.name)].append(arr) phases[identifyPhaseID(arr.phase.name)].append(arr)
# get first P and S onsets from arrivals list # get first P and S onsets from arrivals list
arrival_time_p = 0 estFirstP = 0
arrival_time_s = 0 estFirstS = 0
if len(phases['P']) > 0: if len(phases['P']) > 0:
arrP, arrival_time_p = min([(arr, arr.time) for arr in phases['P']], key=lambda t: t[1]) arrP, estFirstP = min([(arr, arr.time) for arr in phases['P']], key=lambda t: t[1])
if len(phases['S']) > 0: if len(phases['S']) > 0:
arrS, arrival_time_s = min([(arr, arr.time) for arr in phases['S']], key=lambda t: t[1]) arrS, estFirstS = min([(arr, arr.time) for arr in phases['S']], key=lambda t: t[1])
print('autopick: estimated first arrivals for P: {} s, S:{} s after event' print('autopick: estimated first arrivals for P: {} s, S:{} s after event'
' origin time using TauPy'.format(arrival_time_p, arrival_time_s)) ' origin time using TauPy'.format(estFirstP, estFirstS))
return arrival_time_p, arrival_time_s return estFirstP, estFirstS
def exit_taupy(): def exit_taupy():
"""If taupy failed to calculate theoretical starttimes, picking continues. """If taupy failed to calculate theoretical starttimes, picking continues.
@@ -481,13 +477,10 @@ class AutopickStation(object):
raise AttributeError('No source origins given!') raise AttributeError('No source origins given!')
arrivals = create_arrivals(self.metadata, self.origin, self.pickparams["taup_model"]) arrivals = create_arrivals(self.metadata, self.origin, self.pickparams["taup_model"])
arrival_P, arrival_S = first_PS_onsets(arrivals) estFirstP, estFirstS = first_PS_onsets(arrivals)
self.estFirstP = (self.origin[0].time + arrival_P) - self.ztrace.stats.starttime
# modifiy pstart and pstop relative to estimated first P arrival (relative to station time axis) # modifiy pstart and pstop relative to estimated first P arrival (relative to station time axis)
self.pickparams["pstart"] += self.estFirstP self.pickparams["pstart"] += (self.origin[0].time + estFirstP) - self.ztrace.stats.starttime
self.pickparams["pstop"] += self.estFirstP self.pickparams["pstop"] += (self.origin[0].time + estFirstP) - self.ztrace.stats.starttime
print('autopick: CF calculation times respectively:' print('autopick: CF calculation times respectively:'
' pstart: {} s, pstop: {} s'.format(self.pickparams["pstart"], self.pickparams["pstop"])) ' pstart: {} s, pstop: {} s'.format(self.pickparams["pstart"], self.pickparams["pstop"]))
# make sure pstart and pstop are inside the starttime/endtime of vertical trace # make sure pstart and pstop are inside the starttime/endtime of vertical trace
@@ -498,10 +491,9 @@ class AutopickStation(object):
# for the two horizontal components take earliest and latest time to make sure that the s onset is not clipped # for the two horizontal components take earliest and latest time to make sure that the s onset is not clipped
# if start and endtime of horizontal traces differ, the s windowsize will automatically increase # if start and endtime of horizontal traces differ, the s windowsize will automatically increase
trace_s_start = min([self.etrace.stats.starttime, self.ntrace.stats.starttime]) trace_s_start = min([self.etrace.stats.starttime, self.ntrace.stats.starttime])
self.estFirstS = (self.origin[0].time + arrival_S) - trace_s_start
# modifiy sstart and sstop relative to estimated first S arrival (relative to station time axis) # modifiy sstart and sstop relative to estimated first S arrival (relative to station time axis)
self.pickparams["sstart"] += self.estFirstS self.pickparams["sstart"] += (self.origin[0].time + estFirstS) - trace_s_start
self.pickparams["sstop"] += self.estFirstS self.pickparams["sstop"] += (self.origin[0].time + estFirstS) - trace_s_start
print('autopick: CF calculation times respectively:' print('autopick: CF calculation times respectively:'
' sstart: {} s, sstop: {} s'.format(self.pickparams["sstart"], self.pickparams["sstop"])) ' sstart: {} s, sstop: {} s'.format(self.pickparams["sstart"], self.pickparams["sstop"]))
# make sure pstart and pstop are inside the starttime/endtime of horizontal traces # make sure pstart and pstop are inside the starttime/endtime of horizontal traces
@@ -617,12 +609,6 @@ class AutopickStation(object):
# plot tapered trace filtered with bpz2 filter settings # plot tapered trace filtered with bpz2 filter settings
ax1.plot(tdata, self.tr_filt_z_bpz2.data / max(self.tr_filt_z_bpz2.data), color=linecolor, linewidth=0.7, ax1.plot(tdata, self.tr_filt_z_bpz2.data / max(self.tr_filt_z_bpz2.data), color=linecolor, linewidth=0.7,
label='Data') label='Data')
# plot pickwindows for P
pstart, pstop = self.pickparams['pstart'], self.pickparams['pstop']
if pstart is not None and pstop is not None:
ax1.axvspan(pstart, pstop, color='r', alpha=0.1, zorder=0, label='P window')
if self.estFirstP is not None:
ax1.axvline(self.estFirstP, ls='dashed', color='r', alpha=0.4, label='TauPy estimate')
if self.p_results.weight < 4: if self.p_results.weight < 4:
# plot CF of initial onset (HOScf or ARZcf) # plot CF of initial onset (HOScf or ARZcf)
ax1.plot(self.cf1.getTimeArray(), self.cf1.getCF() / max(self.cf1.getCF()), 'b', label='CF1') ax1.plot(self.cf1.getTimeArray(), self.cf1.getCF() / max(self.cf1.getCF()), 'b', label='CF1')
@@ -727,15 +713,6 @@ class AutopickStation(object):
ax3.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [-1.3, -1.3], 'g', linewidth=2) ax3.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [-1.3, -1.3], 'g', linewidth=2)
ax3.plot([self.s_results.lpp, self.s_results.lpp], [-1.1, 1.1], 'g--', label='lpp') ax3.plot([self.s_results.lpp, self.s_results.lpp], [-1.1, 1.1], 'g--', label='lpp')
ax3.plot([self.s_results.epp, self.s_results.epp], [-1.1, 1.1], 'g--', label='epp') ax3.plot([self.s_results.epp, self.s_results.epp], [-1.1, 1.1], 'g--', label='epp')
# plot pickwindows for S
sstart, sstop = self.pickparams['sstart'], self.pickparams['sstop']
if sstart is not None and sstop is not None:
for axis in [ax2, ax3]:
axis.axvspan(sstart, sstop, color='b', alpha=0.1, zorder=0, label='S window')
if self.estFirstS is not None:
axis.axvline(self.estFirstS, ls='dashed', color='b', alpha=0.4, label='TauPy estimate')
ax3.legend(loc=1) ax3.legend(loc=1)
ax3.set_yticks([]) ax3.set_yticks([])
ax3.set_ylim([-1.5, 1.5]) ax3.set_ylim([-1.5, 1.5])
@@ -858,21 +835,14 @@ class AutopickStation(object):
self.cf1 = None self.cf1 = None
assert isinstance(self.cf1, CharacteristicFunction), 'cf1 is not set correctly: maybe the algorithm name ({})' \ assert isinstance(self.cf1, CharacteristicFunction), 'cf1 is not set correctly: maybe the algorithm name ({})' \
' is corrupted'.format(self.pickparams["algoP"]) ' is corrupted'.format(self.pickparams["algoP"])
# get the original waveform stream from first CF class cut to identical length as CF for plotting
cut_ogstream = self.cf1.getDataArray(self.cf1.getCut())
# MP: Rename to cf_stream for further use of z_copy and to prevent chaos when z_copy suddenly becomes a cf
# stream and later again a waveform stream
cf_stream = z_copy.copy()
cf_stream[0].data = self.cf1.getCF()
# calculate AIC cf from first cf (either HOS or ARZ) # calculate AIC cf from first cf (either HOS or ARZ)
aiccf = AICcf(cf_stream, cuttimes) z_copy[0].data = self.cf1.getCF()
aiccf = AICcf(z_copy, cuttimes)
# get preliminary onset time from AIC-CF # get preliminary onset time from AIC-CF
self.set_current_figure('aicFig') self.set_current_figure('aicFig')
aicpick = AICPicker(aiccf, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot, aicpick = AICPicker(aiccf, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot,
Tsmooth=self.pickparams["aictsmooth"], fig=self.current_figure, Tsmooth=self.pickparams["aictsmooth"], fig=self.current_figure,
linecolor=self.current_linecolor, ogstream=cut_ogstream) linecolor=self.current_linecolor)
# save aicpick for plotting later # save aicpick for plotting later
self.p_data.aicpick = aicpick self.p_data.aicpick = aicpick
# add pstart and pstop to aic plot # add pstart and pstop to aic plot
@@ -885,7 +855,7 @@ class AutopickStation(object):
label='P stop') label='P stop')
ax.legend(loc=1) ax.legend(loc=1)
Pflag = self._pick_p_quality_control(aicpick, cf_stream, tr_filt) Pflag = self._pick_p_quality_control(aicpick, z_copy, tr_filt)
# go on with processing if AIC onset passes quality control # go on with processing if AIC onset passes quality control
slope = aicpick.getSlope() slope = aicpick.getSlope()
if not slope: slope = 0 if not slope: slope = 0
@@ -924,7 +894,7 @@ class AutopickStation(object):
refPpick = PragPicker(self.cf2, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot, refPpick = PragPicker(self.cf2, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot,
self.pickparams["ausP"], self.pickparams["ausP"],
self.pickparams["tsmoothP"], aicpick.getpick(), self.current_figure, self.pickparams["tsmoothP"], aicpick.getpick(), self.current_figure,
self.current_linecolor, ogstream=cut_ogstream) self.current_linecolor)
# save PragPicker result for plotting # save PragPicker result for plotting
self.p_data.refPpick = refPpick self.p_data.refPpick = refPpick
self.p_results.mpp = refPpick.getpick() self.p_results.mpp = refPpick.getpick()
@@ -1176,14 +1146,11 @@ class AutopickStation(object):
# calculate AIC cf # calculate AIC cf
haiccf = self._calculate_aic_cf_s_pick(cuttimesh) haiccf = self._calculate_aic_cf_s_pick(cuttimesh)
# get the original waveform stream cut to identical length as CF for plotting
ogstream = haiccf.getDataArray(haiccf.getCut())
# get preliminary onset time from AIC cf # get preliminary onset time from AIC cf
self.set_current_figure('aicARHfig') self.set_current_figure('aicARHfig')
aicarhpick = AICPicker(haiccf, self.pickparams["tsnrh"], self.pickparams["pickwinS"], self.iplot, aicarhpick = AICPicker(haiccf, self.pickparams["tsnrh"], self.pickparams["pickwinS"], self.iplot,
Tsmooth=self.pickparams["aictsmoothS"], fig=self.current_figure, Tsmooth=self.pickparams["aictsmoothS"], fig=self.current_figure,
linecolor=self.current_linecolor, ogstream=ogstream) linecolor=self.current_linecolor)
# save pick for later plotting # save pick for later plotting
self.aicarhpick = aicarhpick self.aicarhpick = aicarhpick
+43 -15
View File
@@ -17,11 +17,7 @@ autoregressive prediction: application ot local and regional distances, Geophys.
:author: MAGS2 EP3 working group :author: MAGS2 EP3 working group
""" """
import numpy as np import numpy as np
try: from scipy import signal
from scipy.signal import tukey
except ImportError:
from scipy.signal.windows import tukey
from obspy.core import Stream from obspy.core import Stream
from pylot.core.pick.utils import PickingFailedException from pylot.core.pick.utils import PickingFailedException
@@ -60,7 +56,7 @@ class CharacteristicFunction(object):
self.setOrder(order) self.setOrder(order)
self.setFnoise(fnoise) self.setFnoise(fnoise)
self.setARdetStep(t2) self.setARdetStep(t2)
self.calcCF() self.calcCF(self.getDataArray())
self.arpara = np.array([]) self.arpara = np.array([])
self.xpred = np.array([]) self.xpred = np.array([])
@@ -212,15 +208,17 @@ class CharacteristicFunction(object):
data = self.orig_data.copy() data = self.orig_data.copy()
return data return data
def calcCF(self): def calcCF(self, data=None):
pass self.cf = data
class AICcf(CharacteristicFunction): class AICcf(CharacteristicFunction):
def calcCF(self): def calcCF(self, data):
""" """
Function to calculate the Akaike Information Criterion (AIC) after Maeda (1985). Function to calculate the Akaike Information Criterion (AIC) after Maeda (1985).
:param data: data, time series (whether seismogram or CF)
:type data: tuple
:return: AIC function :return: AIC function
:rtype: :rtype:
""" """
@@ -229,7 +227,7 @@ class AICcf(CharacteristicFunction):
ind = np.where(~np.isnan(xnp))[0] ind = np.where(~np.isnan(xnp))[0]
if ind.size: if ind.size:
xnp[:ind[0]] = xnp[ind[0]] xnp[:ind[0]] = xnp[ind[0]]
xnp = tukey(len(xnp), alpha=0.05) * xnp xnp = signal.tukey(len(xnp), alpha=0.05) * xnp
xnp = xnp - np.mean(xnp) xnp = xnp - np.mean(xnp)
datlen = len(xnp) datlen = len(xnp)
k = np.arange(1, datlen) k = np.arange(1, datlen)
@@ -258,11 +256,13 @@ class HOScf(CharacteristicFunction):
""" """
super(HOScf, self).__init__(data, cut, pickparams["tlta"], pickparams["hosorder"]) super(HOScf, self).__init__(data, cut, pickparams["tlta"], pickparams["hosorder"])
def calcCF(self): def calcCF(self, data):
""" """
Function to calculate skewness (statistics of order 3) or kurtosis Function to calculate skewness (statistics of order 3) or kurtosis
(statistics of order 4), using one long moving window, as published (statistics of order 4), using one long moving window, as published
in Kueperkoch et al. (2010), or order 2, i.e. STA/LTA. in Kueperkoch et al. (2010), or order 2, i.e. STA/LTA.
:param data: data, time series (whether seismogram or CF)
:type data: tuple
:return: HOS cf :return: HOS cf
:rtype: :rtype:
""" """
@@ -277,28 +277,47 @@ class HOScf(CharacteristicFunction):
elif self.getOrder() == 4: # this is kurtosis elif self.getOrder() == 4: # this is kurtosis
y = np.power(xnp, 4) y = np.power(xnp, 4)
y1 = np.power(xnp, 2) y1 = np.power(xnp, 2)
elif self.getOrder() == 2: # this is variance, used for STA/LTA processing
y = np.power(xnp, 2)
y1 = np.power(xnp, 2)
# Initialisation # Initialisation
# t2: long term moving window # t2: long term moving window
ilta = int(round(self.getTime2() / self.getIncrement())) ilta = int(round(self.getTime2() / self.getIncrement()))
ista = int(round((self.getTime2() / 10) / self.getIncrement())) # TODO: still hard coded!!
lta = y[0] lta = y[0]
lta1 = y1[0] lta1 = y1[0]
sta = y[0]
# moving windows # moving windows
LTA = np.zeros(len(xnp)) LTA = np.zeros(len(xnp))
STA = np.zeros(len(xnp))
for j in range(0, len(xnp)): for j in range(0, len(xnp)):
if j < 4: if j < 4:
LTA[j] = 0 LTA[j] = 0
STA[j] = 0
elif j <= ista and self.getOrder() == 2:
lta = (y[j] + lta * (j - 1)) / j
if self.getOrder() == 2:
sta = (y[j] + sta * (j - 1)) / j
# elif j < 4:
elif j <= ilta: elif j <= ilta:
lta = (y[j] + lta * (j - 1)) / j lta = (y[j] + lta * (j - 1)) / j
lta1 = (y1[j] + lta1 * (j - 1)) / j lta1 = (y1[j] + lta1 * (j - 1)) / j
if self.getOrder() == 2:
sta = (y[j] - y[j - ista]) / ista + sta
else: else:
lta = (y[j] - y[j - ilta]) / ilta + lta lta = (y[j] - y[j - ilta]) / ilta + lta
lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1 lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1
if self.getOrder() == 2:
sta = (y[j] - y[j - ista]) / ista + sta
# define LTA # define LTA
if self.getOrder() == 3: if self.getOrder() == 3:
LTA[j] = lta / np.power(lta1, 1.5) LTA[j] = lta / np.power(lta1, 1.5)
elif self.getOrder() == 4: elif self.getOrder() == 4:
LTA[j] = lta / np.power(lta1, 2) LTA[j] = lta / np.power(lta1, 2)
else:
LTA[j] = lta
STA[j] = sta
# remove NaN's with first not-NaN-value, # remove NaN's with first not-NaN-value,
# so autopicker doesnt pick discontinuity at start of the trace # so autopicker doesnt pick discontinuity at start of the trace
@@ -307,7 +326,10 @@ class HOScf(CharacteristicFunction):
first = ind[0] first = ind[0]
LTA[:first] = LTA[first] LTA[:first] = LTA[first]
self.cf = LTA if self.getOrder() > 2:
self.cf = LTA
else: # order 2 means STA/LTA!
self.cf = STA / LTA
self.xcf = x self.xcf = x
@@ -317,10 +339,12 @@ class ARZcf(CharacteristicFunction):
super(ARZcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Parorder"], super(ARZcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Parorder"],
fnoise=pickparams["addnoise"]) fnoise=pickparams["addnoise"])
def calcCF(self): def calcCF(self, data):
""" """
function used to calculate the AR prediction error from a single vertical trace. Can be used to pick function used to calculate the AR prediction error from a single vertical trace. Can be used to pick
P onsets. P onsets.
:param data:
:type data: ~obspy.core.stream.Stream
:return: ARZ cf :return: ARZ cf
:rtype: :rtype:
""" """
@@ -451,12 +475,14 @@ class ARHcf(CharacteristicFunction):
super(ARHcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"], super(ARHcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
fnoise=pickparams["addnoise"]) fnoise=pickparams["addnoise"])
def calcCF(self): def calcCF(self, data):
""" """
Function to calculate a characteristic function using autoregressive modelling of the waveform of Function to calculate a characteristic function using autoregressive modelling of the waveform of
both horizontal traces. both horizontal traces.
The waveform is predicted in a moving time window using the calculated AR parameters. The difference The waveform is predicted in a moving time window using the calculated AR parameters. The difference
between the predicted and the actual waveform servers as a characteristic function. between the predicted and the actual waveform servers as a characteristic function.
:param data: wavefor stream
:type data: ~obspy.core.stream.Stream
:return: ARH cf :return: ARH cf
:rtype: :rtype:
""" """
@@ -605,12 +631,14 @@ class AR3Ccf(CharacteristicFunction):
super(AR3Ccf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"], super(AR3Ccf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
fnoise=pickparams["addnoise"]) fnoise=pickparams["addnoise"])
def calcCF(self): def calcCF(self, data):
""" """
Function to calculate a characteristic function using autoregressive modelling of the waveform of Function to calculate a characteristic function using autoregressive modelling of the waveform of
all three traces. all three traces.
The waveform is predicted in a moving time window using the calculated AR parameters. The difference The waveform is predicted in a moving time window using the calculated AR parameters. The difference
between the predicted and the actual waveform servers as a characteristic function between the predicted and the actual waveform servers as a characteristic function
:param data: stream holding all three traces
:type data: ~obspy.core.stream.Stream
:return: AR3C cf :return: AR3C cf
:rtype: :rtype:
""" """
+14 -15
View File
@@ -37,8 +37,7 @@ class AutoPicker(object):
warnings.simplefilter('ignore') warnings.simplefilter('ignore')
def __init__(self, cf, TSNR, PickWindow, iplot=0, aus=None, Tsmooth=None, Pick1=None, def __init__(self, cf, TSNR, PickWindow, iplot=0, aus=None, Tsmooth=None, Pick1=None, fig=None, linecolor='k'):
fig=None, linecolor='k', ogstream=None):
""" """
Create AutoPicker object Create AutoPicker object
:param cf: characteristic function, on which the picking algorithm is applied :param cf: characteristic function, on which the picking algorithm is applied
@@ -60,15 +59,12 @@ class AutoPicker(object):
:type fig: `~matplotlib.figure.Figure` :type fig: `~matplotlib.figure.Figure`
:param linecolor: matplotlib line color string :param linecolor: matplotlib line color string
:type linecolor: str :type linecolor: str
:param ogstream: original stream (waveform), e.g. for plotting purposes
:type ogstream: `~obspy.core.stream.Stream`
""" """
assert isinstance(cf, CharacteristicFunction), "%s is not a CharacteristicFunction object" % str(cf) assert isinstance(cf, CharacteristicFunction), "%s is not a CharacteristicFunction object" % str(cf)
self._linecolor = linecolor self._linecolor = linecolor
self._pickcolor_p = 'b' self._pickcolor_p = 'b'
self.cf = cf.getCF() self.cf = cf.getCF()
self.ogstream = ogstream
self.Tcf = cf.getTimeArray() self.Tcf = cf.getTimeArray()
self.Data = cf.getXCF() self.Data = cf.getXCF()
self.dt = cf.getIncrement() self.dt = cf.getIncrement()
@@ -177,7 +173,7 @@ class AICPicker(AutoPicker):
nn = np.isnan(self.cf) nn = np.isnan(self.cf)
if len(nn) > 1: if len(nn) > 1:
self.cf[nn] = 0 self.cf[nn] = 0
# taper AIC-CF to get rid of side maxima # taper AIC-CF to get rid off side maxima
tap = np.hanning(len(self.cf)) tap = np.hanning(len(self.cf))
aic = tap * self.cf + max(abs(self.cf)) aic = tap * self.cf + max(abs(self.cf))
# smooth AIC-CF # smooth AIC-CF
@@ -320,7 +316,16 @@ class AICPicker(AutoPicker):
plt.close(fig) plt.close(fig)
return return
iislope = islope[0][0:imax + 1] iislope = islope[0][0:imax + 1]
dataslope = self.Data[0].data[iislope] # MP MP change slope calculation
# get all maxima of aicsmooth
iaicmaxima = argrelmax(aicsmooth)[0]
# get first index of maximum after pickindex (indices saved in iaicmaxima)
aicmax = iaicmaxima[np.where(iaicmaxima > pickindex)[0]]
if len(aicmax) > 0:
iaicmax = aicmax[0]
else:
iaicmax = -1
dataslope = aicsmooth[pickindex: iaicmax]
# calculate slope as polynomal fit of order 1 # calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope), 1) xslope = np.arange(0, len(dataslope), 1)
try: try:
@@ -331,7 +336,7 @@ class AICPicker(AutoPicker):
else: else:
self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0]) self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
# normalize slope to maximum of cf to make it unit independent # normalize slope to maximum of cf to make it unit independent
self.slope /= self.Data[0].data[icfmax] self.slope /= aicsmooth[iaicmax]
except Exception as e: except Exception as e:
print("AICPicker: Problems with data fitting! {}".format(e)) print("AICPicker: Problems with data fitting! {}".format(e))
@@ -351,12 +356,6 @@ class AICPicker(AutoPicker):
self.Tcf = self.Tcf[0:len(self.Tcf) - 1] self.Tcf = self.Tcf[0:len(self.Tcf) - 1]
ax1.plot(self.Tcf, cf / max(cf), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data') ax1.plot(self.Tcf, cf / max(cf), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
ax1.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF') ax1.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
# plot the original waveform also for evaluation of the CF and pick
if self.ogstream:
data = self.ogstream[0].data
if len(data) == len(self.Tcf):
ax1.plot(self.Tcf, 0.5 * data / max(data), 'k', label='Seismogram', alpha=0.3, zorder=0,
lw=0.5)
if self.Pick is not None: if self.Pick is not None:
ax1.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b', linewidth=2, label='AIC-Pick') ax1.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b', linewidth=2, label='AIC-Pick')
ax1.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime) ax1.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
@@ -377,7 +376,7 @@ class AICPicker(AutoPicker):
label='Signal Window') label='Signal Window')
ax2.axvspan(self.Tcf[iislope[0]], self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0, ax2.axvspan(self.Tcf[iislope[0]], self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0,
label='Slope Window') label='Slope Window')
ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2, ax2.plot(self.Tcf[pickindex: iaicmax], datafit, 'g', linewidth=2,
label='Slope') # MP MP changed temporarily! label='Slope') # MP MP changed temporarily!
if self.slope is not None: if self.slope is not None:
+5 -18
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, common_range 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'):
@@ -828,22 +828,14 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
if len(X) > 1: if len(X) > 1:
# all three components available # all three components available
# make sure, all components have equal lengths # make sure, all components have equal lengths
earliest_starttime = min(tr.stats.starttime for tr in X) ilen = min([len(X[0].data), len(X[1].data), len(X[2].data)])
cuttimes = common_range(X) x1 = X[0][0:ilen]
X = X.slice(cuttimes[0], cuttimes[1]) x2 = X[1][0:ilen]
x1, x2, x3 = X[:3] x3 = X[2][0:ilen]
if not (len(x1) == len(x2) == len(x3)):
raise PickingFailedException('checksignallength: unequal lengths of components!')
# get RMS trace # get RMS trace
rms = np.sqrt((np.power(x1, 2) + np.power(x2, 2) + np.power(x3, 2)) / 3) rms = np.sqrt((np.power(x1, 2) + np.power(x2, 2) + np.power(x3, 2)) / 3)
ilen = len(rms)
dt = earliest_starttime - X[0].stats.starttime
pick -= dt
else: else:
x1 = X[0].data x1 = X[0].data
x2 = x3 = None
ilen = len(x1) ilen = len(x1)
rms = abs(x1) rms = abs(x1)
@@ -882,10 +874,6 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
fig._tight = True fig._tight = True
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
ax.plot(t, rms, color=linecolor, linewidth=0.7, label='RMS Data') ax.plot(t, rms, color=linecolor, linewidth=0.7, label='RMS Data')
ax.plot(t, x1, 'k', alpha=0.3, lw=0.3, zorder=0)
if x2 is not None and x3 is not None:
ax.plot(t, x2, 'r', alpha=0.3, lw=0.3, zorder=0)
ax.plot(t, x3, 'g', alpha=0.3, lw=0.3, zorder=0)
ax.axvspan(t[inoise[0]], t[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window') ax.axvspan(t[inoise[0]], t[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
ax.axvspan(t[isignal[0]], t[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window') ax.axvspan(t[isignal[0]], t[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window')
ax.plot([t[isignal[0]], t[isignal[len(isignal) - 1]]], ax.plot([t[isignal[0]], t[isignal[len(isignal) - 1]]],
@@ -895,7 +883,6 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
ax.set_xlabel('Time [s] since %s' % X[0].stats.starttime) ax.set_xlabel('Time [s] since %s' % X[0].stats.starttime)
ax.set_ylabel('Counts') ax.set_ylabel('Counts')
ax.set_title('Check for Signal Length, Station %s' % X[0].stats.station) ax.set_title('Check for Signal Length, Station %s' % X[0].stats.station)
ax.set_xlim(pickparams["pstart"], pickparams["pstop"])
ax.set_yticks([]) ax.set_yticks([])
if plt_flag == 1: if plt_flag == 1:
fig.show() fig.show()
+107 -225
View File
@@ -1,23 +1,18 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import copy
import traceback import traceback
import cartopy.crs as ccrs import cartopy.crs as ccrs
import cartopy.feature as cf import cartopy.feature as cf
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib import matplotlib
import matplotlib.patheffects as PathEffects import matplotlib.patheffects as PathEffects
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import obspy import obspy
from PySide2 import QtWidgets, QtGui from PySide2 import QtWidgets
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from mpl_toolkits.axes_grid1.inset_locator import inset_axes from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from obspy import UTCDateTime
from pylot.core.io.data import Data
from pylot.core.io.phases import picksdict_from_picks
from pylot.core.util.utils import identifyPhaseID from pylot.core.util.utils import identifyPhaseID
from scipy.interpolate import griddata from scipy.interpolate import griddata
@@ -29,10 +24,10 @@ matplotlib.use('Qt5Agg')
class MplCanvas(FigureCanvas): class MplCanvas(FigureCanvas):
def __init__(self, extern_axes=None, projection=None, width=15, height=5, dpi=100): def __init__(self, parent=None, extern_axes=None, width=5, height=4, dpi=100):
if extern_axes is None: if extern_axes is None:
self.fig = plt.figure(figsize=(width, height), dpi=dpi) self.fig = plt.figure(figsize=(width, height), dpi=dpi)
self.axes = self.fig.add_subplot(111, projection=projection) self.axes = self.fig.add_subplot(111)
else: else:
self.fig = extern_axes.figure self.fig = extern_axes.figure
self.axes = extern_axes self.axes = extern_axes
@@ -43,7 +38,6 @@ class MplCanvas(FigureCanvas):
class Array_map(QtWidgets.QWidget): class Array_map(QtWidgets.QWidget):
def __init__(self, parent, metadata, parameter=None, axes=None, annotate=True, pointsize=25., def __init__(self, parent, metadata, parameter=None, axes=None, annotate=True, pointsize=25.,
linewidth=1.5, width=5e6, height=2e6): linewidth=1.5, width=5e6, height=2e6):
QtWidgets.QWidget.__init__(self, parent=parent) QtWidgets.QWidget.__init__(self, parent=parent)
assert (parameter is not None or parent is not None), 'either parent or parameter has to be set' assert (parameter is not None or parent is not None), 'either parent or parameter has to be set'
@@ -61,36 +55,28 @@ class Array_map(QtWidgets.QWidget):
self.uncertainties = None self.uncertainties = None
self.autopicks_dict = None self.autopicks_dict = None
self.hybrids_dict = None self.hybrids_dict = None
self.picks_dict_ref = None
self.autopicks_dict_ref = None
self.eventLoc = None self.eventLoc = None
self.parameter = parameter if parameter else parent._inputs self.parameter = parameter if parameter else parent._inputs
self.picks_rel = {} self.picks_rel = {}
self.picks_rel_mean_corrected = {}
self.marked_stations = [] self.marked_stations = []
self.highlighted_stations = [] self.highlighted_stations = []
# call functions to draw everything # call functions to draw everything
self.projection = ccrs.PlateCarree()
self.init_graphics() self.init_graphics()
self.ax = self.canvas.axes
self.ax.set_adjustable('datalim')
self.init_stations() self.init_stations()
self.init_crtpyMap() self.init_crtpyMap()
self.init_map() self.init_map()
# set original map limits to fall back on when home button is pressed # set original map limits to fall back on when home button is pressed
self.org_xlim = self.ax.get_xlim() self.org_xlim = self.canvas.axes.get_xlim()
self.org_ylim = self.ax.get_ylim() self.org_ylim = self.canvas.axes.get_ylim()
# initial map without event # initial map without event
self.ax.set_xlim(self.org_xlim[0], self.org_xlim[1]) self.canvas.axes.set_xlim(self.org_xlim[0], self.org_xlim[1])
self.ax.set_ylim(self.org_ylim[0], self.org_ylim[1]) self.canvas.axes.set_ylim(self.org_ylim[0], self.org_ylim[1])
self._style = None if not hasattr(parent, '_style') else parent._style self._style = None if not hasattr(parent, '_style') else parent._style
def init_map(self): def init_map(self):
self.init_colormap() self.init_colormap()
self.connectSignals() self.connectSignals()
@@ -103,26 +89,23 @@ class Array_map(QtWidgets.QWidget):
# initialize figure elements # initialize figure elements
if self.extern_plot_axes is None: if self.extern_plot_axes is None:
self.canvas = MplCanvas(projection=self.projection) self.canvas = MplCanvas(self)
self.plotWidget = FigureCanvas(self.canvas.fig)
else: else:
self.canvas = MplCanvas(extern_axes=self.extern_plot_axes) self.canvas = MplCanvas(self, extern_axes=self.extern_plot_axes)
self.plotWidget = FigureCanvas(self.canvas.fig)
self.plotWidget = self.canvas
# initialize GUI elements # initialize GUI elements
self.status_label = QtWidgets.QLabel() self.status_label = QtWidgets.QLabel()
self.map_reset_button = QtWidgets.QPushButton('Reset Map View') self.map_reset_button = QtWidgets.QPushButton('Reset Map View')
self.save_map_button = QtWidgets.QPushButton('Save Map') self.save_map_button = QtWidgets.QPushButton('Save Map')
self.load_reference_picks_button = QtWidgets.QPushButton('Load reference picks.')
self.go2eq_button = QtWidgets.QPushButton('Go to Event Location') self.go2eq_button = QtWidgets.QPushButton('Go to Event Location')
self.subtract_mean_cb = QtWidgets.QCheckBox('Subtract mean')
self.subtract_ref_cb = QtWidgets.QCheckBox('Subtract reference onsets')
self.main_box = QtWidgets.QVBoxLayout() self.main_box = QtWidgets.QVBoxLayout()
self.setLayout(self.main_box) self.setLayout(self.main_box)
self.top_row = QtWidgets.QHBoxLayout() self.top_row = QtWidgets.QHBoxLayout()
self.main_box.addLayout(self.top_row, 0) self.main_box.addLayout(self.top_row, 1)
self.comboBox_phase = QtWidgets.QComboBox() self.comboBox_phase = QtWidgets.QComboBox()
self.comboBox_phase.insertItem(0, 'P') self.comboBox_phase.insertItem(0, 'P')
@@ -155,17 +138,14 @@ class Array_map(QtWidgets.QWidget):
self.top_row.addWidget(self.auto_refresh_box) self.top_row.addWidget(self.auto_refresh_box)
self.top_row.addWidget(self.refresh_button) self.top_row.addWidget(self.refresh_button)
self.main_box.addWidget(self.plotWidget, 10) self.main_box.addWidget(self.plotWidget, 1)
self.bot_row = QtWidgets.QHBoxLayout() self.bot_row = QtWidgets.QHBoxLayout()
self.main_box.addLayout(self.bot_row, 0) self.main_box.addLayout(self.bot_row, 0.3)
self.bot_row.addWidget(QtWidgets.QLabel(''), 5) self.bot_row.addWidget(QtWidgets.QLabel(''), 5)
self.bot_row.addWidget(self.map_reset_button, 2) self.bot_row.addWidget(self.map_reset_button, 2)
self.bot_row.addWidget(self.go2eq_button, 2) self.bot_row.addWidget(self.go2eq_button, 2)
self.bot_row.addWidget(self.save_map_button, 2) self.bot_row.addWidget(self.save_map_button, 2)
self.bot_row.addWidget(self.load_reference_picks_button, 2)
self.bot_row.addWidget(self.subtract_mean_cb, 0)
self.bot_row.addWidget(self.subtract_ref_cb, 0)
self.bot_row.addWidget(self.status_label, 5) self.bot_row.addWidget(self.status_label, 5)
def init_colormap(self): def init_colormap(self):
@@ -173,12 +153,14 @@ class Array_map(QtWidgets.QWidget):
self.init_lat_lon_grid() self.init_lat_lon_grid()
def init_crtpyMap(self): def init_crtpyMap(self):
self.ax.add_feature(cf.LAND) self.canvas.axes.cla()
self.ax.add_feature(cf.OCEAN) self.canvas.axes = plt.axes(projection=ccrs.PlateCarree())
self.ax.add_feature(cf.COASTLINE, linewidth=1, edgecolor='gray') self.canvas.axes.add_feature(cf.LAND)
self.ax.add_feature(cf.BORDERS, alpha=0.7) self.canvas.axes.add_feature(cf.OCEAN)
self.ax.add_feature(cf.LAKES, alpha=0.7) self.canvas.axes.add_feature(cf.COASTLINE, linewidth=1, edgecolor='gray')
self.ax.add_feature(cf.RIVERS, linewidth=1) self.canvas.axes.add_feature(cf.BORDERS, alpha=0.7)
self.canvas.axes.add_feature(cf.LAKES, alpha=0.7)
self.canvas.axes.add_feature(cf.RIVERS, linewidth=1)
# parallels and meridians # parallels and meridians
self.add_merid_paral() self.add_merid_paral()
@@ -186,8 +168,12 @@ class Array_map(QtWidgets.QWidget):
self.canvas.fig.tight_layout() self.canvas.fig.tight_layout()
def add_merid_paral(self): def add_merid_paral(self):
self.gridlines = self.ax.gridlines(draw_labels=False, alpha=0.6, color='gray', self.gridlines = self.canvas.axes.gridlines(draw_labels=False, alpha=0.6, color='gray',
linewidth=self.linewidth / 2, zorder=7, crs=ccrs.PlateCarree()) linewidth=self.linewidth / 2, zorder=7)
# TODO: current cartopy version does not support label removal. Devs are working on it.
# Should be fixed in coming cartopy versions
# self.gridlines.xformatter = LONGITUDE_FORMATTER
# self.gridlines.yformatter = LATITUDE_FORMATTER
def remove_merid_paral(self): def remove_merid_paral(self):
if len(self.gridlines.xline_artists): if len(self.gridlines.xline_artists):
@@ -195,24 +181,24 @@ class Array_map(QtWidgets.QWidget):
self.gridlines.yline_artists[0].remove() self.gridlines.yline_artists[0].remove()
def org_map_view(self): def org_map_view(self):
self.ax.set_xlim(self.org_xlim[0], self.org_xlim[1]) self.canvas.axes.set_xlim(self.org_xlim[0], self.org_xlim[1])
self.ax.set_ylim(self.org_ylim[0], self.org_ylim[1]) self.canvas.axes.set_ylim(self.org_ylim[0], self.org_ylim[1])
# parallels and meridians # parallels and meridians
#self.remove_merid_paral() self.remove_merid_paral()
#self.add_merid_paral() self.add_merid_paral()
self.canvas.draw_idle() self.canvas.axes.figure.canvas.draw_idle()
def go2eq(self): def go2eq(self):
if self.eventLoc: if self.eventLoc:
lats, lons = self.eventLoc lats, lons = self.eventLoc
self.ax.set_xlim(lons - 10, lons + 10) self.canvas.axes.set_xlim(lons - 10, lons + 10)
self.ax.set_ylim(lats - 5, lats + 5) self.canvas.axes.set_ylim(lats - 5, lats + 5)
# parallels and meridians # parallels and meridians
#self.remove_merid_paral() self.remove_merid_paral()
#self.add_merid_paral() self.add_merid_paral()
self.canvas.draw_idle() self.canvas.axes.figure.canvas.draw_idle()
else: else:
self.status_label.setText('No event information available') self.status_label.setText('No event information available')
@@ -226,9 +212,6 @@ class Array_map(QtWidgets.QWidget):
self.map_reset_button.clicked.connect(self.org_map_view) self.map_reset_button.clicked.connect(self.org_map_view)
self.go2eq_button.clicked.connect(self.go2eq) self.go2eq_button.clicked.connect(self.go2eq)
self.save_map_button.clicked.connect(self.saveFigure) self.save_map_button.clicked.connect(self.saveFigure)
self.load_reference_picks_button.clicked.connect(self.load_reference_picks)
self.subtract_mean_cb.stateChanged.connect(self.toggle_subtract_mean)
self.subtract_ref_cb.stateChanged.connect(self.toggle_subtract_ref)
self.plotWidget.mpl_connect('motion_notify_event', self.mouse_moved) self.plotWidget.mpl_connect('motion_notify_event', self.mouse_moved)
self.plotWidget.mpl_connect('scroll_event', self.mouse_scroll) self.plotWidget.mpl_connect('scroll_event', self.mouse_scroll)
@@ -237,32 +220,21 @@ class Array_map(QtWidgets.QWidget):
# set mouse events ----------------------------------------------------- # set mouse events -----------------------------------------------------
def mouse_moved(self, event): def mouse_moved(self, event):
if not event.inaxes == self.ax: if not event.inaxes == self.canvas.axes:
return return
else:
cont, inds = self.sc.contains(event)
lat = event.ydata lat = event.ydata
lon = event.xdata lon = event.xdata
text = f'Longitude: {lon:3.3f}, Latitude: {lat:3.3f}' self.status_label.setText('Latitude: {:3.5f}, Longitude: {:3.5f}'.format(lat, lon))
if cont:
indices = inds['ind']
text += ' | Station: ' if len(indices) == 1 else ' | Stations: '
text += ' - '.join([self._station_onpick_ids[index] for index in indices[:5]])
if len(indices) > 5:
text += '...'
self.status_label.setText(text)
def mouse_scroll(self, event): def mouse_scroll(self, event):
if not event.inaxes == self.ax: if not event.inaxes == self.canvas.axes:
return return
zoom = {'up': 1. / 2., 'down': 2.} zoom = {'up': 1. / 2., 'down': 2.}
if event.button in zoom: if event.button in zoom:
xlim = self.ax.get_xlim() xlim = self.canvas.axes.get_xlim()
ylim = self.ax.get_ylim() ylim = self.canvas.axes.get_ylim()
x, y = event.xdata, event.ydata x, y = event.xdata, event.ydata
@@ -274,24 +246,24 @@ class Array_map(QtWidgets.QWidget):
yb = y - 0.5 * ydiff yb = y - 0.5 * ydiff
yt = y + 0.5 * ydiff yt = y + 0.5 * ydiff
self.ax.set_xlim(xl, xr) self.canvas.axes.set_xlim(xl, xr)
self.ax.set_ylim(yb, yt) self.canvas.axes.set_ylim(yb, yt)
# parallels and meridians # parallels and meridians
#self.remove_merid_paral() self.remove_merid_paral()
#self.add_merid_paral() self.add_merid_paral()
self.ax.figure.canvas.draw_idle() self.canvas.axes.figure.canvas.draw_idle()
def mouseLeftPress(self, event): def mouseLeftPress(self, event):
if not event.inaxes == self.ax: if not event.inaxes == self.canvas.axes:
return return
self.map_x = event.xdata self.map_x = event.xdata
self.map_y = event.ydata self.map_y = event.ydata
self.map_xlim = self.ax.get_xlim() self.map_xlim = self.canvas.axes.get_xlim()
self.map_ylim = self.ax.get_ylim() self.map_ylim = self.canvas.axes.get_ylim()
def mouseLeftRelease(self, event): def mouseLeftRelease(self, event):
if not event.inaxes == self.ax: if not event.inaxes == self.canvas.axes:
return return
new_x = event.xdata new_x = event.xdata
new_y = event.ydata new_y = event.ydata
@@ -299,13 +271,13 @@ class Array_map(QtWidgets.QWidget):
dx = new_x - self.map_x dx = new_x - self.map_x
dy = new_y - self.map_y dy = new_y - self.map_y
self.ax.set_xlim((self.map_xlim[0] - dx, self.map_xlim[1] - dx)) self.canvas.axes.set_xlim((self.map_xlim[0] - dx, self.map_xlim[1] - dx))
self.ax.set_ylim(self.map_ylim[0] - dy, self.map_ylim[1] - dy) self.canvas.axes.set_ylim(self.map_ylim[0] - dy, self.map_ylim[1] - dy)
# parallels and meridians # parallels and meridians
#self.remove_merid_paral() self.remove_merid_paral()
#self.add_merid_paral() self.add_merid_paral()
self.ax.figure.canvas.draw_idle() self.canvas.axes.figure.canvas.draw_idle()
def onpick(self, event): def onpick(self, event):
btn_msg = {1: ' in selection. Aborted', 2: ' to delete a pick on. Aborted', 3: ' to display info.'} btn_msg = {1: ' in selection. Aborted', 2: ' to delete a pick on. Aborted', 3: ' to display info.'}
@@ -385,11 +357,14 @@ class Array_map(QtWidgets.QWidget):
def get_max_from_stations(self, key): def get_max_from_stations(self, key):
return self._from_dict(max, key) return self._from_dict(max, key)
def get_selected_pick_type(self): def get_min_from_picks(self):
return self.comboBox_am.currentText().split(' ')[0] return min(self.picks_rel.values())
def get_max_from_picks(self):
return max(self.picks_rel.values())
def current_picks_dict(self): def current_picks_dict(self):
picktype = self.get_selected_pick_type() picktype = self.comboBox_am.currentText().split(' ')[0]
auto_manu = {'auto': self.autopicks_dict, auto_manu = {'auto': self.autopicks_dict,
'manual': self.picks_dict, 'manual': self.picks_dict,
'hybrid': self.hybrids_dict} 'hybrid': self.hybrids_dict}
@@ -432,37 +407,22 @@ class Array_map(QtWidgets.QWidget):
print('Cannot display pick for station {}. Reason: {}'.format(station_name, e)) print('Cannot display pick for station {}. Reason: {}'.format(station_name, e))
return picks, uncertainties return picks, uncertainties
def get_picks_rel(picks, func=min): def get_picks_rel(picks):
picks_rel = {} picks_rel = {}
picks_utc = [] picks_utc = []
for pick in picks.values(): for pick in picks.values():
if type(pick) is UTCDateTime: if type(pick) is obspy.core.utcdatetime.UTCDateTime:
picks_utc.append(pick.timestamp) picks_utc.append(pick)
if picks_utc: if picks_utc:
self._reference_picktime = UTCDateTime(func(picks_utc)) self._earliest_picktime = min(picks_utc)
for st_id, pick in picks.items(): for st_id, pick in picks.items():
if type(pick) is UTCDateTime: if type(pick) is obspy.core.utcdatetime.UTCDateTime:
pick -= self._reference_picktime pick -= self._earliest_picktime
picks_rel[st_id] = pick picks_rel[st_id] = pick
return picks_rel return picks_rel
def get_picks_rel_mean_corr(picks):
return get_picks_rel(picks, func=np.nanmean)
self.picks, self.uncertainties = get_picks(self.stations_dict) self.picks, self.uncertainties = get_picks(self.stations_dict)
self.picks_rel = get_picks_rel(self.picks) self.picks_rel = get_picks_rel(self.picks)
self.picks_rel_mean_corrected = get_picks_rel_mean_corr(self.picks)
def toggle_subtract_mean(self):
if self.subtract_mean_cb.isChecked():
cmap = 'seismic'
else:
cmap = 'viridis'
self.cmaps_box.setCurrentIndex(self.cmaps_box.findText(cmap))
self._refresh_drawings()
def toggle_subtract_ref(self):
self._refresh_drawings()
def init_lat_lon_dimensions(self): def init_lat_lon_dimensions(self):
# init minimum and maximum lon and lat dimensions # init minimum and maximum lon and lat dimensions
@@ -476,10 +436,7 @@ class Array_map(QtWidgets.QWidget):
self.longrid, self.latgrid = np.meshgrid(lonaxis, lataxis) self.longrid, self.latgrid = np.meshgrid(lonaxis, lataxis)
def init_picksgrid(self): def init_picksgrid(self):
rval = self.get_picks_lat_lon() picks, uncertainties, lats, lons = self.get_picks_lat_lon()
if not rval:
return
picks, uncertainties, lats, lons = rval
try: try:
self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid), method='linear') self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid), method='linear')
except Exception as e: except Exception as e:
@@ -496,21 +453,15 @@ class Array_map(QtWidgets.QWidget):
return stations, latitudes, longitudes return stations, latitudes, longitudes
def get_picks_lat_lon(self): def get_picks_lat_lon(self):
picks_rel = self.picks_rel_mean_corrected if self.subtract_mean_cb.isChecked() else self.picks_rel
picks_rel = self.picks_rel_ref_corrected if self.subtract_ref_cb.isChecked() else picks_rel
if not picks_rel:
return
picks = [] picks = []
uncertainties = [] uncertainties = []
latitudes = [] latitudes = []
longitudes = [] longitudes = []
for nwst_id, pick in picks_rel.items(): for st_id, pick in self.picks_rel.items():
if nwst_id not in self.uncertainties or nwst_id not in self.stations_dict:
continue
picks.append(pick) picks.append(pick)
uncertainties.append(self.uncertainties.get(nwst_id)) uncertainties.append(self.uncertainties.get(st_id))
latitudes.append(self.stations_dict[nwst_id]['latitude']) latitudes.append(self.stations_dict[st_id]['latitude'])
longitudes.append(self.stations_dict[nwst_id]['longitude']) longitudes.append(self.stations_dict[st_id]['longitude'])
return picks, uncertainties, latitudes, longitudes return picks, uncertainties, latitudes, longitudes
# plotting ----------------------------------------------------- # plotting -----------------------------------------------------
@@ -518,20 +469,13 @@ class Array_map(QtWidgets.QWidget):
stat_dict = self.stations_dict['{}.{}'.format(network, station)] stat_dict = self.stations_dict['{}.{}'.format(network, station)]
lat = stat_dict['latitude'] lat = stat_dict['latitude']
lon = stat_dict['longitude'] lon = stat_dict['longitude']
self.highlighted_stations.append(self.ax.scatter(lon, lat, s=self.pointsize, edgecolors=color, self.highlighted_stations.append(self.canvas.axes.scatter(lon, lat, s=self.pointsize, edgecolors=color,
facecolors='none', zorder=12, facecolors='none', zorder=12,
transform=ccrs.PlateCarree(), label='deleted')) transform=ccrs.PlateCarree(), label='deleted'))
def openPickDlg(self, ind): def openPickDlg(self, ind):
try: wfdata = self._parent.get_data().getWFData()
wfdata = self._parent.get_data().getWFData() wfdata_comp = self._parent.get_data().getWFDataComp()
except AttributeError:
QtWidgets.QMessageBox.warning(
self, "PyLoT Warning",
"No waveform data found. Check if they were already loaded in Waveform plot tab."
)
return
wfdata_comp = self._parent.get_data().getAltWFdata()
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
@@ -546,6 +490,7 @@ class Array_map(QtWidgets.QWidget):
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,
model=self.parameter.get('taup_model'),
event=pyl_mw.get_current_event()) event=pyl_mw.get_current_event())
except Exception as e: except Exception as e:
message = 'Could not generate Plot for station {st}.\n {er}'.format(st=station, er=e) message = 'Could not generate Plot for station {st}.\n {er}'.format(st=station, er=e)
@@ -573,19 +518,12 @@ class Array_map(QtWidgets.QWidget):
print(message, e) print(message, e)
print(traceback.format_exc()) print(traceback.format_exc())
def draw_contour_filled(self, nlevel=51): def draw_contour_filled(self, nlevel=50):
if self.subtract_mean_cb.isChecked(): levels = np.linspace(self.get_min_from_picks(), self.get_max_from_picks(), nlevel)
abs_max = self.get_residuals_absmax()
levels = np.linspace(-abs_max, abs_max, nlevel)
else:
levels = np.linspace(min(self.picks_rel.values()), max(self.picks_rel.values()), nlevel)
self.contourf = self.ax.contourf(self.longrid, self.latgrid, self.picksgrid_active, levels, self.contourf = self.canvas.axes.contourf(self.longrid, self.latgrid, self.picksgrid_active, levels,
linewidths=self.linewidth * 5, transform=ccrs.PlateCarree(), linewidths=self.linewidth * 5, transform=ccrs.PlateCarree(),
alpha=0.4, zorder=8, cmap=self.get_colormap()) alpha=0.4, zorder=8, cmap=self.get_colormap())
def get_residuals_absmax(self):
return np.max(np.absolute(list(self.picks_rel_mean_corrected.values())))
def get_colormap(self): def get_colormap(self):
return plt.get_cmap(self.cmaps_box.currentText()) return plt.get_cmap(self.cmaps_box.currentText())
@@ -593,24 +531,21 @@ class Array_map(QtWidgets.QWidget):
def scatter_all_stations(self): def scatter_all_stations(self):
stations, lats, lons = self.get_st_lat_lon_for_plot() stations, lats, lons = self.get_st_lat_lon_for_plot()
self.sc = self.ax.scatter(lons, lats, s=self.pointsize * 3, facecolor='none', marker='.', self.sc = self.canvas.axes.scatter(lons, lats, s=self.pointsize * 3, facecolor='none', marker='.',
zorder=10, picker=True, edgecolor='0.5', label='Not Picked', zorder=10, picker=True, edgecolor='0.5', label='Not Picked',
transform=ccrs.PlateCarree()) transform=ccrs.PlateCarree())
self.cid = self.plotWidget.mpl_connect('pick_event', self.onpick) self.cid = self.plotWidget.mpl_connect('pick_event', self.onpick)
self._station_onpick_ids = stations self._station_onpick_ids = stations
if self.eventLoc: if self.eventLoc:
lats, lons = self.eventLoc lats, lons = self.eventLoc
self.sc_event = self.ax.scatter(lons, lats, s=5 * self.pointsize, facecolor='red', zorder=11, self.sc_event = self.canvas.axes.scatter(lons, lats, s=5 * self.pointsize, facecolor='red', zorder=11,
label='Event (might be outside map region)', marker='*', label='Event (might be outside map region)', marker='*',
edgecolors='black', edgecolors='black',
transform=ccrs.PlateCarree()) transform=ccrs.PlateCarree())
def scatter_picked_stations(self): def scatter_picked_stations(self):
rval = self.get_picks_lat_lon() picks, uncertainties, lats, lons = self.get_picks_lat_lon()
if not rval:
return
picks, uncertainties, lats, lons = rval
if len(lons) < 1 and len(lats) < 1: if len(lons) < 1 and len(lats) < 1:
return return
@@ -620,13 +555,8 @@ class Array_map(QtWidgets.QWidget):
for uncertainty in uncertainties]) for uncertainty in uncertainties])
cmap = self.get_colormap() cmap = self.get_colormap()
self.sc_picked = self.canvas.axes.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap,
vmin = vmax = None c=picks, zorder=11, label='Picked', transform=ccrs.PlateCarree())
if self.subtract_mean_cb.isChecked():
vmin, vmax = -self.get_residuals_absmax(), self.get_residuals_absmax()
self.sc_picked = self.ax.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap, vmin=vmin, vmax=vmax,
c=picks, zorder=11, label='Picked', transform=ccrs.PlateCarree())
def annotate_ax(self): def annotate_ax(self):
self.annotations = [] self.annotations = []
@@ -644,20 +574,20 @@ class Array_map(QtWidgets.QWidget):
if st in self.marked_stations: if st in self.marked_stations:
color = 'red' color = 'red'
self.annotations.append( self.annotations.append(
self.ax.annotate(f'{st}', xy=(x + 0.003, y + 0.003), fontsize=self.pointsize / 4., self.canvas.axes.annotate(' %s' % st, xy=(x + 0.003, y + 0.003), fontsize=self.pointsize / 4.,
fontweight='semibold', color=color, alpha=0.8, fontweight='semibold', color=color, alpha=0.8,
transform=ccrs.PlateCarree(), zorder=14, transform=ccrs.PlateCarree(), zorder=14,
path_effects=[PathEffects.withStroke( path_effects=[PathEffects.withStroke(
linewidth=self.pointsize / 15., foreground='k')])) linewidth=self.pointsize / 15., foreground='k')]))
self.legend = self.ax.legend(loc=1, framealpha=1) self.legend = self.canvas.axes.legend(loc=1, framealpha=1)
self.legend.set_zorder(100) self.legend.set_zorder(100)
self.legend.get_frame().set_facecolor((1, 1, 1, 0.95)) self.legend.get_frame().set_facecolor((1, 1, 1, 0.95))
def add_cbar(self, label): def add_cbar(self, label):
self.cbax_bg = inset_axes(self.ax, width="6%", height="75%", loc=5) self.cbax_bg = inset_axes(self.canvas.axes, width="6%", height="75%", loc=5)
cbax = inset_axes(self.ax, width='2%', height='70%', loc=5) cbax = inset_axes(self.canvas.axes, width='2%', height='70%', loc=5)
cbar = self.ax.figure.colorbar(self.sc_picked, cax=cbax) cbar = self.canvas.axes.figure.colorbar(self.sc_picked, cax=cbax)
cbar.set_label(label) cbar.set_label(label)
cbax.yaxis.tick_left() cbax.yaxis.tick_left()
cbax.yaxis.set_label_position('left') cbax.yaxis.set_label_position('left')
@@ -702,9 +632,7 @@ class Array_map(QtWidgets.QWidget):
if picks_available: if picks_available:
self.scatter_picked_stations() self.scatter_picked_stations()
if hasattr(self, 'sc_picked'): if hasattr(self, 'sc_picked'):
self.cbar = self.add_cbar( self.cbar = self.add_cbar(label='Time relative to first onset ({}) [s]'.format(self._earliest_picktime))
label='Time relative to reference onset ({}) [s]'.format(self._reference_picktime)
)
self.comboBox_phase.setEnabled(True) self.comboBox_phase.setEnabled(True)
else: else:
self.comboBox_phase.setEnabled(False) self.comboBox_phase.setEnabled(False)
@@ -765,52 +693,6 @@ class Array_map(QtWidgets.QWidget):
fname += '.png' fname += '.png'
self.canvas.fig.savefig(fname) self.canvas.fig.savefig(fname)
def load_reference_picks(self):
fname = self._parent.filename_from_action()
data_ref = Data(parent=self._parent, evtdata=str(fname))
evt_ref = data_ref.get_evt_data()
if not evt_ref:
return
picks_ref = evt_ref.picks
if not picks_ref:
return
picksdict_ref = picksdict_from_picks(evt_ref, nwst_id=False)
self.autopicks_dict = picksdict_ref['manual']
self.autopicks_dict_ref = picksdict_ref['auto']
return True
@property
def picks_rel_ref_corrected(self):
picktype = self.get_selected_pick_type()
if picktype in ['auto', 'hybrid']:
picks_dict = self.autopicks_dict_ref
elif picktype == 'manual':
picks_dict = self.autopicks_dict
else:
return
if picks_dict:
return self.subtract_picks(picks_dict)
else:
if self.load_reference_picks():
return self.picks_rel_ref_corrected
def subtract_picks(self, picks_dict):
current_picks = self.current_picks_dict()
subtracted_picks = {}
for station, ps_dict in current_picks.items():
for phase, pick in ps_dict.items():
pick_ref_ps = picks_dict.get(station)
if not pick_ref_ps:
continue
pick_ref = pick_ref_ps.get(phase)
if not pick_ref:
continue
mpp = pick['mpp'] - pick_ref['mpp']
nwst_id = f'{pick["network"]}.{station}'
subtracted_picks[nwst_id] = mpp
return subtracted_picks
def _warn(self, message): def _warn(self, message):
self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Warning, 'Warning', message) self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Warning, 'Warning', message)
self.qmb.show() self.qmb.show()
-6
View File
@@ -27,10 +27,6 @@ class Metadata(object):
# saves which metadata files are from obspy dmt # saves which metadata files are from obspy dmt
self.obspy_dmt_invs = [] self.obspy_dmt_invs = []
if inventory: if inventory:
# make sure that no accidental backslashes mess up the path
if isinstance(inventory, str):
inventory = inventory.replace('\\', '/')
inventory = os.path.abspath(inventory)
if os.path.isdir(inventory): if os.path.isdir(inventory):
self.add_inventory(inventory) self.add_inventory(inventory)
if os.path.isfile(inventory): if os.path.isfile(inventory):
@@ -59,8 +55,6 @@ class Metadata(object):
:type path_to_inventory: str :type path_to_inventory: str
:return: None :return: None
""" """
path_to_inventory = path_to_inventory.replace('\\', '/')
path_to_inventory = os.path.abspath(path_to_inventory)
assert (os.path.isdir(path_to_inventory)), '{} is no directory'.format(path_to_inventory) assert (os.path.isdir(path_to_inventory)), '{} is no directory'.format(path_to_inventory)
if path_to_inventory not in self.inventories: if path_to_inventory not in self.inventories:
self.inventories.append(path_to_inventory) self.inventories.append(path_to_inventory)
+3
View File
@@ -22,11 +22,14 @@ class Event(ObsPyEvent):
:param path: path to event directory :param path: path to event directory
:type path: str :type path: str
""" """
# TODO: remove rootpath and database
self.pylot_id = path.split('/')[-1] self.pylot_id = path.split('/')[-1]
# initialize super class # initialize super class
super(Event, self).__init__(resource_id=ResourceIdentifier('smi:local/' + self.pylot_id)) super(Event, self).__init__(resource_id=ResourceIdentifier('smi:local/' + self.pylot_id))
self.path = path self.path = path
self.database = path.split('/')[-2]
self.datapath = os.path.split(path)[0] # path.split('/')[-3] self.datapath = os.path.split(path)[0] # path.split('/')[-3]
self.rootpath = '/' + os.path.join(*path.split('/')[:-3])
self.pylot_autopicks = {} self.pylot_autopicks = {}
self.pylot_picks = {} self.pylot_picks = {}
self.notes = '' self.notes = ''
+8 -10
View File
@@ -1,7 +1,6 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import os import os
from functools import lru_cache
try: try:
import pyqtgraph as pg import pyqtgraph as pg
@@ -26,14 +25,14 @@ def pick_linestyle_pg(picktype, key):
:return: Qt line style parameters :return: Qt line style parameters
:rtype: :rtype:
""" """
linestyles_manu = {'mpp': (QtCore.Qt.SolidLine, 2), linestyles_manu = {'mpp': (QtCore.Qt.SolidLine, 2.),
'epp': (QtCore.Qt.DashLine, 1), 'epp': (QtCore.Qt.DashLine, 1.),
'lpp': (QtCore.Qt.DashLine, 1), 'lpp': (QtCore.Qt.DashLine, 1.),
'spe': (QtCore.Qt.DashLine, 1)} 'spe': (QtCore.Qt.DashLine, 1.)}
linestyles_auto = {'mpp': (QtCore.Qt.DotLine, 2), linestyles_auto = {'mpp': (QtCore.Qt.DotLine, 2.),
'epp': (QtCore.Qt.DashDotLine, 1), 'epp': (QtCore.Qt.DashDotLine, 1.),
'lpp': (QtCore.Qt.DashDotLine, 1), 'lpp': (QtCore.Qt.DashDotLine, 1.),
'spe': (QtCore.Qt.DashDotLine, 1)} 'spe': (QtCore.Qt.DashDotLine, 1.)}
linestyles = {'manual': linestyles_manu, linestyles = {'manual': linestyles_manu,
'auto': linestyles_auto} 'auto': linestyles_auto}
return linestyles[picktype][key] return linestyles[picktype][key]
@@ -81,7 +80,6 @@ def which(program, parameter):
return None return None
@lru_cache(maxsize=128)
def make_pen(picktype, phase, key, quality): def make_pen(picktype, phase, key, quality):
""" """
Make PyQtGraph.QPen Make PyQtGraph.QPen
+2 -2
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}
+5 -4
View File
@@ -358,6 +358,8 @@ def get_bool(value):
False False
>>> get_bool(None) >>> get_bool(None)
None None
>>> get_bool('Stream')
'Stream'
""" """
if type(value) is bool: if type(value) is bool:
return value return value
@@ -1076,7 +1078,7 @@ def check4rotated(data, metadata=None, verbosity=1):
return wfs_in return wfs_in
# check metadata quality # check metadata quality
t_start = full_range(wfs_in)[0] t_start = full_range(wfs_in)
try: try:
azimuths = [] azimuths = []
dips = [] dips = []
@@ -1084,9 +1086,8 @@ def check4rotated(data, metadata=None, verbosity=1):
azimuths.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 err: except (KeyError, TypeError) as err:
logging.warning(f"Rotating not possible, not all azimuth and dip information " logging.error(f"{type(err)=} occurred: {err=} Rotating not possible, not all azimuth and dip information "
f"available in metadata. Stream remains unchanged.") f"available in metadata. Stream remains unchanged.")
logging.debug(f"Rotating not possible, {err=}, {type(err)=}")
return wfs_in return wfs_in
except Exception as err: except Exception as err:
print(f"Unexpected {err=}, {type(err)=}") print(f"Unexpected {err=}, {type(err)=}")
+12 -14
View File
@@ -8,7 +8,6 @@ import copy
import datetime import datetime
import getpass import getpass
import glob import glob
import logging
import multiprocessing import multiprocessing
import os import os
import subprocess import subprocess
@@ -1871,14 +1870,13 @@ class PickDlg(QDialog):
def __init__(self, parent=None, data=None, data_compare=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, show_comp_data=False, autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None, show_comp_data=False,
event=None, filteroptions=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
self.setAttribute(Qt.WA_DeleteOnClose) self.setAttribute(Qt.WA_DeleteOnClose)
# initialize attributes # initialize attributes
self.parameter = parameter self.parameter = parameter
model = self.parameter.get('taup_model')
self._embedded = embedded self._embedded = embedded
self.showCompData = show_comp_data self.showCompData = show_comp_data
self.station = station self.station = station
@@ -2271,8 +2269,8 @@ class PickDlg(QDialog):
arrivals = func[plot](source_origin.depth, arrivals = func[plot](source_origin.depth,
source_origin.latitude, source_origin.latitude,
source_origin.longitude, source_origin.longitude,
station_coords.get('latitude'), station_coords['latitude'],
station_coords.get('longitude'), station_coords['longitude'],
phases) phases)
self.arrivals = arrivals self.arrivals = arrivals
@@ -3837,7 +3835,7 @@ class TuneAutopicker(QWidget):
self.stb_names = ['aicARHfig', 'refSpick', 'el_S1pick', 'el_S2pick'] self.stb_names = ['aicARHfig', 'refSpick', 'el_S1pick', 'el_S2pick']
def add_parameters(self): def add_parameters(self):
self.paraBox = PylotParameterWidget(self.parameter, parent=self, windowflag=Qt.Widget) self.paraBox = PylotParaBox(self.parameter, parent=self, windowflag=Qt.Widget)
self.paraBox.set_tune_mode(True) self.paraBox.set_tune_mode(True)
self.update_eventID() self.update_eventID()
self.parameter_layout.addWidget(self.paraBox) self.parameter_layout.addWidget(self.paraBox)
@@ -4204,7 +4202,7 @@ class TuneAutopicker(QWidget):
self.qmb.show() self.qmb.show()
class PylotParameterWidget(QtWidgets.QWidget): class PylotParaBox(QtWidgets.QWidget):
accepted = QtCore.Signal(str) accepted = QtCore.Signal(str)
rejected = QtCore.Signal(str) rejected = QtCore.Signal(str)
@@ -4318,11 +4316,6 @@ class PylotParameterWidget(QtWidgets.QWidget):
grid = QtWidgets.QGridLayout() grid = QtWidgets.QGridLayout()
for index1, name in enumerate(parameter_names): for index1, name in enumerate(parameter_names):
if name in ['rootpath', 'database']:
logging.warning(
f'Deprecated parameter loaded: {name}. Check if datapath is still correct in parameter widget.'
)
continue
default_item = self.parameter.get_defaults()[name] default_item = self.parameter.get_defaults()[name]
tooltip = default_item['tooltip'] tooltip = default_item['tooltip']
tooltip += ' | type: {}'.format(default_item['type']) tooltip += ' | type: {}'.format(default_item['type'])
@@ -4892,7 +4885,7 @@ class PropTab(QWidget):
def getValues(self): def getValues(self):
return None return None
def resetValues(self, infile): def resetValues(self, infile=None):
return None return None
@@ -4989,7 +4982,12 @@ class InputsTab(PropTab):
else: else:
index = 2 index = 2
datapath = para.get('datapath') if not para.get('datapath') is None else '' datapath = para.get('datapath') if not para.get('datapath') is None else ''
values = {"data/dataRoot": self.dataDirEdit.setText("%s" % datapath), rootpath = para.get('rootpath') if not para.get('rootpath') is None else ''
database = para.get('database') if not para.get('database') is None else ''
if isinstance(database, int):
database = str(database)
path = os.path.join(os.path.expanduser('~'), rootpath, datapath, database)
values = {"data/dataRoot": self.dataDirEdit.setText("%s" % path),
"user/FullName": self.fullNameEdit.text(), "user/FullName": self.fullNameEdit.text(),
"data/Structure": self.structureSelect.setCurrentIndex(index), "data/Structure": self.structureSelect.setCurrentIndex(index),
"tstart": self.tstartBox.setValue(0), "tstart": self.tstartBox.setValue(0),
-2
View File
@@ -1,2 +0,0 @@
# -*- coding: utf-8 -*-
#
@@ -1,90 +0,0 @@
############################# correlation parameters #####################################
# min_corr_stacking: minimum correlation coefficient for building beam trace
# min_corr_export: minimum correlation coefficient for pick export
# min_stack: minimum number of stations for building beam trace
# t_before: correlation window before pick
# t_after: correlation window after pick#
# cc_maxlag: maximum shift for initial correlation
# cc_maxlag2: maximum shift for second (final) correlation (also for calculating pick uncertainty)
# initial_pick_outlier_threshold: (hopefully) threshold for excluding large outliers of initial (AIC) picks
# export_threshold: automatically exclude all onsets which deviate more than this threshold from corrected taup onsets
# min_picks_export: minimum number of correlated picks for export
# min_picks_autopylot: minimum number of reference autopicks picks to continue with event
# check_RMS: do RMS check to search for restitution errors (very experimental)
# use_taupy_onsets: use taupy onsets as reference picks instead of external picks
# station_list: use the following stations as reference for stacking
# use_stacked_trace: use existing stacked trace if found (spare re-computation)
# data_dir: obspyDMT data subdirectory (e.g. 'raw', 'processed')
# pickfile_extension: use quakeML files (PyLoT output) with the following extension, e.g. '_autopylot' for pickfiles
# such as 'PyLoT_20170501_141822_autopylot.xml'
logging: info
pick_phases: ['P', 'S']
# P-phase
P:
min_corr_stacking: 0.8
min_corr_export: 0.6
min_stack: 20
t_before: 30.
t_after: 50.
cc_maxlag: 50.
cc_maxlag2: 5.
initial_pick_outlier_threshold: 30.
export_threshold: 2.5
min_picks_export: 100
min_picks_autopylot: 50
check_RMS: True
use_taupy_onsets: False
station_list: ['HU.MORH', 'HU.TIH', 'OX.FUSE', 'OX.BAD']
use_stacked_trace: False
data_dir: 'processed'
pickfile_extension: '_autopylot'
dt_stacking: [250, 250]
# filter for first correlation (rough)
filter_options:
freqmax: 0.5
freqmin: 0.03
# filter for second correlation (fine)
filter_options_final:
freqmax: 0.5
freqmin: 0.03
filter_type: bandpass
sampfreq: 20.0
# S-phase
S:
min_corr_stacking: 0.7
min_corr_export: 0.6
min_stack: 20
t_before: 60.
t_after: 60.
cc_maxlag: 100.
cc_maxlag2: 25.
initial_pick_outlier_threshold: 30.
export_threshold: 5.0
min_picks_export: 200
min_picks_autopylot: 50
check_RMS: True
use_taupy_onsets: False
station_list: ['HU.MORH','HU.TIH', 'OX.FUSE', 'OX.BAD']
use_stacked_trace: False
data_dir: 'processed'
pickfile_extension: '_autopylot'
dt_stacking: [250, 250]
# filter for first correlation (rough)
filter_options:
freqmax: 0.1
freqmin: 0.01
# filter for second correlation (fine)
filter_options_final:
freqmax: 0.2
freqmin: 0.01
filter_type: bandpass
sampfreq: 20.0
File diff suppressed because it is too large Load Diff
@@ -1,40 +0,0 @@
#!/bin/bash
#ulimit -s 8192
#ulimit -v $(ulimit -v | awk '{printf("%d",$1*0.95)}')
#ulimit -v
#655360
source /opt/anaconda3/etc/profile.d/conda.sh
conda activate pylot_311
NSLOTS=20
#qsub -l low -cwd -l "os=*stretch" -pe smp 40 submit_pick_corr_correction.sh
#$ -l low
#$ -l h_vmem=6G
#$ -cwd
#$ -pe smp 20
#$ -N corr_pick
export PYTHONPATH="$PYTHONPATH:/home/marcel/git/pylot_tools/"
export PYTHONPATH="$PYTHONPATH:/home/marcel/git/"
export PYTHONPATH="$PYTHONPATH:/home/marcel/git/pylot/"
#export MKL_NUM_THREADS=${NSLOTS:=1}
#export NUMEXPR_NUM_THREADS=${NSLOTS:=1}
#export OMP_NUM_THREADS=${NSLOTS:=1}
#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 0 -istop 100
#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 100 -istop 200
#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M6.0-6.5' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 0 -istop 100
#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 100 -istop 200
#python pick_correlation_correction.py 'H:\sciebo\dmt_database' 'H:\Sciebo\dmt_database\pylot_alparray_mantle_corr_S_0.01-0.2.in' -pd -n 4 -t
pylot_infile='/home/marcel/.pylot/pylot_alparray_syn_fwi_mk6_it3.in'
#pylot_infile='/home/marcel/.pylot/pylot_adriaarray_corr_P_and_S.in'
# THIS SCRIPT SHOLD BE CALLED BY "submit_to_grid_engine.py" using the following line:
python pick_correlation_correction.py $1 $pylot_infile -pd -n ${NSLOTS:=1} -istart $2 --params 'parameters_fwi_mk6_it3.yaml'
#--event_blacklist eventlist.txt
@@ -1,23 +0,0 @@
#!/usr/bin/env python
import subprocess
fnames = [
('/data/AlpArray_Data/dmt_database_synth_model_mk6_it3_no_rotation', 0),
]
#fnames = [('/data/AlpArray_Data/dmt_database_mantle_0.01-0.2_SKS-phase', 0),
# ('/data/AlpArray_Data/dmt_database_mantle_0.01-0.2_S-phase', 0),]
####
script_location = '/home/marcel/VersionCtrl/git/code_base/correlation_picker/submit_pick_corr_correction.sh'
####
for fnin, istart in fnames:
input_cmds = f'qsub -q low.q@minos15,low.q@minos14,low.q@minos13,low.q@minos12,low.q@minos11 {script_location} {fnin} {istart}'
print(input_cmds)
print(subprocess.check_output(input_cmds.split()))
-61
View File
@@ -1,61 +0,0 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import glob
import json
from obspy import read_events
from pylot.core.util.dataprocessing import Metadata
from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
def get_event_obspy_dmt(eventdir):
event_pkl_file = os.path.join(eventdir, 'info', 'event.pkl')
if not os.path.exists(event_pkl_file):
raise IOError('Could not find event path for event: {}'.format(eventdir))
event = qml_from_obspyDMT(event_pkl_file)
return event
def get_event_pylot(eventdir, extension=''):
event_id = get_event_id(eventdir)
filename = os.path.join(eventdir, 'PyLoT_{}{}.xml'.format(event_id, extension))
if not os.path.isfile(filename):
return
cat = read_events(filename)
return cat[0]
def get_event_id(eventdir):
event_id = os.path.split(eventdir)[-1]
return event_id
def get_picks(eventdir, extension=''):
event_id = get_event_id(eventdir)
filename = 'PyLoT_{}{}.xml'
filename = filename.format(event_id, extension)
fpath = os.path.join(eventdir, filename)
fpaths = glob.glob(fpath)
if len(fpaths) == 1:
cat = read_events(fpaths[0])
picks = cat[0].picks
return picks
elif len(fpaths) == 0:
print('get_picks: File not found: {}'.format(fpath))
return
print(f'WARNING: Ambiguous pick file specification. Found the following pick files {fpaths}\nFilemask: {fpath}')
return
def write_json(object, fname):
with open(fname, 'w') as outfile:
json.dump(object, outfile, sort_keys=True, indent=4)
def get_metadata(eventdir):
metadata_path = os.path.join(eventdir, 'resp')
metadata = Metadata(inventory=metadata_path, verbosity=0)
return metadata
+12 -7
View File
@@ -1,7 +1,12 @@
Cartopy==0.23.0 # This file may be used to create an environment using:
joblib==1.4.2 # $ conda create --name <env> --file <this file>
obspy==1.4.1 # platform: win-64
pyaml==24.7.0 cartopy=0.20.2
pyqtgraph==0.13.7 matplotlib-base=3.3.4
PySide2==5.15.8 numpy=1.22.3
pytest==8.3.2 obspy=1.3.0
pyqtgraph=0.12.4
pyside2=5.13.2
python=3.8.12
qt=5.12.9
scipy=1.8.0
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
@@ -1,99 +0,0 @@
%This is a parameter input file for PyLoT/autoPyLoT.
%All main and special settings regarding data handling
%and picking are to be set here!
%Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings#
dmt_database_test #datapath# %data path
20171010_063224.a #eventID# %event ID for single event processing (* for all events found in database)
#invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#NLLoc settings#
None #nllocbin# %path to NLLoc executable
None #nllocroot# %root of NLLoc-processing directory
None #phasefile# %name of autoPyLoT-output phase file for NLLoc
None #ctrfile# %name of autoPyLoT-output control file for NLLoc
ttime #ttpatter# %pattern of NLLoc ttimes from grid
AUTOLOC_nlloc #outpatter# %pattern of NLLoc-output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#parameters for seismic moment estimation#
3530.0 #vp# %average P-wave velocity
2500.0 #rho# %average rock density [kg/m^3]
300.0 0.8 #Qp# %quality factor for P waves (Qp*f^a); list(Qp, a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#settings local magnitude#
1.0 1.0 1.0 #WAscaling# %Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm] If zeros are set, original Richter magnitude is calculated!
1.0 1.0 #magscaling# %Scaling relation for derived local magnitude [a*Ml+b]. If zeros are set, no scaling of network magnitude is applied!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#filter settings#
0.03 0.03 #minfreq# %Lower filter frequency [P, S]
0.5 0.5 #maxfreq# %Upper filter frequency [P, S]
4 4 #filter_order# %filter order [P, S]
bandpass bandpass #filter_type# %filter type (bandpass, bandstop, lowpass, highpass) [P, S]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#common settings picker#
global #extent# %extent of array ("local", "regional" or "global")
-100.0 #pstart# %start time [s] for calculating CF for P-picking (if TauPy: seconds relative to estimated onset)
50.0 #pstop# %end time [s] for calculating CF for P-picking (if TauPy: seconds relative to estimated onset)
-50.0 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking
50.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
ak135 #taup_model# %Define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
P,Pdiff,S,SKS #taup_phases# %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
0.03 0.5 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
0.01 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
0.03 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
0.01 0.5 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz]
#special settings for calculating CF#
%!!Edit the following only if you know what you are doing!!%
#Z-component#
HOS #algoP# %choose algorithm for P-onset determination (HOS, ARZ, or AR3)
300.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s]
4 #hosorder# %for HOS-picker, order of Higher Order Statistics
2 #Parorder# %for AR-picker, order of AR process of Z-component
16.0 #tdet1z# %for AR-picker, length of AR determination window [s] for Z-component, 1st pick
10.0 #tpred1z# %for AR-picker, length of AR prediction window [s] for Z-component, 1st pick
12.0 #tdet2z# %for AR-picker, length of AR determination window [s] for Z-component, 2nd pick
6.0 #tpred2z# %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick
0.001 #addnoise# %add noise to seismogram for stable AR prediction
60.0 5.0 20.0 12.0 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]
50.0 #pickwinP# %for initial AIC pick, length of P-pick window [s]
30.0 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)
2.0 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s]
2.0 #tsmoothP# %for HOS/AR, take average of samples in this time window for smoothing CF [s]
0.006 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P)
2.0 #nfacP# %for HOS/AR, noise factor for noise level determination (P)
#H-components#
ARH #algoS# %choose algorithm for S-onset determination (ARH or AR3)
12.0 #tdet1h# %for HOS/AR, length of AR-determination window [s], H-components, 1st pick
6.0 #tpred1h# %for HOS/AR, length of AR-prediction window [s], H-components, 1st pick
8.0 #tdet2h# %for HOS/AR, length of AR-determinaton window [s], H-components, 2nd pick
4.0 #tpred2h# %for HOS/AR, length of AR-prediction window [s], H-components, 2nd pick
4 #Sarorder# %for AR-picker, order of AR process of H-components
100.0 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
195.0 #pickwinS# %for initial AIC pick, length of S-pick window [s]
60.0 10.0 30.0 12.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]
22.0 #aictsmoothS# %for AIC-picker, take average of samples in this time window for smoothing of AIC-function [s]
20.0 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S)
0.001 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S)
2.0 #nfacS# %for AR-picker, noise factor for noise level determination (S)
#first-motion picker#
1 #minfmweight# %minimum required P weight for first-motion determination
3.0 #minFMSNR# %miniumum required SNR for first-motion determination
10.0 #fmpickwin# %pick window [s] around P onset for calculating zero crossings
#quality assessment#
0.1 0.2 0.4 0.8 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P
4.0 8.0 16.0 32.0 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S
0.005 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected
1.1 #minAICPSNR# %below this SNR the initial P pick is rejected
0.002 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected
1.3 #minAICSSNR# %below this SNR the initial S pick is rejected
20.0 #minsiglength# %length of signal part for which amplitudes must exceed noiselevel [s]
1.0 #noisefactor# %noiselevel*noisefactor=threshold
10.0 #minpercent# %required percentage of amplitudes exceeding threshold
0.1 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude
100.0 #mdttolerance# %maximum allowed deviation of P picks from median [s]
50.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram
25.0 #jackfactor# %pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor
-67
View File
@@ -1,67 +0,0 @@
import os
import pytest
from obspy import read_events
from autoPyLoT import autoPyLoT
class TestAutopickerGlobal():
def init(self):
self.params_infile = 'pylot_alparray_mantle_corr_stack_0.03-0.5.in'
self.test_event_dir = 'dmt_database_test'
self.fname_outfile_xml = os.path.join(
self.test_event_dir, '20171010_063224.a', 'PyLoT_20171010_063224.a_autopylot.xml'
)
# check if the input files exist
if not os.path.isfile(self.params_infile):
print(f'Test input file {os.path.abspath(self.params_infile)} not found.')
return False
if not os.path.exists(self.test_event_dir):
print(
f'Test event directory not found at location "{os.path.abspath(self.test_event_dir)}". '
f'Make sure to load it from the website first.'
)
return False
return True
def test_autopicker(self):
assert self.init(), 'Initialization failed due to missing input files.'
# check for output file in test directory and remove it if necessary
if os.path.isfile(self.fname_outfile_xml):
os.remove(self.fname_outfile_xml)
autoPyLoT(inputfile=self.params_infile, eventid='20171010_063224.a', obspyDMT_wfpath='processed')
# test for different known output files if they are identical or not
compare_pickfiles(self.fname_outfile_xml, 'PyLoT_20171010_063224.a_autopylot.xml', True)
compare_pickfiles(self.fname_outfile_xml, 'PyLoT_20171010_063224.a_saved_from_GUI.xml', True)
compare_pickfiles(self.fname_outfile_xml, 'PyLoT_20171010_063224.a_corrected_taup_times_0.03-0.5_P.xml', False)
def compare_pickfiles(pickfile1: str, pickfile2: str, samefile: bool = True) -> None:
"""
Compare the pick times and errors from two pick files.
Parameters:
pickfile1 (str): The path to the first pick file.
pickfile2 (str): The path to the second pick file.
samefile (bool): A flag indicating whether the two files are expected to be the same. Defaults to True.
Returns:
None
"""
cat1 = read_events(pickfile1)
cat2 = read_events(pickfile2)
picks1 = sorted(cat1[0].picks, key=lambda pick: str(pick.waveform_id))
picks2 = sorted(cat2[0].picks, key=lambda pick: str(pick.waveform_id))
pick_times1 = [pick.time for pick in picks1]
pick_times2 = [pick.time for pick in picks2]
pick_terrs1 = [pick.time_errors for pick in picks1]
pick_terrs2 = [pick.time_errors for pick in picks2]
# check if times and errors are identical or not depending on the samefile flag
assert (pick_times1 == pick_times2) is samefile, 'Pick times error'
assert (pick_terrs1 == pick_terrs2) is samefile, 'Pick time errors errors'
@@ -4,8 +4,10 @@
%Parameters are optimized for %extent data sets! %Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings# #main settings#
/home/darius/alparray/waveforms_used #datapath# %data path /home/darius #rootpath# %project path
e0093.173.16 #eventID# %event ID for single event processing (* for all events found in datapath) alparray #datapath# %data path
waveforms_used #database# %name of data base
e0093.173.16 #eventID# %event ID for single event processing (* for all events found in database)
/home/darius/alparray/metadata #invdir# %full path to inventory or dataless-seed file /home/darius/alparray/metadata #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output True #apverbose# %choose 'True' or 'False' for terminal output
@@ -41,7 +43,6 @@ global #extent# %extent of a
875.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking 875.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
False #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF False #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
IASP91 #taup_model# %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6 IASP91 #taup_model# %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
P,Pdiff,S,Sdiff #taup_phases# %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
0.01 0.1 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 0.01 0.1 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
0.001 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] 0.001 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
0.01 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] 0.01 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
@@ -4,8 +4,10 @@
%Parameters are optimized for %extent data sets! %Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings# #main settings#
/home/darius/alparray/waveforms_used #datapath# %data path /home/darius #rootpath# %project path
e0093.173.16 #eventID# %event ID for single event processing (* for all events found in datapath) alparray #datapath# %data path
waveforms_used #database# %name of data base
e0093.173.16 #eventID# %event ID for single event processing (* for all events found in database)
/home/darius/alparray/metadata #invdir# %full path to inventory or dataless-seed file /home/darius/alparray/metadata #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output True #apverbose# %choose 'True' or 'False' for terminal output
@@ -41,7 +43,6 @@ global #extent# %extent of a
875.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking 875.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
IASP91 #taup_model# %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6 IASP91 #taup_model# %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
P,Pdiff,S,Sdiff #taup_phases# %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
0.01 0.1 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 0.01 0.1 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
0.001 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] 0.001 0.5 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
0.01 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] 0.01 0.5 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
@@ -1,7 +1,6 @@
import os import os
import sys import sys
import unittest import unittest
import pytest
import obspy import obspy
from obspy import UTCDateTime from obspy import UTCDateTime
@@ -25,7 +24,7 @@ class HidePrints:
if self.hide: if self.hide:
self._original_stdout = sys.stdout self._original_stdout = sys.stdout
devnull = open(os.devnull, "w") devnull = open(os.devnull, "w")
#sys.stdout = devnull sys.stdout = devnull
def __exit__(self, exc_type, exc_val, exc_tb): def __exit__(self, exc_type, exc_val, exc_tb):
if self.hide: if self.hide:
@@ -106,6 +105,7 @@ class TestAutopickStation(unittest.TestCase):
# show complete diff when difference in results dictionaries are found # show complete diff when difference in results dictionaries are found
self.maxDiff = None self.maxDiff = None
# @skip("Works")
def test_autopickstation_taupy_disabled_gra1(self): def test_autopickstation_taupy_disabled_gra1(self):
expected = { expected = {
'P': {'picker': 'auto', 'snrdb': 15.405649120980094, 'weight': 0, 'Mo': None, 'marked': [], 'Mw': None, 'P': {'picker': 'auto', 'snrdb': 15.405649120980094, 'weight': 0, 'Mo': None, 'marked': [], 'Mw': None,
@@ -121,8 +121,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.gra1, pickparam=self.pickparam_taupy_disabled, result, station = autopickstation(wfstream=self.gra1, pickparam=self.pickparam_taupy_disabled,
metadata=(None, None)) metadata=(None, None))
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('GRA1', station) self.assertEqual('GRA1', station)
def test_autopickstation_taupy_enabled_gra1(self): def test_autopickstation_taupy_enabled_gra1(self):
@@ -140,8 +140,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.gra1, pickparam=self.pickparam_taupy_enabled, result, station = autopickstation(wfstream=self.gra1, pickparam=self.pickparam_taupy_enabled,
metadata=self.metadata, origin=self.origin) metadata=self.metadata, origin=self.origin)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('GRA1', station) self.assertEqual('GRA1', station)
def test_autopickstation_taupy_disabled_gra2(self): def test_autopickstation_taupy_disabled_gra2(self):
@@ -157,8 +157,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.gra2, pickparam=self.pickparam_taupy_disabled, result, station = autopickstation(wfstream=self.gra2, pickparam=self.pickparam_taupy_disabled,
metadata=(None, None)) metadata=(None, None))
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('GRA2', station) self.assertEqual('GRA2', station)
def test_autopickstation_taupy_enabled_gra2(self): def test_autopickstation_taupy_enabled_gra2(self):
@@ -175,8 +175,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.gra2, pickparam=self.pickparam_taupy_enabled, result, station = autopickstation(wfstream=self.gra2, pickparam=self.pickparam_taupy_enabled,
metadata=self.metadata, origin=self.origin) metadata=self.metadata, origin=self.origin)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('GRA2', station) self.assertEqual('GRA2', station)
def test_autopickstation_taupy_disabled_ech(self): def test_autopickstation_taupy_disabled_ech(self):
@@ -190,8 +190,8 @@ class TestAutopickStation(unittest.TestCase):
'fm': None, 'spe': None, 'channel': u'LHE'}} 'fm': None, 'spe': None, 'channel': u'LHE'}}
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.ech, pickparam=self.pickparam_taupy_disabled) result, station = autopickstation(wfstream=self.ech, pickparam=self.pickparam_taupy_disabled)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('ECH', station) self.assertEqual('ECH', station)
def test_autopickstation_taupy_enabled_ech(self): def test_autopickstation_taupy_enabled_ech(self):
@@ -208,8 +208,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.ech, pickparam=self.pickparam_taupy_enabled, result, station = autopickstation(wfstream=self.ech, pickparam=self.pickparam_taupy_enabled,
metadata=self.metadata, origin=self.origin) metadata=self.metadata, origin=self.origin)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('ECH', station) self.assertEqual('ECH', station)
def test_autopickstation_taupy_disabled_fiesa(self): def test_autopickstation_taupy_disabled_fiesa(self):
@@ -224,8 +224,8 @@ class TestAutopickStation(unittest.TestCase):
'fm': None, 'spe': None, 'channel': u'LHE'}} 'fm': None, 'spe': None, 'channel': u'LHE'}}
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.fiesa, pickparam=self.pickparam_taupy_disabled) result, station = autopickstation(wfstream=self.fiesa, pickparam=self.pickparam_taupy_disabled)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('FIESA', station) self.assertEqual('FIESA', station)
def test_autopickstation_taupy_enabled_fiesa(self): def test_autopickstation_taupy_enabled_fiesa(self):
@@ -242,8 +242,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.fiesa, pickparam=self.pickparam_taupy_enabled, result, station = autopickstation(wfstream=self.fiesa, pickparam=self.pickparam_taupy_enabled,
metadata=self.metadata, origin=self.origin) metadata=self.metadata, origin=self.origin)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
compare_dicts(expected=expected['S'], result=result['S'], hint='S-') self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
self.assertEqual('FIESA', station) self.assertEqual('FIESA', station)
def test_autopickstation_gra1_z_comp_missing(self): def test_autopickstation_gra1_z_comp_missing(self):
@@ -271,9 +271,8 @@ class TestAutopickStation(unittest.TestCase):
'fm': 'N', 'channel': None}} 'fm': 'N', 'channel': None}}
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=wfstream, pickparam=self.pickparam_taupy_disabled, result, station = autopickstation(wfstream=wfstream, pickparam=self.pickparam_taupy_disabled,
metadata=(None, None), iplot=2) metadata=(None, None))
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertEqual(expected, result)
compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
self.assertEqual('GRA1', station) self.assertEqual('GRA1', station)
def test_autopickstation_a106_taupy_enabled(self): def test_autopickstation_a106_taupy_enabled(self):
@@ -291,9 +290,7 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.a106, pickparam=self.pickparam_taupy_enabled, result, station = autopickstation(wfstream=self.a106, pickparam=self.pickparam_taupy_enabled,
metadata=self.metadata, origin=self.origin) metadata=self.metadata, origin=self.origin)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertEqual(expected, result)
compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
def test_autopickstation_station_missing_in_metadata(self): def test_autopickstation_station_missing_in_metadata(self):
"""This station is not in the metadata, but Taupy is enabled. Taupy should exit cleanly and modify the starttime """This station is not in the metadata, but Taupy is enabled. Taupy should exit cleanly and modify the starttime
@@ -314,37 +311,8 @@ class TestAutopickStation(unittest.TestCase):
with HidePrints(): with HidePrints():
result, station = autopickstation(wfstream=self.a005a, pickparam=self.pickparam_taupy_enabled, result, station = autopickstation(wfstream=self.a005a, pickparam=self.pickparam_taupy_enabled,
metadata=self.metadata, origin=self.origin) metadata=self.metadata, origin=self.origin)
compare_dicts(expected=expected['P'], result=result['P'], hint='P-') self.assertEqual(expected, result)
compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
def run_dict_comparison(result, expected):
for key, expected_value in expected.items():
if isinstance(expected_value, dict):
run_dict_comparison(result[key], expected[key])
else:
res = result[key]
if isinstance(res, UTCDateTime) and isinstance(expected_value, UTCDateTime):
res = res.timestamp
expected_value = expected_value.timestamp
assert expected_value == pytest.approx(res), f'{key}: {expected_value} != {res}'
def compare_dicts(result, expected, hint=''):
try:
run_dict_comparison(result, expected)
except AssertionError:
raise AssertionError(f'{hint}Dictionaries not equal.'
f'\n\n<<Expected>>\n{pretty_print_dict(expected)}'
f'\n\n<<Result>>\n{pretty_print_dict(result)}')
def pretty_print_dict(dct):
retstr = ''
for key, value in sorted(dct.items(), key=lambda x: x[0]):
retstr += f"{key} : {value}\n"
return retstr
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
@@ -1,76 +0,0 @@
import pytest
from obspy import read, Trace, UTCDateTime
from pylot.correlation.pick_correlation_correction import XCorrPickCorrection
class TestXCorrPickCorrection():
def setup(self):
self.make_test_traces()
self.make_test_picks()
self.t_before = 2.
self.t_after = 2.
self.cc_maxlag = 0.5
def make_test_traces(self):
# take first trace of test Stream from obspy
tr1 = read()[0]
# filter trace
tr1.filter('bandpass', freqmin=1, freqmax=20)
# make a copy and shift the copy by 0.1 s
tr2 = tr1.copy()
tr2.stats.starttime += 0.1
self.trace1 = tr1
self.trace2 = tr2
def make_test_picks(self):
# create an artificial reference pick on reference trace (trace1) and another one on the 0.1 s shifted trace
self.tpick1 = UTCDateTime('2009-08-24T00:20:07.7')
# shift the second pick by 0.2 s, the correction should be around 0.1 s now
self.tpick2 = self.tpick1 + 0.2
def test_slice_trace_okay(self):
self.setup()
xcpc = XCorrPickCorrection(UTCDateTime(), Trace(), UTCDateTime(), Trace(),
t_before=self.t_before, t_after=self.t_after, cc_maxlag=self.cc_maxlag)
test_trace = self.trace1
pick_time = self.tpick2
sliced_trace = xcpc.slice_trace(test_trace, pick_time)
assert ((sliced_trace.stats.starttime == pick_time - self.t_before - self.cc_maxlag / 2)
and (sliced_trace.stats.endtime == pick_time + self.t_after + self.cc_maxlag / 2))
def test_slice_trace_fails(self):
self.setup()
test_trace = self.trace1
pick_time = self.tpick1
with pytest.raises(ValueError):
xcpc = XCorrPickCorrection(UTCDateTime(), Trace(), UTCDateTime(), Trace(),
t_before=self.t_before + 20, t_after=self.t_after, cc_maxlag=self.cc_maxlag)
xcpc.slice_trace(test_trace, pick_time)
with pytest.raises(ValueError):
xcpc = XCorrPickCorrection(UTCDateTime(), Trace(), UTCDateTime(), Trace(),
t_before=self.t_before, t_after=self.t_after + 50, cc_maxlag=self.cc_maxlag)
xcpc.slice_trace(test_trace, pick_time)
def test_cross_correlation(self):
self.setup()
# create XCorrPickCorrection object
xcpc = XCorrPickCorrection(self.tpick1, self.trace1, self.tpick2, self.trace2, t_before=self.t_before,
t_after=self.t_after, cc_maxlag=self.cc_maxlag)
# execute correlation
correction, cc_max, uncert, fwfm = xcpc.cross_correlation(False, '', '')
# define awaited test result
test_result = (-0.09983091718314982, 0.9578431835689154, 0.0015285160561610929, 0.03625786256084631)
# check results
assert pytest.approx(test_result, rel=1e-6) == (correction, cc_max, uncert, fwfm)