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

Merge branch 'develop' into filterOptions

Browse Source
This commit is contained in:
Marcel Paffrath 2017-06-22 11:26:02 +02:00
commit 97ee03b443
17 changed files with 1258 additions and 430 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

668
QtPyLoT.py
View File

File diff suppressed because it is too large Load Diff

84
autoPyLoT.py
View File

@ -16,9 +16,9 @@ import pylot.core.loc.focmec as focmec
import pylot.core.loc.hash as hash import pylot.core.loc.hash as hash
import pylot.core.loc.nll as nll import pylot.core.loc.nll as nll
#from PySide.QtGui import QWidget, QInputDialog #from PySide.QtGui import QWidget, QInputDialog
from pylot.core.analysis.magnitude import MomentMagnitude, RichterMagnitude from pylot.core.analysis.magnitude import MomentMagnitude, LocalMagnitude
from pylot.core.io.data import Data from pylot.core.io.data import Data
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import PylotParameter
from pylot.core.pick.autopick import autopickevent, iteratepicker from pylot.core.pick.autopick import autopickevent, iteratepicker
from pylot.core.util.dataprocessing import restitute_data, read_metadata, \ from pylot.core.util.dataprocessing import restitute_data, read_metadata, \
remove_underscores remove_underscores
@ -35,7 +35,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
:param inputfile: path to the input file containing all parameter :param inputfile: path to the input file containing all parameter
information for automatic picking (for formatting details, see. information for automatic picking (for formatting details, see.
`~pylot.core.io.inputs.AutoPickParameter` `~pylot.core.io.inputs.PylotParameter`
:type inputfile: str :type inputfile: str
:return: :return:
@ -71,13 +71,13 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
if not parameter: if not parameter:
if inputfile: if inputfile:
parameter = AutoPickParameter(inputfile) parameter = PylotParameter(inputfile)
iplot = parameter['iplot'] iplot = parameter['iplot']
else: else:
print('No parameters set and no input file given. Choose either of both.') print('No parameters set and no input file given. Choose either of both.')
return return
else: else:
if not type(parameter) == AutoPickParameter: if not type(parameter) == PylotParameter:
print('Wrong input type for parameter: {}'.format(type(parameter))) print('Wrong input type for parameter: {}'.format(type(parameter)))
return return
if inputfile: if inputfile:
@ -130,21 +130,35 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
print("!!No source parameter estimation possible!!") print("!!No source parameter estimation possible!!")
print(" !!! ") print(" !!! ")
datapath = datastructure.expandDataPath() if not input_dict:
if fnames == 'None' and not parameter.hasParam('eventID'): # started in production mode
# multiple event processing datapath = datastructure.expandDataPath()
# read each event in database if fnames == 'None' and not parameter['eventID']:
events = [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)] # multiple event processing
elif fnames == 'None' and parameter.hasParam('eventID'): # read each event in database
# single event processing events = [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]
events = glob.glob(os.path.join(datapath, parameter.get('eventID'))) elif fnames == 'None' and parameter['eventID']:
# single event processing
events = glob.glob(os.path.join(datapath, parameter.get('eventID')))
else:
# autoPyLoT was initialized from GUI
events = []
events.append(eventid)
evID = os.path.split(eventid)[-1]
locflag = 2
else: else:
# autoPyLoT was initialized from GUI # started in tune mode
datapath = os.path.join(parameter['rootpath'],
parameter['datapath'])
events = [] events = []
events.append(eventid) events.append(os.path.join(datapath,
evID = os.path.split(eventid)[-1] parameter['database'],
locflag = 2 parameter['eventID']))
if not events:
print('autoPyLoT: No events given. Return!')
return
for event in events: for event in events:
if fnames == 'None': if fnames == 'None':
data.setWFData(glob.glob(os.path.join(datapath, event, '*'))) data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
@ -238,9 +252,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
for station, props in moment_mag.moment_props.items(): for station, props in moment_mag.moment_props.items():
picks[station]['P'].update(props) picks[station]['P'].update(props)
evt = moment_mag.updated_event() evt = moment_mag.updated_event()
local_mag = RichterMagnitude(corr_dat, evt, local_mag = LocalMagnitude(corr_dat, evt,
parameter.get('sstop'), True,\ parameter.get('sstop'), parameter.get('WAscaling'), \
iplot) True, iplot)
for station, amplitude in local_mag.amplitudes.items(): for station, amplitude in local_mag.amplitudes.items():
picks[station]['S']['Ao'] = amplitude.generic_amplitude picks[station]['S']['Ao'] = amplitude.generic_amplitude
evt = local_mag.updated_event() evt = local_mag.updated_event()
@ -296,9 +310,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
for station, props in moment_mag.moment_props.items(): for station, props in moment_mag.moment_props.items():
picks[station]['P'].update(props) picks[station]['P'].update(props)
evt = moment_mag.updated_event() evt = moment_mag.updated_event()
local_mag = RichterMagnitude(corr_dat, evt, local_mag = LocalMagnitude(corr_dat, evt,
parameter.get('sstop'), True, \ parameter.get('sstop'), parameter.get('WAscaling'), \
iplot) True, iplot)
for station, amplitude in local_mag.amplitudes.items(): for station, amplitude in local_mag.amplitudes.items():
picks[station]['S']['Ao'] = amplitude.generic_amplitude picks[station]['S']['Ao'] = amplitude.generic_amplitude
evt = local_mag.updated_event() evt = local_mag.updated_event()
@ -357,8 +371,15 @@ if __name__ == "__main__":
# parse arguments # parse arguments
parser = argparse.ArgumentParser( parser = argparse.ArgumentParser(
description='''autoPyLoT automatically picks phase onset times using higher order statistics, description='''autoPyLoT automatically picks phase onset times using higher order statistics,
autoregressive prediction and AIC''') autoregressive prediction and AIC followed by locating the seismic events using
NLLoc''')
#parser.add_argument('-d', '-D', '--input_dict', type=str,
# action='store',
# help='''optional, dictionary containing processing parameters''')
#parser.add_argument('-p', '-P', '--parameter', type=str,
# action='store',
# help='''parameter file, default=None''')
parser.add_argument('-i', '-I', '--inputfile', type=str, parser.add_argument('-i', '-I', '--inputfile', type=str,
action='store', action='store',
help='''full path to the file containing the input help='''full path to the file containing the input
@ -369,17 +390,14 @@ if __name__ == "__main__":
parser.add_argument('-e', '-E', '--eventid', type=str, parser.add_argument('-e', '-E', '--eventid', type=str,
action='store', action='store',
help='''optional, event path incl. event ID''') help='''optional, event path incl. event ID''')
# parser.add_argument('-p', '-P', '--plot', action='store',
# help='show interactive plots')
parser.add_argument('-s', '-S', '--spath', type=str, parser.add_argument('-s', '-S', '--spath', type=str,
action='store', action='store',
help='''optional, save path for autoPyLoT output''') help='''optional, save path for autoPyLoT output''')
parser.add_argument('-v', '-V', '--version', action='version', #parser.add_argument('-v', '-V', '--version', action='version',
version='autoPyLoT ' + __version__, # version='autoPyLoT ' + __version__,
help='show version information and exit') # help='show version information and exit')
cla = parser.parse_args() cla = parser.parse_args()
picks, mainFig = autoPyLoT(inputfile=str(cla.inputfile), picks = autoPyLoT(inputfile=str(cla.inputfile), fnames=str(cla.fnames),
fnames=str(cla.fnames), eventid=str(cla.eventid), eventid=str(cla.eventid), savepath=str(cla.spath))
savepath=str(cla.spath))

2
pylot/RELEASE-VERSION
View File

@ -1 +1 @@
62fa-dirty f91e1-dirty

40
pylot/core/analysis/magnitude.py
View File

@ -17,10 +17,17 @@ from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, \
from pylot.core.util.utils import common_range, fit_curve from pylot.core.util.utils import common_range, fit_curve
def richter_magnitude_scaling(delta): def richter_magnitude_scaling(delta):
relation = np.loadtxt(os.path.join(os.path.expanduser('~'), distance = np.array([0, 10, 20, 25, 30, 35,40, 45, 50, 60, 70, 75, 85, 90, 100, 110,
'.pylot', 'richter_scaling.data')) 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 230, 240, 250,
260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380,
390, 400, 430, 470, 510, 560, 600, 700, 800, 900, 1000])
richter_scaling = np.array([1.4, 1.5, 1.7, 1.9, 2.1, 2.3, 2.4, 2.5, 2.6, 2.8, 2.8, 2.9,
2.9, 3.0, 3.1, 3.1, 3.2, 3.2, 3.3, 3.3, 3.4, 3.4, 3.5, 3.5,
3.6, 3.7, 3.7, 3.8, 3.8, 3.9, 3.9, 4.0, 4.0, 4.1, 4.2, 4.2,
4.2, 4.2, 4.3, 4.3, 4.3, 4.4, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9,
5.1, 5.2, 5.4, 5.5, 5.7])
# prepare spline interpolation to calculate return value # prepare spline interpolation to calculate return value
func, params = fit_curve(relation[:, 0], relation[:, 1]) func, params = fit_curve(distance, richter_scaling)
return func(delta, params) return func(delta, params)
@ -31,10 +38,10 @@ class Magnitude(object):
def __init__(self, stream, event, verbosity=False, iplot=0): def __init__(self, stream, event, verbosity=False, iplot=0):
self._type = "M" self._type = "M"
self._stream = stream
self._plot_flag = iplot self._plot_flag = iplot
self._verbosity = verbosity self._verbosity = verbosity
self._event = event self._event = event
self._stream = stream
self._magnitudes = dict() self._magnitudes = dict()
def __str__(self): def __str__(self):
@ -129,7 +136,7 @@ class Magnitude(object):
return None return None
class RichterMagnitude(Magnitude): class LocalMagnitude(Magnitude):
""" """
Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- Method to derive peak-to-peak amplitude as seen on a Wood-Anderson-
seismograph. Has to be derived from instrument corrected traces! seismograph. Has to be derived from instrument corrected traces!
@ -146,10 +153,11 @@ class RichterMagnitude(Magnitude):
_amplitudes = dict() _amplitudes = dict()
def __init__(self, stream, event, calc_win, verbosity=False, iplot=0): def __init__(self, stream, event, calc_win, wascaling=None, verbosity=False, iplot=0):
super(RichterMagnitude, self).__init__(stream, event, verbosity, iplot) super(LocalMagnitude, self).__init__(stream, event, verbosity, iplot)
self._calc_win = calc_win self._calc_win = calc_win
self._wascaling = wascaling
self._type = 'ML' self._type = 'ML'
self.calc() self.calc()
@ -161,6 +169,10 @@ class RichterMagnitude(Magnitude):
def calc_win(self, value): def calc_win(self, value):
self._calc_win = value self._calc_win = value
@property
def wascaling(self):
return self._wascaling
@property @property
def amplitudes(self): def amplitudes(self):
return self._amplitudes return self._amplitudes
@ -244,10 +256,16 @@ class RichterMagnitude(Magnitude):
self.event.amplitudes.append(amplitude) self.event.amplitudes.append(amplitude)
self.amplitudes = (station, amplitude) self.amplitudes = (station, amplitude)
# using standard Gutenberg-Richter relation # using standard Gutenberg-Richter relation
# TODO make the ML calculation more flexible by allowing # or scale WA amplitude with given scaling relation
# use of custom relation functions if self.wascaling == None:
magnitude = ope.StationMagnitude( print("Calculating original Richter magnitude ...")
mag=np.log10(a0) + richter_magnitude_scaling(delta)) magnitude = ope.StationMagnitude(mag=np.log10(a0) \
+ richter_magnitude_scaling(delta))
else:
print("Calculating scaled local magnitude ...")
magnitude = ope.StationMagnitude(mag=np.log10(a0) \
+ self.wascaling[0] * np.log10(delta) + self.wascaling[1]
* delta + self.wascaling[2])
magnitude.origin_id = self.origin_id magnitude.origin_id = self.origin_id
magnitude.waveform_id = pick.waveform_id magnitude.waveform_id = pick.waveform_id
magnitude.amplitude_id = amplitude.resource_id magnitude.amplitude_id = amplitude.resource_id

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

@ -5,6 +5,7 @@ import copy
import os import os
from obspy import read_events from obspy import read_events
from obspy.core import read, Stream, UTCDateTime from obspy.core import read, Stream, UTCDateTime
from obspy.io.sac import SacIOError
from obspy.core.event import Event from obspy.core.event import Event
from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \ from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \
picksdict_from_pilot, merge_picks picksdict_from_pilot, merge_picks
@ -230,6 +231,8 @@ class Data(object):
self.wfdata += read(fname, format='GSE2') self.wfdata += read(fname, format='GSE2')
except Exception as e: except Exception as e:
warnmsg += '{0}\n{1}\n'.format(fname, e) warnmsg += '{0}\n{1}\n'.format(fname, e)
except SacIOError as se:
warnmsg += '{0}\n{1}\n'.format(fname, se)
if warnmsg: if warnmsg:
warnmsg = 'WARNING: unable to read\n' + warnmsg warnmsg = 'WARNING: unable to read\n' + warnmsg
print(warnmsg) print(warnmsg)

13
pylot/core/io/default_parameters.py
View File

@ -275,7 +275,15 @@ defaults = {'rootpath': {'type': str,
'wdttolerance': {'type': float, 'wdttolerance': {'type': float,
'tooltip': 'maximum allowed deviation from Wadati-diagram', 'tooltip': 'maximum allowed deviation from Wadati-diagram',
'value': 1.0} 'value': 1.0},
'WAscaling': {'type': (float, float, float),
'tooltip': 'Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm]',
'value': (1.0, 1.0, 1.0)},
'magscaling': {'type': (float, float),
'tooltip': 'Scaling relation for derived local magnitude [a*Ml+b]',
'value': (1.0, 1.0)}
} }
settings_main={ settings_main={
@ -298,6 +306,9 @@ settings_main={
'vp', 'vp',
'rho', 'rho',
'Qp'], 'Qp'],
'localmag':[
'WAscaling',
'magscaling'],
'pick':[ 'pick':[
'extent', 'extent',
'pstart', 'pstart',

11
pylot/core/io/inputs.py
View File

@ -4,9 +4,9 @@
from pylot.core.util.errors import ParameterError from pylot.core.util.errors import ParameterError
import default_parameters import default_parameters
class AutoPickParameter(object): class PylotParameter(object):
''' '''
AutoPickParameters is a parameter type object capable to read and/or write PylotParameter is a parameter type object capable to read and/or write
parameter ASCII. parameter ASCII.
:param fn str: Filename of the input file :param fn str: Filename of the input file
@ -78,7 +78,7 @@ class AutoPickParameter(object):
# String representation of the object # String representation of the object
def __repr__(self): def __repr__(self):
return "AutoPickParameter('%s')" % self.__filename return "PylotParameter('%s')" % self.__filename
# Boolean test # Boolean test
def __nonzero__(self): def __nonzero__(self):
@ -140,7 +140,8 @@ class AutoPickParameter(object):
all_names += self.get_main_para_names()['dirs'] all_names += self.get_main_para_names()['dirs']
all_names += self.get_main_para_names()['nlloc'] all_names += self.get_main_para_names()['nlloc']
all_names += self.get_main_para_names()['smoment'] all_names += self.get_main_para_names()['smoment']
all_names += self.get_main_para_names()['pick'] all_names += self.get_main_para_names()['localmag']
all_names += self.get_main_para_names()['pick']
all_names += self.get_special_para_names()['z'] all_names += self.get_special_para_names()['z']
all_names += self.get_special_para_names()['h'] all_names += self.get_special_para_names()['h']
all_names += self.get_special_para_names()['fm'] all_names += self.get_special_para_names()['fm']
@ -234,6 +235,8 @@ class AutoPickParameter(object):
'NLLoc settings', seperator) 'NLLoc settings', seperator)
self.write_section(fid_out, self.get_main_para_names()['smoment'], self.write_section(fid_out, self.get_main_para_names()['smoment'],
'parameters for seismic moment estimation', seperator) 'parameters for seismic moment estimation', seperator)
self.write_section(fid_out, self.get_main_para_names()['localmag'],
'settings local magnitude', seperator)
self.write_section(fid_out, self.get_main_para_names()['pick'], self.write_section(fid_out, self.get_main_para_names()['pick'],
'common settings picker', seperator) 'common settings picker', seperator)
fid_out.write(('#special settings for calculating CF#\n'+ fid_out.write(('#special settings for calculating CF#\n'+

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

@ -8,7 +8,7 @@ import scipy.io as sio
import warnings import warnings
from obspy.core import UTCDateTime from obspy.core import UTCDateTime
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import PylotParameter
from pylot.core.io.location import create_arrival, create_event, \ from pylot.core.io.location import create_arrival, create_event, \
create_magnitude, create_origin, create_pick create_magnitude, create_origin, create_pick
from pylot.core.pick.utils import select_for_phase from pylot.core.pick.utils import select_for_phase
@ -116,7 +116,7 @@ def picksdict_from_pilot(fn):
picks = dict() picks = dict()
phases_pilot = sio.loadmat(fn) phases_pilot = sio.loadmat(fn)
stations = stations_from_pilot(phases_pilot['stat']) stations = stations_from_pilot(phases_pilot['stat'])
params = AutoPickParameter(TIMEERROR_DEFAULTS) params = PylotParameter(TIMEERROR_DEFAULTS)
timeerrors = dict(P=params.get('timeerrorsP'), timeerrors = dict(P=params.get('timeerrorsP'),
S=params.get('timeerrorsS')) S=params.get('timeerrorsS'))
for n, station in enumerate(stations): for n, station in enumerate(stations):
@ -295,14 +295,14 @@ def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0
def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None, verbosity=0): def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None, verbosity=0):
from obspy import read from obspy import read
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import PylotParameter
from pylot.core.pick.utils import earllatepicker from pylot.core.pick.utils import earllatepicker
if fn_param is None: if fn_param is None:
import pylot.core.util.defaults as defaults import pylot.core.util.defaults as defaults
fn_param = defaults.AUTOMATIC_DEFAULTS fn_param = defaults.AUTOMATIC_DEFAULTS
default = AutoPickParameter(fn_param, verbosity) default = PylotParameter(fn_param, verbosity)
search_base = os.path.join(root_dir, db_dir, event_id) search_base = os.path.join(root_dir, db_dir, event_id)
phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat')) phases_file = glob.glob(os.path.join(search_base, 'PHASES.mat'))
@ -358,7 +358,7 @@ def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None
default.get('tsnrz' if phase == 'P' else 'tsnrh'), default.get('tsnrz' if phase == 'P' else 'tsnrh'),
Pick1=rel_pick, Pick1=rel_pick,
iplot=None, iplot=None,
stealth_mode=True) verbosity=0)
if epp is None or lpp is None: if epp is None or lpp is None:
continue continue
epp = stime + epp epp = stime + epp

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

@ -11,7 +11,7 @@ function conglomerate utils.
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import PylotParameter
from pylot.core.pick.picker import AICPicker, PragPicker from pylot.core.pick.picker import AICPicker, PragPicker
from pylot.core.pick.charfuns import CharacteristicFunction from pylot.core.pick.charfuns import CharacteristicFunction
from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
@ -81,7 +81,7 @@ def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None):
:param pickparam: container of picking parameters from input file, :param pickparam: container of picking parameters from input file,
usually autoPyLoT.in usually autoPyLoT.in
:type pickparam: AutoPickParameter :type pickparam: PylotParameter
:param verbose: :param verbose:
:type verbose: bool :type verbose: bool

51
pylot/core/pick/picker.py
View File

@ -223,26 +223,31 @@ class AICPicker(AutoPicker):
# find maximum within slope determination window # find maximum within slope determination window
# 'cause slope should be calculated up to first local minimum only! # 'cause slope should be calculated up to first local minimum only!
imax = np.argmax(self.Data[0].data[islope]) imax = np.argmax(self.Data[0].data[islope])
if imax == 0: iislope = islope[0][0:imax]
print('AICPicker: Maximum for slope determination right at the beginning of the window!') if len(iislope) <= 2:
print('Choose longer slope determination window!') # calculate slope from initial onset to maximum of AIC function
if self.iplot > 1: print("AICPicker: Not enough data samples left for slope calculation!")
if not self.fig: print("Calculating slope from initial onset to maximum of AIC function ...")
fig = plt.figure() #self.iplot) ### WHY? MP MP imax = np.argmax(aicsmooth[islope])
else: if imax == 0:
fig = self.fig print("AICPicker: Maximum for slope determination right at the beginning of the window!")
ax = fig.add_subplot(111) print("Choose longer slope determination window!")
x = self.Data[0].data if self.iplot > 1:
ax.plot(self.Tcf, x / max(x), 'k', legend='(HOS-/AR-) Data') if not self.fig:
ax.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', legend='Smoothed AIC-CF') fig = plt.figure() #self.iplot) ### WHY? MP MP
ax.legend() else:
ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime) fig = self.fig
ax.set_yticks([]) ax = fig.add_subplot(111)
ax.set_title(self.Data[0].stats.station) x = self.Data[0].data
return ax.plot(self.Tcf, x / max(x), 'k', label='(HOS-/AR-) Data')
ax.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
islope = islope[0][0:imax] ax.legend()
dataslope = self.Data[0].data[islope] ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax.set_yticks([])
ax.set_title(self.Data[0].stats.station)
return
iislope = islope[0][0:imax]
dataslope = self.Data[0].data[iislope]
# 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)
P = np.polyfit(xslope, dataslope, 1) P = np.polyfit(xslope, dataslope, 1)
@ -276,12 +281,12 @@ class AICPicker(AutoPicker):
ax2.plot(self.Tcf, x, 'k', label='Data') ax2.plot(self.Tcf, x, 'k', label='Data')
ax1.axvspan(self.Tcf[inoise[0]],self.Tcf[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window') ax1.axvspan(self.Tcf[inoise[0]],self.Tcf[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
ax1.axvspan(self.Tcf[isignal[0]],self.Tcf[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window') ax1.axvspan(self.Tcf[isignal[0]],self.Tcf[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window')
ax1.axvspan(self.Tcf[islope[0]],self.Tcf[islope[-1]], color='g', alpha=0.2, lw=0, label='Slope Window') ax1.axvspan(self.Tcf[iislope[0]],self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0, label='Slope Window')
ax2.axvspan(self.Tcf[inoise[0]],self.Tcf[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window') ax2.axvspan(self.Tcf[inoise[0]],self.Tcf[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
ax2.axvspan(self.Tcf[isignal[0]],self.Tcf[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window') ax2.axvspan(self.Tcf[isignal[0]],self.Tcf[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window')
ax2.axvspan(self.Tcf[islope[0]],self.Tcf[islope[-1]], color='g', alpha=0.2, lw=0, label='Slope Window') ax2.axvspan(self.Tcf[iislope[0]],self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0, label='Slope Window')
ax2.plot(self.Tcf[islope], datafit, 'g', linewidth=2, label='Slope') ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2, label='Slope')
ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope)) self.SNR, self.slope))

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

@ -14,7 +14,7 @@ import numpy as np
from obspy.core import Stream, UTCDateTime from obspy.core import Stream, UTCDateTime
def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False, fig=None): def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, verbosity=1, fig=None):
''' '''
Function to derive earliest and latest possible pick after Diehl & Kissling (2009) Function to derive earliest and latest possible pick after Diehl & Kissling (2009)
as reasonable uncertainties. Latest possible pick is based on noise level, as reasonable uncertainties. Latest possible pick is based on noise level,
@ -40,10 +40,16 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False, fig=Non
assert isinstance(X, Stream), "%s is not a stream object" % str(X) assert isinstance(X, Stream), "%s is not a stream object" % str(X)
if verbosity == 2:
print('earllatepicker:')
print('nfac:', nfac)
print('Init pick:', Pick1)
print('TSNR (T_noise, T_gap, T_signal):', TSNR)
LPick = None LPick = None
EPick = None EPick = None
PickError = None PickError = None
if stealth_mode is False: if verbosity:
print('earllatepicker: Get earliest and latest possible pick' print('earllatepicker: Get earliest and latest possible pick'
' relative to most likely pick ...') ' relative to most likely pick ...')
@ -57,11 +63,18 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False, fig=Non
x = x - np.mean(x[inoise]) x = x - np.mean(x[inoise])
# calculate noise level # calculate noise level
nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac
if verbosity == 2:
print('x:', x)
print('t:', t)
print('x_inoise:', x[inoise])
print('x_isignal:', x[isignal])
print('nlevel:', nlevel)
# get time where signal exceeds nlevel # get time where signal exceeds nlevel
ilup, = np.where(x[isignal] > nlevel) ilup, = np.where(x[isignal] > nlevel)
ildown, = np.where(x[isignal] < -nlevel) ildown, = np.where(x[isignal] < -nlevel)
if not ilup.size and not ildown.size: if not ilup.size and not ildown.size:
if stealth_mode is False: if verbosity:
print ("earllatepicker: Signal lower than noise level!\n" print ("earllatepicker: Signal lower than noise level!\n"
"Skip this trace!") "Skip this trace!")
return LPick, EPick, PickError return LPick, EPick, PickError
@ -78,7 +91,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False, fig=Non
# if EPick stays NaN the signal window size will be doubled # if EPick stays NaN the signal window size will be doubled
while np.isnan(EPick): while np.isnan(EPick):
if count > 0: if count > 0:
if stealth_mode is False: if verbosity:
print("\nearllatepicker: Doubled signal window size %s time(s) " print("\nearllatepicker: Doubled signal window size %s time(s) "
"because of NaN for earliest pick." % count) "because of NaN for earliest pick." % count)
isigDoubleWinStart = pis[-1] + 1 isigDoubleWinStart = pis[-1] + 1
@ -87,7 +100,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False, fig=Non
if (isigDoubleWinStart + len(pis)) < X[0].data.size: if (isigDoubleWinStart + len(pis)) < X[0].data.size:
pis = np.concatenate((pis, isignalDoubleWin)) pis = np.concatenate((pis, isignalDoubleWin))
else: else:
if stealth_mode is False: if verbosity:
print("Could not double signal window. Index out of bounds.") print("Could not double signal window. Index out of bounds.")
break break
count += 1 count += 1

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

@ -169,7 +169,8 @@ def read_metadata(path_to_inventory):
dlfile = list() dlfile = list()
invfile = list() invfile = list()
respfile = list() respfile = list()
inv = dict(dless=dlfile, xml=invfile, resp=respfile) # possible file extensions specified here:
inv = dict(dless=dlfile, xml=invfile, resp=respfile, dseed=dlfile)
if os.path.isfile(path_to_inventory): if os.path.isfile(path_to_inventory):
ext = os.path.splitext(path_to_inventory)[1].split('.')[1] ext = os.path.splitext(path_to_inventory)[1].split('.')[1]
inv[ext] += [path_to_inventory] inv[ext] += [path_to_inventory]

66
pylot/core/util/defaults.py
View File

@ -56,27 +56,55 @@ OUTPUTFORMATS = {'.xml': 'QUAKEML',
LOCTOOLS = dict(nll=nll, hyposat=hyposat, velest=velest, hypo71=hypo71, hypodd=hypodd) LOCTOOLS = dict(nll=nll, hyposat=hyposat, velest=velest, hypo71=hypo71, hypodd=hypodd)
class SetChannelComponents: class SetChannelComponents(object):
def getdefaultCompPosition(): def __init__(self):
self.setDefaultCompPosition()
def setDefaultCompPosition(self):
# default component order # default component order
CompPosition_Map = dict(Z=2, N=1, E=0) self.compPosition_Map = dict(Z=2, N=1, E=0)
CompPosition_Map['1'] = 1 self.compName_Map = {'3': 'Z',
CompPosition_Map['2'] = 0 '1': 'N',
CompPosition_Map['3'] = 2 '2': 'E'}
CompName_Map = dict(Z='3', N='1', E='2')
CompName_Map['1'] = str(1) def _getCurrentPosition(self, component):
CompName_Map['2'] = str(2) for key, value in self.compName_Map.items():
CompName_Map['3'] = str(3) if value == component:
return CompPosition_Map, CompName_Map return key, value
errMsg = 'getCurrentPosition: Could not find former position of component {}.'.format(component)
raise ValueError(errMsg)
CompPosition_Map, CompName_Map = getdefaultCompPosition() def _switch(self, component, component_alter):
# Without switching, multiple definitions of the same alter_comp are possible
old_alter_comp, _ = self._getCurrentPosition(component)
old_comp = self.compName_Map[component_alter]
if not old_alter_comp == component_alter and not old_comp == component:
self.compName_Map[old_alter_comp] = old_comp
print('switch: Automatically switched component {} to {}'.format(old_alter_comp, old_comp))
def setCompPosition(self, component, position): def setCompPosition(self, component_alter, component, switch=True):
self.CompPosition_Map[component] = position component_alter = str(component_alter)
self.CompName_Map[component] = str(position) if not component_alter in self.compName_Map.keys():
errMsg='setCompPosition: Unrecognized alternative component {}. Expecting one of {}.'
raise ValueError(errMsg.format(component_alter, self.compName_Map.keys()))
if not component in self.compPosition_Map.keys():
errMsg='setCompPosition: Unrecognized target component {}. Expecting one of {}.'
raise ValueError(errMsg.format(component, self.compPosition_Map.keys()))
print('setCompPosition: set component {} to {}'.format(component_alter, component))
if switch:
self._switch(component, component_alter)
self.compName_Map[component_alter] = component
def getCompPosition(self, component): def getCompPosition(self, component):
self.comppos = self.CompPosition_Map[component] return self._getCurrentPosition(component)[0]
self.compname = self.CompName_Map[component]
return self.comppos, self.compname def getPlotPosition(self, component):
component = str(component)
if component in self.compPosition_Map.keys():
return self.compPosition_Map[component]
elif component in self.compName_Map.keys():
return self.compPosition_Map[self.compName_Map[component]]
else:
errMsg='getCompPosition: Unrecognized component {}. Expecting one of {} or {}.'
raise ValueError(errMsg.format(component, self.compPosition_Map.keys(), self.compName_Map.keys()))

12
pylot/core/util/map_projection.py
View File

@ -43,7 +43,7 @@ class map_projection(QtGui.QWidget):
return return
data = self._parent.get_data().getWFData() data = self._parent.get_data().getWFData()
for index in ind: for index in ind:
station=str(self.station_names[index]) station=str(self.station_names[index].split('.')[-1])
try: try:
pickDlg = PickDlg(self, parameter=self._parent._inputs, pickDlg = PickDlg(self, parameter=self._parent._inputs,
data=data.select(station=station), data=data.select(station=station),
@ -51,7 +51,7 @@ class map_projection(QtGui.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))
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)
self._warn(message) self._warn(message)
print(message, e) print(message, e)
return return
@ -120,8 +120,9 @@ class map_projection(QtGui.QWidget):
lon=[] lon=[]
for station in parser.stations: for station in parser.stations:
station_name=station[0].station_call_letters station_name=station[0].station_call_letters
network=station[0].network_code
if not station_name in station_names: if not station_name in station_names:
station_names.append(station_name) station_names.append(network+'.'+station_name)
lat.append(station[0].latitude) lat.append(station[0].latitude)
lon.append(station[0].longitude) lon.append(station[0].longitude)
return station_names, lat, lon return station_names, lat, lon
@ -137,6 +138,7 @@ class map_projection(QtGui.QWidget):
picks=[] picks=[]
for station in station_names: for station in station_names:
try: try:
station=station.split('.')[-1]
picks.append(self.picks_dict[station][phase]['mpp']) picks.append(self.picks_dict[station][phase]['mpp'])
except: except:
picks.append(np.nan) picks.append(np.nan)
@ -234,9 +236,9 @@ class map_projection(QtGui.QWidget):
self.picks_no_nan, (self.latgrid, self.longrid), method='linear') ################## self.picks_no_nan, (self.latgrid, self.longrid), method='linear') ##################
def draw_contour_filled(self, nlevel='50'): def draw_contour_filled(self, nlevel='50'):
levels = np.linspace(min(self.picks_rel), max(self.picks_rel), nlevel) levels = np.linspace(min(self.picks_no_nan), max(self.picks_no_nan), nlevel)
self.contourf = self.basemap.contourf(self.longrid, self.latgrid, self.picksgrid_no_nan, self.contourf = self.basemap.contourf(self.longrid, self.latgrid, self.picksgrid_no_nan,
levels, latlon=True, zorder=9) levels, latlon=True, zorder=9, alpha=0.5)
def scatter_all_stations(self): def scatter_all_stations(self):
self.sc = self.basemap.scatter(self.lon, self.lat, s=50, facecolor='none', latlon=True, self.sc = self.basemap.scatter(self.lon, self.lat, s=50, facecolor='none', latlon=True,

6
pylot/core/util/thread.py
View File

@ -1,5 +1,5 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import sys import sys, os
from PySide.QtCore import QThread, Signal, Qt from PySide.QtCore import QThread, Signal, Qt
from PySide.QtGui import QDialog, QProgressBar, QLabel, QHBoxLayout from PySide.QtGui import QDialog, QProgressBar, QLabel, QHBoxLayout
@ -64,7 +64,9 @@ class Thread(QThread):
except Exception as e: except Exception as e:
self._executed = False self._executed = False
self._executedError = e self._executedError = e
print(e) exc_type, exc_obj, exc_tb = sys.exc_info()
fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
print('Exception: {}, file: {}, line: {}'.format(exc_type, fname, exc_tb.tb_lineno))
sys.stdout = sys.__stdout__ sys.stdout = sys.__stdout__
def __del__(self): def __del__(self):

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

@ -10,7 +10,7 @@ import re
import warnings import warnings
import subprocess import subprocess
from obspy import UTCDateTime, read from obspy import UTCDateTime, read
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import PylotParameter
def _pickle_method(m): def _pickle_method(m):
@ -497,7 +497,7 @@ def which(program, infile=None):
bpath = os.path.join(os.path.expanduser('~'), '.pylot', infile) bpath = os.path.join(os.path.expanduser('~'), '.pylot', infile)
if os.path.exists(bpath): if os.path.exists(bpath):
nllocpath = ":" + AutoPickParameter(bpath).get('nllocbin') nllocpath = ":" + PylotParameter(bpath).get('nllocbin')
os.environ['PATH'] += nllocpath os.environ['PATH'] += nllocpath
except ImportError as e: except ImportError as e:
print(e.message) print(e.message)

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

File diff suppressed because it is too large Load Diff