Merge remote-tracking branch 'origin/develop' into feature/refactor

This commit is contained in:
Darius Arnold
2018-07-06 09:49:52 +02:00
19 changed files with 1382 additions and 640 deletions

View File

@@ -6,6 +6,7 @@ Revised/extended summer 2017.
:author: Ludger Küperkoch / MAGS2 EP3 working group
"""
import matplotlib.pyplot as plt
import numpy as np
import obspy.core.event as ope
@@ -121,7 +122,13 @@ class Magnitude(object):
def net_magnitude(self, magscaling=None):
if self:
if magscaling is not None and str(magscaling) is not '[0.0, 0.0]':
if magscaling == None:
scaling = False
elif magscaling[0] != 0 and magscaling[1] != 0:
scaling = False
else:
scaling = True
if scaling:
# scaling necessary
print("Scaling network magnitude ...")
mag = ope.Magnitude(
@@ -140,7 +147,6 @@ class Magnitude(object):
station_count=len(self.magnitudes),
azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap)
return mag
return None
class LocalMagnitude(Magnitude):
@@ -219,7 +225,7 @@ class LocalMagnitude(Magnitude):
sqH = np.sqrt(power_sum)
# get time array
th = np.arange(0, len(sqH) * dt, dt)
th=np.arange(0, st[0].stats.npts/st[0].stats.sampling_rate, st[0].stats.delta)
# get maximum peak within pick window
iwin = getsignalwin(th, t0 - stime, self.calc_win)
ii = min([iwin[len(iwin) - 1], len(th)])

View File

@@ -15,7 +15,7 @@ from pylot.core.util.event import Event
from pylot.core.util.utils import fnConstructor, full_range, remove_underscores, check4gaps, check4doubled, \
check4rotated, trim_station_components
import pylot.core.loc.velest as velest
from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
class Data(object):
"""
@@ -43,7 +43,7 @@ class Data(object):
elif isinstance(evtdata, dict):
evt = readPILOTEvent(**evtdata)
evtdata = evt
elif isinstance(evtdata, str):
elif type(evtdata) in [str, unicode]:
try:
cat = read_events(evtdata)
if len(cat) is not 1:
@@ -60,6 +60,8 @@ class Data(object):
raise NotImplementedError('PILOT location information '
'read support not yet '
'implemeted.')
elif 'event.pkl' in evtdata:
evtdata = qml_from_obspyDMT(evtdata)
else:
raise e
else:
@@ -72,6 +74,7 @@ class Data(object):
self.wforiginal = None
self.cuttimes = None
self.dirty = False
self.processed = None
def __str__(self):
return str(self.wfdata)
@@ -368,7 +371,7 @@ class Data(object):
data.filter(**kwargs)
self.dirty = True
def setWFData(self, fnames, checkRotated=False, metadata=None, obspy_dmt=False):
def setWFData(self, fnames, fnames_syn=None, checkRotated=False, metadata=None):
"""
Clear current waveform data and set given waveform data
:param fnames: waveform data names to append
@@ -377,30 +380,31 @@ class Data(object):
self.wfdata = Stream()
self.wforiginal = None
self.wfsyn = Stream()
wffnames = None
wffnames_syn = None
wfdir = 'processed' if 'processed' in [fname.split('/')[-1] for fname in fnames] else 'raw'
if obspy_dmt:
for fpath in fnames:
if fpath.endswith(wfdir):
wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
if 'syngine' in fpath.split('/')[-1]:
wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
else:
wffnames = fnames
if wffnames is not None:
self.appendWFData(wffnames)
if wffnames_syn is not None:
self.appendWFData(wffnames_syn, synthetic=True)
# if obspy_dmt:
# wfdir = 'raw'
# self.processed = False
# for fname in fnames:
# if fname.endswith('processed'):
# wfdir = 'processed'
# self.processed = True
# break
# for fpath in fnames:
# if fpath.endswith(wfdir):
# wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
# if 'syngine' in fpath.split('/')[-1]:
# wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
# else:
# wffnames = fnames
if fnames is not None:
self.appendWFData(fnames)
if fnames_syn is not None:
self.appendWFData(fnames_syn, synthetic=True)
else:
return False
# various pre-processing steps:
# remove possible underscores in station names
self.wfdata = remove_underscores(self.wfdata)
# check for gaps and doubled channels
check4gaps(self.wfdata)
check4doubled(self.wfdata)
# check for stations with rotated components
if checkRotated and metadata is not None:
self.wfdata = check4rotated(self.wfdata, metadata, verbosity=0)
@@ -441,7 +445,7 @@ class Data(object):
except SacIOError as se:
warnmsg += '{0}\n{1}\n'.format(fname, se)
if warnmsg:
warnmsg = 'WARNING: unable to read\n' + warnmsg
warnmsg = 'WARNING in appendWFData: unable to read waveform data\n' + warnmsg
print(warnmsg)
def getWFData(self):

View File

@@ -100,9 +100,9 @@ def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None,
# quality control
# median check and jackknife on P-onset times
jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, 1, fig_dict_wadatijack)
jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, iplot, fig_dict_wadatijack)
# check S-P times (Wadati)
wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, 1, fig_dict_wadatijack)
wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, iplot, fig_dict_wadatijack)
return wadationsets
@@ -1057,8 +1057,64 @@ def nautopickstation(wfstream, pickparam, verbose=False,
# read your pylot.in for details!
plt_flag = 0
# get picking parameter dictionaries
p_params, s_params, first_motion_params, signal_length_params = get_pickparams(pickparam)
# special parameters for P picking
algoP = pickparam.get('algoP')
pstart = pickparam.get('pstart')
pstop = pickparam.get('pstop')
thosmw = pickparam.get('tlta')
tsnrz = pickparam.get('tsnrz')
hosorder = pickparam.get('hosorder')
bpz1 = pickparam.get('bpz1')
bpz2 = pickparam.get('bpz2')
pickwinP = pickparam.get('pickwinP')
aictsmoothP = pickparam.get('aictsmooth')
tsmoothP = pickparam.get('tsmoothP')
ausP = pickparam.get('ausP')
nfacP = pickparam.get('nfacP')
tpred1z = pickparam.get('tpred1z')
tdet1z = pickparam.get('tdet1z')
tpred2z = pickparam.get('tpred2z')
tdet2z = pickparam.get('tdet2z')
Parorder = pickparam.get('Parorder')
addnoise = pickparam.get('addnoise')
Precalcwin = pickparam.get('Precalcwin')
minAICPslope = pickparam.get('minAICPslope')
minAICPSNR = pickparam.get('minAICPSNR')
timeerrorsP = pickparam.get('timeerrorsP')
# special parameters for S picking
algoS = pickparam.get('algoS')
sstart = pickparam.get('sstart')
sstop = pickparam.get('sstop')
use_taup = real_Bool(pickparam.get('use_taup'))
taup_model = pickparam.get('taup_model')
bph1 = pickparam.get('bph1')
bph2 = pickparam.get('bph2')
tsnrh = pickparam.get('tsnrh')
pickwinS = pickparam.get('pickwinS')
tpred1h = pickparam.get('tpred1h')
tdet1h = pickparam.get('tdet1h')
tpred2h = pickparam.get('tpred2h')
tdet2h = pickparam.get('tdet2h')
Sarorder = pickparam.get('Sarorder')
aictsmoothS = pickparam.get('aictsmoothS')
tsmoothS = pickparam.get('tsmoothS')
ausS = pickparam.get('ausS')
minAICSslope = pickparam.get('minAICSslope')
minAICSSNR = pickparam.get('minAICSSNR')
Srecalcwin = pickparam.get('Srecalcwin')
nfacS = pickparam.get('nfacS')
timeerrorsS = pickparam.get('timeerrorsS')
# parameters for first-motion determination
minFMSNR = pickparam.get('minFMSNR')
fmpickwin = pickparam.get('fmpickwin')
minfmweight = pickparam.get('minfmweight')
# parameters for checking signal length
minsiglength = pickparam.get('minsiglength')
minpercent = pickparam.get('minpercent')
nfacsl = pickparam.get('noisefactor')
# parameter to check for spuriously picked S onset
zfac = pickparam.get('zfac')
# path to inventory-, dataless- or resp-files
# initialize output
Pweight = 4 # weight for P onset
@@ -1159,6 +1215,8 @@ def nautopickstation(wfstream, pickparam, verbose=False,
if p_params['use_taup'] is True:
Lc = np.inf
print('autopickstation: use_taup flag active.')
if not metadata:
metadata = [None, None]
if not metadata[1]:
print('Warning: Could not use TauPy to estimate onsets as there are no metadata given.')
else:
@@ -1207,7 +1265,7 @@ def nautopickstation(wfstream, pickparam, verbose=False,
return
Ldiff = Lwf - abs(Lc)
if Ldiff < 0 or pstop <= pstart:
if Ldiff <= 0 or pstop <= pstart or pstop - pstart <= thosmw:
msg = 'autopickstation: Cutting times are too large for actual ' \
'waveform!\nUsing entire waveform instead!'
if verbose: print(msg)
@@ -1345,8 +1403,8 @@ def nautopickstation(wfstream, pickparam, verbose=False,
elif p_params['algoP'] == 'ARZ':
# calculate ARZ-CF using subclass ARZcf of class
# CharcteristicFunction
cf2 = ARZcf(z_copy, cuttimes2, p_params['tpred1z'], p_params['Parorder'], p_params['tdet1z'],
p_params['addnoise']) # instance of ARZcf
cf2 = ARZcf(z_copy, cuttimes2, tpred2z, Parorder, tdet2z,
addnoise) # instance of ARZcf
##############################################################
# get refined onset time from CF2 using class Picker
assert isinstance(cf2, CharacteristicFunction), 'cf2 is not set ' \
@@ -1401,7 +1459,8 @@ def nautopickstation(wfstream, pickparam, verbose=False,
msg = "autopickstation: P-weight: {0}, " \
"SNR: {1}, SNR[dB]: {2}, Polarity: {3}".format(Pweight, SNRP, SNRPdB, FM)
print(msg)
msg = 'autopickstation: Refined P-Pick: {} s | P-Error: {} s'.format(mpickP, Perror)
msg = 'autopickstation: Refind P-Pick: {} s | P-Error: {} s'.format(zdat[0].stats.starttime \
+ mpickP, Perror)
print(msg)
Sflag = 1
@@ -1667,7 +1726,8 @@ def nautopickstation(wfstream, pickparam, verbose=False,
lpickS = lpick[ipick]
Serror = pickerr[ipick]
msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(mpickS, Serror)
msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(hdat[0].stats.starttime \
+ mpickS, Serror)
print(msg)
# get SNR

View File

@@ -78,7 +78,7 @@ class CharacteristicFunction(object):
return self.cut
def setCut(self, cut):
self.cut = cut
self.cut = (int(cut[0]), int(cut[1]))
def getTime1(self):
return self.t1

View File

@@ -233,12 +233,41 @@ class AICPicker(AutoPicker):
# quality assessment using SNR and slope from CF
if self.Pick is not None:
self.Data[0].data = check_counts_ms(self.Data[0].data)
# get noise window
inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1])
# check, if these are counts or m/s, important for slope estimation!
# this is quick and dirty, better solution? #Todo wtf
if max(self.Data[0].data < 1e-3) and max(self.Data[0].data >= 1e-6):
self.Data[0].data = self.Data[0].data * 1000000.
elif max(self.Data[0].data < 1e-6):
self.Data[0].data = self.Data[0].data * 1e13
# get signal window
isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2])
if len(isignal) == 0:
return
ii = min([isignal[len(isignal) - 1], len(self.Tcf)])
isignal = isignal[0:ii]
try:
self.Data[0].data[isignal]
except IndexError as e:
msg = "Time series out of bounds! {}".format(e)
print(msg)
return
# calculate SNR from CF
# subtract beginning of Tcf to get sample index in self.Data, which
# is cut out off trace
self.SNR = getSNR(self.Data, self.TSNR, self.Pick-self.Tcf[0])[0] # TODO Check wether this yields similar results
self.SNR = max(abs(self.Data[0].data[isignal])) / \
abs(np.mean(self.Data[0].data[inoise]))
# calculate slope from CF after initial pick
# get slope window
tslope = self.TSNR[3] # slope determination window
tsafety = self.TSNR[1] # safety gap, AIC is usually a little bit too late
if tsafety >= 0:
islope = np.where((self.Tcf <= min([self.Pick + tslope + tsafety, self.Tcf[-1]])) \
& (self.Tcf >= self.Pick)) # TODO: put this in a seperate function like getsignalwin
else:
islope = np.where((self.Tcf <= min([self.Pick + tslope, self.Tcf[-1]])) \
& (self.Tcf >= self.Pick + tsafety)) # TODO: put this in a seperate function like getsignalwin
# find maximum within slope determination window
# 'cause slope should be calculated up to first local minimum only!
try:
self.slope, iislope, datafit = calcSlope(self.Data, self.aicsmooth, self.Tcf, self.Pick, self.TSNR)
except Exception:
@@ -264,6 +293,56 @@ class AICPicker(AutoPicker):
pass
plt.close(fig)
return
if len(dataslope) < 2:
print('No or not enough data in slope window found!')
return
try:
imaxs, = argrelmax(dataslope)
imaxs.size
imax = imaxs[0]
except ValueError as e:
print(e, 'picker: argrelmax not working!')
imax = np.argmax(dataslope)
iislope = islope[0][0:imax + 1]
if len(iislope) < 2:
# calculate slope from initial onset to maximum of AIC function
print("AICPicker: Not enough data samples left for slope calculation!")
print("Calculating slope from initial onset to maximum of AIC function ...")
imax = np.argmax(aicsmooth[islope[0][0:-1]])
if imax == 0:
print("AICPicker: Maximum for slope determination right at the beginning of the window!")
print("Choose longer slope determination window!")
if self.iplot > 1:
if self.fig == None or self.fig == 'None':
fig = plt.figure()
plt_flag = 1
else:
fig = self.fig
ax = fig.add_subplot(111)
x = self.Data[0].data
ax.plot(self.Tcf, x / max(x), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
ax.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
ax.legend(loc=1)
ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax.set_yticks([])
ax.set_title(self.Data[0].stats.station)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
plt.close(fig)
return
iislope = islope[0][0:imax+1]
dataslope = self.Data[0].data[iislope]
# calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope), 1)
P = np.polyfit(xslope, dataslope, 1)
datafit = np.polyval(P, xslope)
if datafit[0] >= datafit[-1]:
print('AICPicker: Negative slope, bad onset skipped!')
else:
self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
else:
self.SNR = -1
self.slope = -1
@@ -305,24 +384,23 @@ class AICPicker(AutoPicker):
label='Slope Window')
ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2, label='Slope')
ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope))
if self.slope is not None:
ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope))
else:
ax1.set_title('Station %s, SNR=%7.2f' % (self.Data[0].stats.station, self.SNR))
ax2.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax2.set_ylabel('Counts')
ax2.set_yticks([])
ax2.legend(loc=1)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
plt.close(fig)
else:
ax1.set_title(self.Data[0].stats.station)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
plt.close(fig)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
plt.close(fig)
if self.Pick == None:
print('AICPicker: Could not find minimum, picking window too short?')

View File

@@ -0,0 +1,458 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import obspy
from PySide import QtGui
from mpl_toolkits.basemap import Basemap
from matplotlib.figure import Figure
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from pylot.core.util.widgets import PickDlg, PylotCanvas
from scipy.interpolate import griddata
plt.interactive(False)
class Array_map(QtGui.QWidget):
def __init__(self, parent, figure=None):
'''
Create a map of the array.
:param parent: PyLoT Mainwindow class
:param figure:
'''
QtGui.QWidget.__init__(self)
self._parent = parent
self.metadata_type = parent.metadata[0]
self.metadata = parent.metadata[1]
self.picks = None
self.picks_dict = None
self.autopicks_dict = None
self.eventLoc = None
self.figure = figure
self.init_graphics()
self.init_stations()
self.init_basemap(resolution='l')
self.init_map()
self._style = parent._style
# self.show()
def init_map(self):
self.init_lat_lon_dimensions()
self.init_lat_lon_grid()
self.init_x_y_dimensions()
self.connectSignals()
self.draw_everything()
self.canvas.setZoomBorders2content()
def onpick(self, event):
ind = event.ind
button = event.mouseevent.button
if ind == [] or not button == 1:
return
data = self._parent.get_data().getWFData()
for index in ind:
station = str(self._station_onpick_ids[index].split('.')[-1])
try:
data = data.select(station=station)
if not data:
self._warn('No data for station {}'.format(station))
return
pickDlg = PickDlg(self._parent, parameter=self._parent._inputs,
data=data,
station=station,
picks=self._parent.get_current_event().getPick(station),
autopicks=self._parent.get_current_event().getAutopick(station),
filteroptions=self._parent.filteroptions)
except Exception as e:
message = 'Could not generate Plot for station {st}.\n {er}'.format(st=station, er=e)
self._warn(message)
print(message, e)
return
pyl_mw = self._parent
try:
if pickDlg.exec_():
pyl_mw.setDirty(True)
pyl_mw.update_status('picks accepted ({0})'.format(station))
replot = pyl_mw.get_current_event().setPick(station, pickDlg.getPicks())
self._refresh_drawings()
if replot:
pyl_mw.plotWaveformData()
pyl_mw.drawPicks()
pyl_mw.draw()
else:
pyl_mw.drawPicks(station)
pyl_mw.draw()
else:
pyl_mw.update_status('picks discarded ({0})'.format(station))
except Exception as e:
message = 'Could not save picks for station {st}.\n{er}'.format(st=station, er=e)
self._warn(message)
print(message, e)
def connectSignals(self):
self.comboBox_phase.currentIndexChanged.connect(self._refresh_drawings)
self.comboBox_am.currentIndexChanged.connect(self._refresh_drawings)
self.canvas.mpl_connect('motion_notify_event', self.mouse_moved)
#self.zoom_id = self.basemap.ax.figure.canvas.mpl_connect('scroll_event', self.zoom)
def _from_dict(self, function, key):
return function(self.stations_dict.values(), key=lambda x: x[key])[key]
def get_min_from_stations(self, key):
return self._from_dict(min, key)
def get_max_from_stations(self, key):
return self._from_dict(max, key)
def get_min_from_picks(self):
return min(self.picks_rel.values())
def get_max_from_picks(self):
return max(self.picks_rel.values())
def mouse_moved(self, event):
if not event.inaxes == self.main_ax:
return
x = event.xdata
y = event.ydata
lat, lon = self.basemap(x, y, inverse=True)
self.status_label.setText('Latitude: {}, Longitude: {}'.format(lat, lon))
def current_picks_dict(self):
picktype = self.comboBox_am.currentText()
auto_manu = {'auto': self.autopicks_dict,
'manual': self.picks_dict}
return auto_manu[picktype]
def init_graphics(self):
if not self.figure:
self.figure = Figure()
self.status_label = QtGui.QLabel()
self.main_ax = self.figure.add_subplot(111)
self.canvas = PylotCanvas(self.figure, parent=self._parent, multicursor=True,
panZoomX=False, panZoomY=False)
self.main_box = QtGui.QVBoxLayout()
self.setLayout(self.main_box)
self.top_row = QtGui.QHBoxLayout()
self.main_box.addLayout(self.top_row, 1)
self.comboBox_phase = QtGui.QComboBox()
self.comboBox_phase.insertItem(0, 'P')
self.comboBox_phase.insertItem(1, 'S')
self.comboBox_am = QtGui.QComboBox()
self.comboBox_am.insertItem(0, 'auto')
self.comboBox_am.insertItem(1, 'manual')
self.top_row.addWidget(QtGui.QLabel('Select a phase: '))
self.top_row.addWidget(self.comboBox_phase)
self.top_row.setStretch(1, 1) # set stretch of item 1 to 1
self.top_row.addWidget(QtGui.QLabel('Pick type: '))
self.top_row.addWidget(self.comboBox_am)
self.top_row.setStretch(3, 1) # set stretch of item 1 to 1
self.main_box.addWidget(self.canvas, 1)
self.main_box.addWidget(self.status_label, 0)
def init_stations(self):
def stat_info_from_parser(parser):
stations_dict = {}
for station in parser.stations:
station_name = station[0].station_call_letters
network_name = station[0].network_code
if not station_name in stations_dict.keys():
st_id = network_name + '.' + station_name
stations_dict[st_id] = {'latitude': station[0].latitude,
'longitude': station[0].longitude}
return stations_dict
def stat_info_from_inventory(inventory):
stations_dict = {}
for network in inventory.networks:
for station in network.stations:
station_name = station.code
network_name = network_name.code
if not station_name in stations_dict.keys():
st_id = network_name + '.' + station_name
stations_dict[st_id] = {'latitude': station[0].latitude,
'longitude': station[0].longitude}
return stations_dict
read_stat = {'xml': stat_info_from_inventory,
'dless': stat_info_from_parser}
self.stations_dict = read_stat[self.metadata_type](self.metadata)
self.latmin = self.get_min_from_stations('latitude')
self.lonmin = self.get_min_from_stations('longitude')
self.latmax = self.get_max_from_stations('latitude')
self.lonmax = self.get_max_from_stations('longitude')
def init_picks(self):
def get_picks(station_dict):
picks = {}
phase = self.comboBox_phase.currentText()
for st_id in station_dict.keys():
try:
station_name = st_id.split('.')[-1]
pick = self.current_picks_dict()[station_name][phase]
if pick['picker'] == 'auto':
if pick['weight'] > 3:
continue
picks[st_id] = pick['mpp']
except KeyError:
continue
except Exception as e:
print('Cannot display pick for station {}. Reason: {}'.format(station_name, e))
return picks
def get_picks_rel(picks):
picks_rel = {}
picks_utc = []
for pick in picks.values():
if type(pick) is obspy.core.utcdatetime.UTCDateTime:
picks_utc.append(pick)
self._earliest_picktime = min(picks_utc)
for st_id, pick in picks.items():
if type(pick) is obspy.core.utcdatetime.UTCDateTime:
pick -= self._earliest_picktime
picks_rel[st_id] = pick
return picks_rel
self.picks = get_picks(self.stations_dict)
self.picks_rel = get_picks_rel(self.picks)
def init_lat_lon_dimensions(self):
# init minimum and maximum lon and lat dimensions
self.londim = self.lonmax - self.lonmin
self.latdim = self.latmax - self.latmin
def init_x_y_dimensions(self):
# transformation of lat/lon to ax coordinate system
for st_id, coords in self.stations_dict.items():
lat, lon = coords['latitude'], coords['longitude']
coords['x'], coords['y'] = self.basemap(lon, lat)
self.xdim = self.get_max_from_stations('x') - self.get_min_from_stations('x')
self.ydim = self.get_max_from_stations('y') - self.get_min_from_stations('y')
def init_basemap(self, resolution='l'):
# basemap = Basemap(projection=projection, resolution = resolution, ax=self.main_ax)
width = 5e6
height = 2e6
basemap = Basemap(projection='lcc', resolution=resolution, ax=self.main_ax,
width=width, height=height,
lat_0=(self.latmin + self.latmax) / 2.,
lon_0=(self.lonmin + self.lonmax) / 2.)
# basemap.fillcontinents(color=None, lake_color='aqua',zorder=1)
basemap.drawmapboundary(zorder=2) # fill_color='darkblue')
basemap.shadedrelief(zorder=3)
basemap.drawcountries(zorder=4)
basemap.drawstates(zorder=5)
basemap.drawcoastlines(zorder=6)
# labels = [left,right,top,bottom]
parallels = np.arange(-90, 90, 5.)
parallels_small = np.arange(-90, 90, 2.5)
basemap.drawparallels(parallels_small, linewidth=0.5, zorder=7)
basemap.drawparallels(parallels, zorder=7)
meridians = np.arange(-180, 180, 5.)
meridians_small = np.arange(-180, 180, 2.5)
basemap.drawmeridians(meridians_small, linewidth=0.5, zorder=7)
basemap.drawmeridians(meridians, zorder=7)
self.basemap = basemap
self.figure._tight = True
self.figure.tight_layout()
def init_lat_lon_grid(self, nstep=250):
# create a regular grid to display colormap
lataxis = np.linspace(self.latmin, self.latmax, nstep)
lonaxis = np.linspace(self.lonmin, self.lonmax, nstep)
self.longrid, self.latgrid = np.meshgrid(lonaxis, lataxis)
def init_picksgrid(self):
picks, lats, lons = self.get_picks_lat_lon()
try:
self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid),
method='linear')
except Exception as e:
self._warn('Could not init picksgrid: {}'.format(e))
def get_st_lat_lon_for_plot(self):
stations = []
latitudes = []
longitudes = []
for st_id, coords in self.stations_dict.items():
stations.append(st_id)
latitudes.append(coords['latitude'])
longitudes.append(coords['longitude'])
return stations, latitudes, longitudes
def get_st_x_y_for_plot(self):
stations = []
xs = []
ys = []
for st_id, coords in self.stations_dict.items():
stations.append(st_id)
xs.append(coords['x'])
ys.append(coords['y'])
return stations, xs, ys
def get_picks_lat_lon(self):
picks = []
latitudes = []
longitudes = []
for st_id, pick in self.picks_rel.items():
picks.append(pick)
latitudes.append(self.stations_dict[st_id]['latitude'])
longitudes.append(self.stations_dict[st_id]['longitude'])
return picks, latitudes, longitudes
def draw_contour_filled(self, nlevel='50'):
levels = np.linspace(self.get_min_from_picks(), self.get_max_from_picks(), nlevel)
self.contourf = self.basemap.contourf(self.longrid, self.latgrid, self.picksgrid_active,
levels, latlon=True, zorder=9, alpha=0.5)
def scatter_all_stations(self):
stations, lats, lons = self.get_st_lat_lon_for_plot()
self.sc = self.basemap.scatter(lons, lats, s=50, facecolor='none', latlon=True,
zorder=10, picker=True, edgecolor='m', label='Not Picked')
self.cid = self.canvas.mpl_connect('pick_event', self.onpick)
self._station_onpick_ids = stations
if self.eventLoc:
lats, lons = self.eventLoc
self.sc_event = self.basemap.scatter(lons, lats, s=100, facecolor='red',
latlon=True, zorder=11, label='Event (might be outside map region)')
def scatter_picked_stations(self):
picks, lats, lons = self.get_picks_lat_lon()
# workaround because of an issue with latlon transformation of arrays with len <3
if len(lons) <= 2 and len(lats) <= 2:
self.sc_picked = self.basemap.scatter(lons[0], lats[0], s=50, facecolor='white',
c=picks[0], latlon=True, zorder=11, label='Picked')
if len(lons) == 2 and len(lats) == 2:
self.sc_picked = self.basemap.scatter(lons[1], lats[1], s=50, facecolor='white',
c=picks[1], latlon=True, zorder=11)
else:
self.sc_picked = self.basemap.scatter(lons, lats, s=50, facecolor='white',
c=picks, latlon=True, zorder=11, label='Picked')
def annotate_ax(self):
self.annotations = []
stations, xs, ys = self.get_st_x_y_for_plot()
for st, x, y in zip(stations, xs, ys):
self.annotations.append(self.main_ax.annotate(' %s' % st, xy=(x, y),
fontsize='x-small', color='white', zorder=12))
self.legend = self.main_ax.legend(loc=1)
self.legend.get_frame().set_facecolor((1, 1, 1, 0.75))
def add_cbar(self, label):
self.cbax_bg = inset_axes(self.main_ax, width="6%", height="75%", loc=5)
cbax = inset_axes(self.main_ax, width='2%', height='70%', loc=5)
cbar = self.main_ax.figure.colorbar(self.sc_picked, cax = cbax)
cbar.set_label(label)
cbax.yaxis.tick_left()
cbax.yaxis.set_label_position('left')
for spine in self.cbax_bg.spines.values():
spine.set_visible(False)
self.cbax_bg.yaxis.set_ticks([])
self.cbax_bg.xaxis.set_ticks([])
self.cbax_bg.patch.set_facecolor((1, 1, 1, 0.75))
return cbar
def refresh_drawings(self, picks=None, autopicks=None):
self.picks_dict = picks
self.autopicks_dict = autopicks
self._refresh_drawings()
def _refresh_drawings(self):
self.remove_drawings()
self.draw_everything()
def draw_everything(self):
if self.picks_dict or self.autopicks_dict:
self.init_picks()
if len(self.picks) >= 3:
self.init_picksgrid()
self.draw_contour_filled()
self.scatter_all_stations()
if self.picks_dict or self.autopicks_dict:
self.scatter_picked_stations()
self.cbar = self.add_cbar(label='Time relative to first onset ({}) [s]'.format(self._earliest_picktime))
self.comboBox_phase.setEnabled(True)
else:
self.comboBox_phase.setEnabled(False)
self.annotate_ax()
self.canvas.draw()
def remove_drawings(self):
if hasattr(self, 'cbar'):
self.cbar.remove()
self.cbax_bg.remove()
del (self.cbar, self.cbax_bg)
if hasattr(self, 'sc_picked'):
self.sc_picked.remove()
del (self.sc_picked)
if hasattr(self, 'sc_event'):
self.sc_event.remove()
del (self.sc_event)
if hasattr(self, 'contourf'):
self.remove_contourf()
del (self.contourf)
if hasattr(self, 'cid'):
self.canvas.mpl_disconnect(self.cid)
del (self.cid)
try:
self.sc.remove()
except Exception as e:
print('Warning: could not remove station scatter plot.\nReason: {}'.format(e))
try:
self.legend.remove()
except Exception as e:
print('Warning: could not remove legend. Reason: {}'.format(e))
self.canvas.draw()
def remove_contourf(self):
for item in self.contourf.collections:
item.remove()
def remove_annotations(self):
for annotation in self.annotations:
annotation.remove()
def zoom(self, event):
map = self.basemap
xlim = map.ax.get_xlim()
ylim = map.ax.get_ylim()
x, y = event.xdata, event.ydata
zoom = {'up': 1. / 2.,
'down': 2.}
if not event.xdata or not event.ydata:
return
if event.button in zoom:
factor = zoom[event.button]
xdiff = (xlim[1] - xlim[0]) * factor
xl = x - 0.5 * xdiff
xr = x + 0.5 * xdiff
ydiff = (ylim[1] - ylim[0]) * factor
yb = y - 0.5 * ydiff
yt = y + 0.5 * ydiff
if xl < map.xmin or yb < map.ymin or xr > map.xmax or yt > map.ymax:
xl, xr = map.xmin, map.xmax
yb, yt = map.ymin, map.ymax
map.ax.set_xlim(xl, xr)
map.ax.set_ylim(yb, yt)
map.ax.figure.canvas.draw()
def _warn(self, message):
self.qmb = QtGui.QMessageBox(QtGui.QMessageBox.Icon.Warning,
'Warning', message)
self.qmb.show()

View File

@@ -12,6 +12,186 @@ from pylot.core.util.utils import key_for_set_value, find_in_list, \
remove_underscores, gen_Pool
class Metadata(object):
def __init__(self, inventory=None):
self.inventories = []
if os.path.isdir(inventory):
self.add_inventory(inventory)
if os.path.isfile(inventory):
self.add_inventory_file(inventory)
self.seed_ids = {}
self.inventory_files = {}
def add_inventory(self, path_to_inventory):
# add paths to list of inventories
assert (os.path.isdir(path_to_inventory)), '{} is no directory'.format(path_to_inventory)
if not path_to_inventory in self.inventories:
self.inventories.append(path_to_inventory)
def add_inventory_file(self, path_to_inventory_file):
assert (os.path.isfile(path_to_inventory_file)), '{} is no directory'.format(path_to_inventory_file)
self.add_inventory(os.path.split(path_to_inventory_file)[0])
if not path_to_inventory_file in self.inventory_files.keys():
self.read_single_file(path_to_inventory_file)
def remove_inventory(self, path_to_inventory):
if not path_to_inventory in self.inventories:
print('Path {} not in inventories list.'.format(path_to_inventory))
return
self.inventories.remove(path_to_inventory)
def get_metadata(self, seed_id):
# get metadata for a specific seed_id, if not already read, try to read from inventories
if not seed_id in self.seed_ids.keys():
self._read_inventory_data(seed_id)
# if seed id is not found read all inventories and try to find it there
if not seed_id in self.seed_ids.keys():
print('No data found for seed id {}. Trying to find it in all known inventories...'.format(seed_id))
self.read_all()
for inv_fname, metadata in self.inventory_files.items():
# use get_coordinates to check for seed_id
if metadata['data'].get_coordinates(seed_id):
self.seed_ids[seed_id] = inv_fname
return metadata
print('Could not find metadata for station {}'.format(seed_id))
return None
fname = self.seed_ids[seed_id]
return self.inventory_files[fname]
def read_all(self):
for inventory in self.inventories:
for inv_fname in os.listdir(inventory):
inv_fname = os.path.join(inventory, inv_fname)
if not self.read_single_file(inv_fname):
continue
def read_single_file(self, inv_fname):
if not inv_fname in self.inventory_files.keys():
pass
else:
if not self.inventory_files[inv_fname]:
pass
else:
return
try:
invtype, robj = self._read_metadata_file(inv_fname)
if robj == None:
return
except Exception as e:
print('Could not read file {}'.format(inv_fname))
return
self.inventory_files[inv_fname] = {'invtype': invtype,
'data': robj}
return True
def get_coordinates(self, seed_id):
metadata = self.get_metadata(seed_id)
return metadata['data'].get_coordinates(seed_id)
def get_paz(self, seed_id, time=None):
metadata = self.get_metadata(seed_id)
if metadata['invtype'] in ['dless', 'dseed']:
return metadata['data'].get_paz(seed_id)
elif metadata['invtype'] in ['resp', 'xml']:
if not time:
print('Time needed to extract metadata from station inventory.')
return None
resp = metadata['data'].get_response(seed_id, time)
return resp.get_paz(seed_id)
def _read_inventory_data(self, seed_id=None):
for inventory in self.inventories:
if self._read_metadata_iterator(path_to_inventory=inventory, station_seed_id=seed_id):
return
def _read_metadata_iterator(self, path_to_inventory, station_seed_id):
'''
search for metadata for a specific station iteratively
'''
station, network, location, channel = station_seed_id.split('.')
fnames = glob.glob(os.path.join(path_to_inventory, '*' + station_seed_id + '*'))
if not fnames:
# search for station name in filename
fnames = glob.glob(os.path.join(path_to_inventory, '*' + station + '*'))
if not fnames:
# search for network name in filename
fnames = glob.glob(os.path.join(path_to_inventory, '*' + network + '*'))
if not fnames:
print('Could not find filenames matching station name, network name or seed id')
return
for fname in fnames:
if fname in self.inventory_files.keys():
if self.inventory_files[fname]:
# file already read
continue
invtype, robj = self._read_metadata_file(os.path.join(path_to_inventory, fname))
try:
robj.get_coordinates(station_seed_id)
self.inventory_files[fname] = {'invtype': invtype,
'data': robj}
if station_seed_id in self.seed_ids.keys():
print('WARNING: Overwriting metadata for station {}'.format(station_seed_id))
self.seed_ids[station_seed_id] = fname
return True
except Exception as e:
continue
print('Could not find metadata for station_seed_id {} in path {}'.format(station_seed_id, path_to_inventory))
def _read_metadata_file(self, path_to_inventory_filename):
'''
function reading metadata files (either dataless seed, xml or resp)
:param path_to_inventory_filename:
:return: file type/ending, inventory object (Parser or Inventory)
'''
# functions used to read metadata for different file endings (or file types)
read_functions = {'dless': self._read_dless,
'dseed': self._read_dless,
'xml': self._read_inventory_file,
'resp': self._read_inventory_file}
file_ending = path_to_inventory_filename.split('.')[-1]
if file_ending in read_functions.keys():
robj, exc = read_functions[file_ending](path_to_inventory_filename)
if exc is not None:
raise exc
return file_ending, robj
# in case file endings did not match the above keys, try and error
for file_type in ['dless', 'xml']:
robj, exc = read_functions[file_type](path_to_inventory_filename)
if exc is None:
return file_type, robj
return None, None
def _read_dless(self, path_to_inventory):
exc = None
try:
parser = Parser(path_to_inventory)
except Exception as exc:
parser = None
return parser, exc
def _read_inventory_file(self, path_to_inventory):
exc = None
try:
inv = read_inventory(path_to_inventory)
except Exception as exc:
inv = None
return inv, exc
def time_from_header(header):
"""
Function takes in the second line from a .gse file and takes out the date and time from that line.
@@ -196,6 +376,33 @@ def read_metadata(path_to_inventory):
return invtype, robj
# idea to optimize read_metadata
# def read_metadata_new(path_to_inventory):
# metadata_objects = []
# # read multiple files from directory
# if os.path.isdir(path_to_inventory):
# fnames = os.listdir(path_to_inventory)
# # read single file
# elif os.path.isfile(path_to_inventory):
# fnames = [path_to_inventory]
# else:
# print("Neither dataless-SEED file, inventory-xml file nor "
# "RESP-file found!")
# print("!!WRONG CALCULATION OF SOURCE PARAMETERS!!")
# fnames = []
#
# for fname in fnames:
# path_to_inventory_filename = os.path.join(path_to_inventory, fname)
# try:
# ftype, robj = read_metadata_file(path_to_inventory_filename)
# metadata_objects.append((ftype, robj))
# except Exception as e:
# print('Could not read metadata file {} '
# 'because of the following Exception: {}'.format(path_to_inventory_filename, e))
# return metadata_objects
def restitute_trace(input_tuple):
tr, invtype, inobj, unit, force = input_tuple

View File

@@ -7,6 +7,7 @@ from obspy import UTCDateTime
from obspy.core.event import Event as ObsPyEvent
from obspy.core.event import Origin, ResourceIdentifier
from pylot.core.io.phases import picks_from_picksdict
from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
class Event(ObsPyEvent):
@@ -33,6 +34,7 @@ class Event(ObsPyEvent):
self._testEvent = False
self._refEvent = False
self.get_notes()
self.get_obspy_event_info()
def get_notes_path(self):
"""
@@ -43,6 +45,18 @@ class Event(ObsPyEvent):
notesfile = os.path.join(self.path, 'notes.txt')
return notesfile
def get_obspy_event_info(self):
infile_pickle = os.path.join(self.path, 'info/event.pkl')
if not os.path.isfile(infile_pickle):
return
try:
event_dmt = qml_from_obspyDMT(infile_pickle)
except Exception as e:
print('Could not get obspy event info: {}'.format(e))
return
self.magnitudes = event_dmt.magnitudes
self.origins = event_dmt.origins
def get_notes(self):
"""
set self.note attribute to content of notes file

View File

@@ -1,377 +0,0 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import obspy
from PySide import QtGui
from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
from mpl_toolkits.basemap import Basemap
from pylot.core.util.widgets import PickDlg
from scipy.interpolate import griddata
plt.interactive(False)
class map_projection(QtGui.QWidget):
def __init__(self, parent, figure=None):
'''
:param: picked, can be False, auto, manual
:value: str
'''
QtGui.QWidget.__init__(self)
self._parent = parent
self.metadata = parent.metadata
self.parser = parent.metadata[1]
self.picks = None
self.picks_dict = None
self.eventLoc = None
self.figure = figure
self.init_graphics()
self.init_stations()
self.init_basemap(resolution='l')
self.init_map()
self._style = parent._style
# self.show()
def init_map(self):
self.init_lat_lon_dimensions()
self.init_lat_lon_grid()
self.init_x_y_dimensions()
self.connectSignals()
self.draw_everything()
def onpick(self, event):
ind = event.ind
button = event.mouseevent.button
if ind == [] or not button == 1:
return
data = self._parent.get_data().getWFData()
for index in ind:
station = str(self.station_names[index].split('.')[-1])
try:
pickDlg = PickDlg(self, parameter=self._parent._inputs,
data=data.select(station=station),
station=station,
picks=self._parent.get_current_event().getPick(station),
autopicks=self._parent.get_current_event().getAutopick(station))
except Exception as e:
message = 'Could not generate Plot for station {st}.\n {er}'.format(st=station, er=e)
self._warn(message)
print(message, e)
return
pyl_mw = self._parent
try:
if pickDlg.exec_():
pyl_mw.setDirty(True)
pyl_mw.update_status('picks accepted ({0})'.format(station))
replot = pyl_mw.get_current_event().setPick(station, pickDlg.getPicks())
self._refresh_drawings()
if replot:
pyl_mw.plotWaveformData()
pyl_mw.drawPicks()
pyl_mw.draw()
else:
pyl_mw.drawPicks(station)
pyl_mw.draw()
else:
pyl_mw.update_status('picks discarded ({0})'.format(station))
except Exception as e:
message = 'Could not save picks for station {st}.\n{er}'.format(st=station, er=e)
self._warn(message)
print(message, e)
def connectSignals(self):
self.comboBox_phase.currentIndexChanged.connect(self._refresh_drawings)
self.zoom_id = self.basemap.ax.figure.canvas.mpl_connect('scroll_event', self.zoom)
def init_graphics(self):
if not self.figure:
if not hasattr(self._parent, 'am_figure'):
self.figure = plt.figure()
self.toolbar = NavigationToolbar(self.figure.canvas, self)
else:
self.figure = self._parent.am_figure
self.toolbar = self._parent.am_toolbar
self.main_ax = self.figure.add_subplot(111)
self.canvas = self.figure.canvas
self.main_box = QtGui.QVBoxLayout()
self.setLayout(self.main_box)
self.top_row = QtGui.QHBoxLayout()
self.main_box.addLayout(self.top_row)
self.comboBox_phase = QtGui.QComboBox()
self.comboBox_phase.insertItem(0, 'P')
self.comboBox_phase.insertItem(1, 'S')
self.comboBox_am = QtGui.QComboBox()
self.comboBox_am.insertItem(0, 'auto')
self.comboBox_am.insertItem(1, 'manual')
self.top_row.addWidget(QtGui.QLabel('Select a phase: '))
self.top_row.addWidget(self.comboBox_phase)
self.top_row.setStretch(1, 1) # set stretch of item 1 to 1
self.main_box.addWidget(self.canvas)
self.main_box.addWidget(self.toolbar)
def init_stations(self):
def get_station_names_lat_lon(parser):
station_names = []
lat = []
lon = []
for station in parser.stations:
station_name = station[0].station_call_letters
network = station[0].network_code
if not station_name in station_names:
station_names.append(network + '.' + station_name)
lat.append(station[0].latitude)
lon.append(station[0].longitude)
return station_names, lat, lon
station_names, lat, lon = get_station_names_lat_lon(self.parser)
self.station_names = station_names
self.lat = lat
self.lon = lon
def init_picks(self):
phase = self.comboBox_phase.currentText()
def get_picks(station_names):
picks = []
for station in station_names:
try:
station = station.split('.')[-1]
picks.append(self.picks_dict[station][phase]['mpp'])
except:
picks.append(np.nan)
return picks
def get_picks_rel(picks):
picks_rel = []
picks_utc = []
for pick in picks:
if type(pick) is obspy.core.utcdatetime.UTCDateTime:
picks_utc.append(pick)
minp = min(picks_utc)
for pick in picks:
if type(pick) is obspy.core.utcdatetime.UTCDateTime:
pick -= minp
picks_rel.append(pick)
return picks_rel
self.picks = get_picks(self.station_names)
self.picks_rel = get_picks_rel(self.picks)
def init_picks_active(self):
def remove_nan_picks(picks):
picks_no_nan = []
for pick in picks:
if not np.isnan(pick):
picks_no_nan.append(pick)
return picks_no_nan
self.picks_no_nan = remove_nan_picks(self.picks_rel)
def init_stations_active(self):
def remove_nan_lat_lon(picks, lat, lon):
lat_no_nan = []
lon_no_nan = []
for index, pick in enumerate(picks):
if not np.isnan(pick):
lat_no_nan.append(lat[index])
lon_no_nan.append(lon[index])
return lat_no_nan, lon_no_nan
self.lat_no_nan, self.lon_no_nan = remove_nan_lat_lon(self.picks_rel, self.lat, self.lon)
def init_lat_lon_dimensions(self):
def get_lon_lat_dim(lon, lat):
londim = max(lon) - min(lon)
latdim = max(lat) - min(lat)
return londim, latdim
self.londim, self.latdim = get_lon_lat_dim(self.lon, self.lat)
def init_x_y_dimensions(self):
def get_x_y_dim(x, y):
xdim = max(x) - min(x)
ydim = max(y) - min(y)
return xdim, ydim
self.x, self.y = self.basemap(self.lon, self.lat)
self.xdim, self.ydim = get_x_y_dim(self.x, self.y)
def init_basemap(self, resolution='l'):
# basemap = Basemap(projection=projection, resolution = resolution, ax=self.main_ax)
basemap = Basemap(projection='lcc', resolution=resolution, ax=self.main_ax,
width=5e6, height=2e6,
lat_0=(min(self.lat) + max(self.lat)) / 2.,
lon_0=(min(self.lon) + max(self.lon)) / 2.)
# basemap.fillcontinents(color=None, lake_color='aqua',zorder=1)
basemap.drawmapboundary(zorder=2) # fill_color='darkblue')
basemap.shadedrelief(zorder=3)
basemap.drawcountries(zorder=4)
basemap.drawstates(zorder=5)
basemap.drawcoastlines(zorder=6)
self.basemap = basemap
self.figure.tight_layout()
def init_lat_lon_grid(self):
def get_lat_lon_axis(lat, lon):
steplat = (max(lat) - min(lat)) / 250
steplon = (max(lon) - min(lon)) / 250
lataxis = np.arange(min(lat), max(lat), steplat)
lonaxis = np.arange(min(lon), max(lon), steplon)
return lataxis, lonaxis
def get_lat_lon_grid(lataxis, lonaxis):
longrid, latgrid = np.meshgrid(lonaxis, lataxis)
return latgrid, longrid
self.lataxis, self.lonaxis = get_lat_lon_axis(self.lat, self.lon)
self.latgrid, self.longrid = get_lat_lon_grid(self.lataxis, self.lonaxis)
def init_picksgrid(self):
self.picksgrid_no_nan = griddata((self.lat_no_nan, self.lon_no_nan),
self.picks_no_nan, (self.latgrid, self.longrid),
method='linear') ##################
def draw_contour_filled(self, nlevel='50'):
levels = np.linspace(min(self.picks_no_nan), max(self.picks_no_nan), nlevel)
self.contourf = self.basemap.contourf(self.longrid, self.latgrid, self.picksgrid_no_nan,
levels, latlon=True, zorder=9, alpha=0.5)
def scatter_all_stations(self):
self.sc = self.basemap.scatter(self.lon, self.lat, s=50, facecolor='none', latlon=True,
zorder=10, picker=True, edgecolor='m', label='Not Picked')
self.cid = self.canvas.mpl_connect('pick_event', self.onpick)
if self.eventLoc:
lat, lon = self.eventLoc
self.sc_event = self.basemap.scatter(lon, lat, s=100, facecolor='red',
latlon=True, zorder=11, label='Event (might be outside map region)')
def scatter_picked_stations(self):
lon = self.lon_no_nan
lat = self.lat_no_nan
# workaround because of an issue with latlon transformation of arrays with len <3
if len(lon) <= 2 and len(lat) <= 2:
self.sc_picked = self.basemap.scatter(lon[0], lat[0], s=50, facecolor='white',
c=self.picks_no_nan[0], latlon=True, zorder=11, label='Picked')
if len(lon) == 2 and len(lat) == 2:
self.sc_picked = self.basemap.scatter(lon[1], lat[1], s=50, facecolor='white',
c=self.picks_no_nan[1], latlon=True, zorder=11)
else:
self.sc_picked = self.basemap.scatter(lon, lat, s=50, facecolor='white',
c=self.picks_no_nan, latlon=True, zorder=11, label='Picked')
def annotate_ax(self):
self.annotations = []
for index, name in enumerate(self.station_names):
self.annotations.append(self.main_ax.annotate(' %s' % name, xy=(self.x[index], self.y[index]),
fontsize='x-small', color='white', zorder=12))
self.legend = self.main_ax.legend(loc=1)
def add_cbar(self, label):
cbar = self.main_ax.figure.colorbar(self.sc_picked, fraction=0.025)
cbar.set_label(label)
return cbar
def refresh_drawings(self, picks=None):
self.picks_dict = picks
self._refresh_drawings()
def _refresh_drawings(self):
self.remove_drawings()
self.draw_everything()
def draw_everything(self):
if self.picks_dict:
self.init_picks()
self.init_picks_active()
self.init_stations_active()
if len(self.picks_no_nan) >= 3:
self.init_picksgrid()
self.draw_contour_filled()
self.scatter_all_stations()
if self.picks_dict:
self.scatter_picked_stations()
self.cbar = self.add_cbar(label='Time relative to first onset [s]')
self.comboBox_phase.setEnabled(True)
else:
self.comboBox_phase.setEnabled(False)
self.annotate_ax()
self.canvas.draw()
def remove_drawings(self):
if hasattr(self, 'sc_picked'):
self.sc_picked.remove()
del (self.sc_picked)
if hasattr(self, 'sc_event'):
self.sc_event.remove()
del (self.sc_event)
if hasattr(self, 'cbar'):
self.cbar.remove()
del (self.cbar)
if hasattr(self, 'contourf'):
self.remove_contourf()
del (self.contourf)
if hasattr(self, 'cid'):
self.canvas.mpl_disconnect(self.cid)
del (self.cid)
try:
self.sc.remove()
except Exception as e:
print('Warning: could not remove station scatter plot.\nReason: {}'.format(e))
try:
self.legend.remove()
except Exception as e:
print('Warning: could not remove legend. Reason: {}'.format(e))
self.canvas.draw()
def remove_contourf(self):
for item in self.contourf.collections:
item.remove()
def remove_annotations(self):
for annotation in self.annotations:
annotation.remove()
def zoom(self, event):
map = self.basemap
xlim = map.ax.get_xlim()
ylim = map.ax.get_ylim()
x, y = event.xdata, event.ydata
zoom = {'up': 1. / 2.,
'down': 2.}
if not event.xdata or not event.ydata:
return
if event.button in zoom:
factor = zoom[event.button]
xdiff = (xlim[1] - xlim[0]) * factor
xl = x - 0.5 * xdiff
xr = x + 0.5 * xdiff
ydiff = (ylim[1] - ylim[0]) * factor
yb = y - 0.5 * ydiff
yt = y + 0.5 * ydiff
if xl < map.xmin or yb < map.ymin or xr > map.xmax or yt > map.ymax:
xl, xr = map.xmin, map.xmax
yb, yt = map.ymin, map.ymax
map.ax.set_xlim(xl, xr)
map.ax.set_ylim(yb, yt)
map.ax.figure.canvas.draw()
def _warn(self, message):
self.qmb = QtGui.QMessageBox(QtGui.QMessageBox.Icon.Warning,
'Warning', message)
self.qmb.show()

View File

@@ -25,4 +25,20 @@ def check_obspydmt_eventfolder(folder):
except Exception as e:
return False, e
check_obspydmt_eventfolder('20110311_054623.a')
def qml_from_obspyDMT(path):
import pickle
from obspy.core.event import Event, Magnitude, Origin
if not os.path.exists(path):
return IOError('Could not find Event at {}'.format(path))
infile = open(path, 'rb')
event_dmt = pickle.load(infile)
ev = Event(resource_id=event_dmt['event_id'])
origin = Origin(resource_id=event_dmt['origin_id'], time=event_dmt['datetime'], longitude=event_dmt['longitude'],
latitude=event_dmt['latitude'], depth=event_dmt['depth'])
mag = Magnitude(mag=event_dmt['magnitude'], magnitude_type=event_dmt['magnitude_type'],
origin_id=event_dmt['origin_id'])
ev.magnitudes.append(mag)
ev.origins.append(origin)
return ev

View File

@@ -16,6 +16,9 @@ import time
import numpy as np
import matplotlib
matplotlib.use('QT4Agg')
from matplotlib.figure import Figure
try:
@@ -26,6 +29,7 @@ from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT
from matplotlib.widgets import MultiCursor
from matplotlib.tight_layout import get_renderer, get_subplotspec_list, get_tight_layout_figure
from scipy.signal import argrelmin, argrelmax
from obspy import read
from PySide import QtCore, QtGui
from PySide.QtGui import QAction, QApplication, QCheckBox, QComboBox, \
@@ -447,23 +451,25 @@ class WaveformWidgetPG(QtGui.QWidget):
self.orig_parent = parent
# attribute plotdict is a dictionary connecting position and a name
self.plotdict = dict()
# init labels
self.xlabel = None
self.ylabel = None
self.title = None
# create plot
self.main_layout = QtGui.QVBoxLayout()
self.label_layout = QtGui.QHBoxLayout()
self.status_label = QtGui.QLabel()
self.perm_label = QtGui.QLabel()
self.setLayout(self.main_layout)
self.add_labels()
self.connect_signals()
self.plotWidget = self.pg.PlotWidget(self.parent(), title=title)
self.main_layout.addWidget(self.plotWidget)
self.main_layout.addLayout(self.label_layout)
self.label_layout.addWidget(self.status_label)
self.label_layout.addWidget(self.perm_label)
self.init_labels()
self.activateObspyDMToptions(False)
self.plotWidget.showGrid(x=False, y=True, alpha=0.3)
self.plotWidget.hideAxis('bottom')
self.plotWidget.hideAxis('left')
self.wfstart, self.wfend = 0, 0
self.pen_multicursor = self.pg.mkPen(self.parent()._style['multicursor']['rgba'])
self.pen_linecolor = self.pg.mkPen(self.parent()._style['linecolor']['rgba'])
self.pen_linecolor_highlight = self.pg.mkPen((255, 100, 100, 255))
self.pen_linecolor_syn = self.pg.mkPen((100, 0, 255, 255))
self.reinitMoveProxy()
self._proxy = self.pg.SignalProxy(self.plotWidget.scene().sigMouseMoved, rateLimit=60, slot=self.mouseMoved)
@@ -489,12 +495,59 @@ class WaveformWidgetPG(QtGui.QWidget):
self.vLine.setPos(mousePoint.x())
self.hLine.setPos(mousePoint.y())
def connect_signals(self):
self.qcombo_processed.activated.connect(self.parent().newWF)
self.syn_checkbox.clicked.connect(self.parent().newWF)
def init_labels(self):
self.label_layout.addWidget(self.status_label)
for label in self.perm_labels:
self.label_layout.addWidget(label)
self.label_layout.addWidget(self.syn_checkbox)
self.label_layout.addWidget(self.qcombo_processed)
self.syn_checkbox.setLayoutDirection(Qt.RightToLeft)
self.label_layout.setStretch(0, 4)
self.label_layout.setStretch(1, 0)
self.label_layout.setStretch(2, 0)
self.label_layout.setStretch(3, 0)
self.label_layout.setStretch(4, 3)
self.label_layout.setStretch(5, 1)
def add_labels(self):
self.status_label = QtGui.QLabel()
self.perm_labels = []
for index in range(3):
label = QtGui.QLabel()
self.perm_labels.append(label)
self.qcombo_processed = QtGui.QComboBox()
self.syn_checkbox = QtGui.QCheckBox('synthetics')
self.addQCboxItem('processed', 'green')
self.addQCboxItem('raw', 'black')
#self.perm_qcbox_right.setAlignment(2)
self.setLayout(self.main_layout)
def getPlotDict(self):
return self.plotdict
def setPermText(self, text=None, color='black'):
self.perm_label.setText(text)
self.perm_label.setStyleSheet('color: {}'.format(color))
def activateObspyDMToptions(self, activate):
self.syn_checkbox.setVisible(activate)
self.qcombo_processed.setVisible(activate)
def setPermText(self, number, text=None, color='black'):
if not 0 <= number < len(self.perm_labels):
raise ValueError('No label for number {}'.format(number))
self.perm_labels[number].setText(text)
self.perm_labels[number].setStyleSheet('color: {}'.format(color))
def addQCboxItem(self, text=None, color='black'):
item = QtGui.QStandardItem(text)
model = self.qcombo_processed.model()
model.appendRow(item)
item.setForeground(QtGui.QColor('{}'.format(color)))
def setQCboxItem(self, text):
index = self.qcombo_processed.findText(text)
self.qcombo_processed.setCurrentIndex(index)
def setPlotDict(self, key, value):
self.plotdict[key] = value
@@ -529,6 +582,14 @@ class WaveformWidgetPG(QtGui.QWidget):
else:
st_select = wfdata
gaps = st_select.get_gaps()
if gaps:
merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
st_select.merge()
print('Merged the following stations because of gaps:')
for merged_station in merged:
print(merged_station)
# list containing tuples of network, station, channel (for sorting)
nsc = []
for trace in st_select:
@@ -552,6 +613,8 @@ class WaveformWidgetPG(QtGui.QWidget):
st_syn = wfsyn.select(network=network, station=station, channel=channel)
if st_syn:
trace_syn = st_syn[0].copy()
else:
trace_syn = Trace()
if mapping:
comp = channel[-1]
n = compclass.getPlotPosition(str(comp))
@@ -568,9 +631,11 @@ class WaveformWidgetPG(QtGui.QWidget):
time_ax_syn = prepTimeAxis(stime_syn, trace_syn)
if method == 'fast':
trace.data, time_ax = self.minMax(trace, time_ax)
trace.data, time_ax = self.minMax(trace, time_ax)
if trace_syn:
trace_syn.data, time_ax_syn = self.minMax(trace_syn, time_ax_syn)
if time_ax not in [None, []]:
if len(time_ax) > 0:
if not scaleddata:
trace.detrend('constant')
trace.normalize(np.max(np.abs(trace.data)) * 2)
@@ -581,23 +646,23 @@ class WaveformWidgetPG(QtGui.QWidget):
times = np.array([time for index, time in enumerate(time_ax) if not index % nth_sample])
times_syn = np.array([time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else [])
trace.data = np.array([datum + n for index, datum in enumerate(trace.data) if not index % nth_sample])
trace.data_syn = np.array([datum + n for index, datum in enumerate(trace.data_syn)
trace_syn.data = np.array([datum + n for index, datum in enumerate(trace_syn.data)
if not index % nth_sample] if st_syn else [])
plots.append((times, trace.data,
times_syn, trace.data_syn))
times_syn, trace_syn.data))
self.setPlotDict(n, (station, channel, network))
self.xlabel = 'seconds since {0}'.format(self.wfstart)
self.ylabel = ''
self.setXLims([0, self.wfend - self.wfstart])
self.setYLims([0.5, nmax + 0.5])
return plots
return plots, gaps
def minMax(self, trace, time_ax):
'''
create min/max array for fast plotting (approach based on obspy __plot_min_max function)
:returns data, time_ax
'''
npixel = self.width()
npixel = self.orig_parent.width()
ndata = len(trace.data)
pts_per_pixel = ndata/npixel
if pts_per_pixel < 2:
@@ -1013,6 +1078,14 @@ class PylotCanvas(FigureCanvas):
if mapping:
plot_positions = self.calcPlotPositions(st_select, compclass)
gaps = st_select.get_gaps()
if gaps:
merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
st_select.merge()
print('Merged the following stations because of gaps:')
for merged_station in merged:
print(merged_station)
# list containing tuples of network, station, channel and plot position (for sorting)
nsc = []
for plot_pos, trace in enumerate(st_select):
@@ -1245,9 +1318,10 @@ class PickDlg(QDialog):
def __init__(self, parent=None, data=None, station=None, network=None, picks=None,
autopicks=None, rotate=False, parameter=None, embedded=False, metadata=None,
event=None, filteroptions=None, model='iasp91'):
event=None, filteroptions=None, model='iasp91', wftype=None):
super(PickDlg, self).__init__(parent, 1)
self.orig_parent = parent
self.setAttribute(Qt.WA_DeleteOnClose)
# initialize attributes
self.parameter = parameter
@@ -1256,6 +1330,7 @@ class PickDlg(QDialog):
self.network = network
self.rotate = rotate
self.metadata = metadata
self.wftype = wftype
self.pylot_event = event
self.components = 'ZNE'
self.currentPhase = None
@@ -1305,7 +1380,7 @@ class PickDlg(QDialog):
self.cur_ylim = None
# set attribute holding data
if data is None:
if data is None or not data:
try:
data = parent.get_data().getWFData().copy()
self.data = data.select(station=station)
@@ -1328,7 +1403,7 @@ class PickDlg(QDialog):
# fill compare and scale channels
self.compareChannel.addItem('-', None)
self.scaleChannel.addItem('normalized', None)
self.scaleChannel.addItem('individual', None)
for trace in self.getWFData():
channel = trace.stats.channel
@@ -1345,8 +1420,12 @@ class PickDlg(QDialog):
actionS.setChecked(self.getChannelSettingsS(channel))
# plot data
title = self.getStation()
if self.wftype is not None:
title += ' | ({})'.format(self.wftype)
self.multicompfig.plotWFData(wfdata=self.getWFData(),
title=self.getStation())
title=title)
self.multicompfig.setZoomBorders2content()
@@ -1380,7 +1459,6 @@ class PickDlg(QDialog):
self.setWindowTitle('Pickwindow on station: {}'.format(self.getStation()))
self.setWindowState(QtCore.Qt.WindowMaximized)
self.deleteLater()
def setupUi(self):
menuBar = QtGui.QMenuBar(self)
@@ -1526,7 +1604,7 @@ class PickDlg(QDialog):
_dialtoolbar.addWidget(QtGui.QLabel('Compare to channel: '))
_dialtoolbar.addWidget(self.compareChannel)
_dialtoolbar.addSeparator()
_dialtoolbar.addWidget(QtGui.QLabel('Scale by: '))
_dialtoolbar.addWidget(QtGui.QLabel('Scaling: '))
_dialtoolbar.addWidget(self.scaleChannel)
# layout the innermost widget
@@ -1649,7 +1727,7 @@ class PickDlg(QDialog):
self.arrivalsText.append(ax.text(time_rel, ylims[0], arrival.name, color='0.5'))
def drawArrivalsText(self):
return self.drawArrivals(True)
return self.drawArrivals(textOnly=True)
def refreshArrivalsText(self, event=None):
self.removeArrivalsText()
@@ -2514,6 +2592,9 @@ class PickDlg(QDialog):
filtops_str = transformFilteroptions2String(filtoptions)
title += ' | Filteroptions: {}'.format(filtops_str)
if self.wftype is not None:
title += ' | ({})'.format(self.wftype)
plot_additional = bool(self.compareChannel.currentText())
additional_channel = self.compareChannel.currentText()
scale_channel = self.scaleChannel.currentText()
@@ -2601,6 +2682,7 @@ class PickDlg(QDialog):
phase = 'S'
filter = True
self.plotWFData(phase=phase, filter=filter)
self.drawArrivals()
def resetZoom(self):
ax = self.multicompfig.axes[0]
@@ -2864,7 +2946,7 @@ class AutoPickWidget(MultiEventWidget):
def reinitEvents2plot(self):
for eventID, eventDict in self.events2plot.items():
for widget_key, widget in eventDict.items():
widget.setParent(None)
del(widget)
self.events2plot = {}
self.eventbox.clear()
self.refresh_plot_tabs()
@@ -2946,13 +3028,15 @@ class TuneAutopicker(QWidget):
:type: PyLoT Mainwindow
'''
def __init__(self, parent):
def __init__(self, parent, obspy_dmt=False):
QtGui.QWidget.__init__(self, parent, 1)
self._style = parent._style
self.setWindowTitle('PyLoT - Tune Autopicker')
self.parameter = self.parent()._inputs
self.fig_dict = self.parent().fig_dict
self.data = Data()
self.obspy_dmt = obspy_dmt
self.wftype = None
self.pdlg_widget = None
self.pylot_picks = None
self.init_main_layouts()
@@ -3009,8 +3093,45 @@ class TuneAutopicker(QWidget):
self.eventBox.activated.connect(self.fill_tabs)
self.stationBox.activated.connect(self.fill_tabs)
def catch_station_ids(self):
self.station_ids = {}
eventpath = self.get_current_event_fp()
self.wftype = 'processed' if self.obspy_dmt else ''
wf_path = os.path.join(eventpath, self.wftype)
if not os.path.exists(wf_path) and self.obspy_dmt:
self.wftype = 'raw'
wf_path = os.path.join(eventpath, self.wftype)
for filename in os.listdir(wf_path):
filename = os.path.join(eventpath, self.wftype, filename)
try:
st = read(filename, headonly=True)
except Exception as e:
print('Warning: Could not read file {} as a stream object: {}'.format(filename, e))
continue
for trace in st:
network = trace.stats.network
station = trace.stats.station
location = trace.stats.location
station_id = '{}.{}.{}'.format(network, station, location)
if not station_id in self.station_ids:
self.station_ids[station_id] = []
self.station_ids[station_id].append(filename)
def fill_stationbox(self):
fnames = self.parent().getWFFnames_from_eventbox(eventbox=self.eventBox)
self.stationBox.clear()
model = self.stationBox.model()
self.catch_station_ids()
st_ids_list = list(self.station_ids.keys())
st_ids_list.sort()
for station_id in st_ids_list:
item = QtGui.QStandardItem(station_id)
if station_id.split('.')[1] in self.get_current_event().pylot_picks:
item.setBackground(self.parent()._ref_test_colors['ref'])
model.appendRow(item)
def load_wf_data(self):
fnames = self.station_ids[self.get_current_station_id()]
self.data.setWFData(fnames)
wfdat = self.data.getWFData() # all available streams
# remove possible underscores in station names
@@ -3022,21 +3143,6 @@ class TuneAutopicker(QWidget):
wfdat = check4rotated(wfdat, self.parent().metadata, verbosity=0)
# trim station components to same start value
trim_station_components(wfdat, trim_start=True, trim_end=False)
self.stationBox.clear()
stations = []
for trace in self.data.getWFData():
station = trace.stats.station
network = trace.stats.network
ns_tup = (str(network), str(station))
if not ns_tup in stations:
stations.append(ns_tup)
stations.sort()
model = self.stationBox.model()
for network, station in stations:
item = QtGui.QStandardItem(network + '.' + station)
if station in self.get_current_event().pylot_picks:
item.setBackground(self.parent()._ref_test_colors['ref'])
model.appendRow(item)
def init_figure_tabs(self):
self.figure_tabs = QtGui.QTabWidget()
@@ -3093,10 +3199,14 @@ class TuneAutopicker(QWidget):
def get_current_event_autopicks(self, station):
event = self.get_current_event()
if event.pylot_autopicks:
return event.pylot_autopicks[station]
if station in event.pylot_autopicks:
return event.pylot_autopicks[station]
def get_current_station(self):
return str(self.stationBox.currentText()).split('.')[-1]
return str(self.stationBox.currentText()).split('.')[1]
def get_current_station_id(self):
return str(self.stationBox.currentText())
def gen_tab_widget(self, name, canvas):
widget = QtGui.QWidget()
@@ -3111,17 +3221,19 @@ class TuneAutopicker(QWidget):
self.pdlg_widget.setParent(None)
self.pdlg_widget = None
return
self.load_wf_data()
station = self.get_current_station()
data = self.data.getWFData()
wfdata = self.data.getWFData()
metadata = self.parent().metadata
event = self.get_current_event()
filteroptions = self.parent().filteroptions
self.pickDlg = PickDlg(self.parent(), data=data.select(station=station),
wftype = self.wftype if self.obspy_dmt else ''
self.pickDlg = PickDlg(self.parent(), data=wfdata.select(station=station).copy(),
station=station, parameter=self.parameter,
picks=self.get_current_event_picks(station),
autopicks=self.get_current_event_autopicks(station),
metadata=metadata, event=event, filteroptions=filteroptions,
embedded=True)
embedded=True, wftype=wftype)
self.pickDlg.update_picks.connect(self.picks_from_pickdlg)
self.pickDlg.update_picks.connect(self.fill_eventbox)
self.pickDlg.update_picks.connect(self.fill_stationbox)
@@ -3301,14 +3413,17 @@ class TuneAutopicker(QWidget):
if not station:
self._warn('No station selected')
return
wfpath = self.wftype if self.obspy_dmt else ''
args = {'parameter': self.parameter,
'station': station,
'fnames': 'None',
'eventid': [self.get_current_event_fp()],
'iplot': 2,
'fig_dict': self.fig_dict,
'locflag': 0,
'savexml': False}
'savexml': False,
'obspyDMT_wfpath': wfpath}
event = self.get_current_event()
self.parent().saveData(event, event.path, '.xml')
for key in self.fig_dict.keys():
if not key == 'plot_style':
self.fig_dict[key].clear()
@@ -3364,8 +3479,7 @@ class TuneAutopicker(QWidget):
if hasattr(self, 'pdlg_widget'):
if self.pdlg_widget:
self.pdlg_widget.setParent(None)
# TODO: removing widget by parent deletion raises exception when activating stationbox:
# RuntimeError: Internal C++ object (PylotCanvas) already deleted.
del(self.pdlg_widget)
if hasattr(self, 'overview'):
self.overview.setParent(None)
if hasattr(self, 'p_tabs'):

View File

@@ -179,9 +179,9 @@ background-color:transparent;
QTabWidget::pane{
background-color:rgba(0, 0, 0, 255);
border-style:solid;
border-color:rgba(245, 245, 245, 255);
border-width:1px;
margin: 0px, 0px, 0px, 0px;
padding: 0px;
border-width:0px;
}
QTabWidget::tab{

View File

@@ -178,9 +178,9 @@ background-color:transparent;
QTabWidget::pane{
background-color:rgba(70, 70, 80, 255);
border-style:solid;
border-color:rgba(70, 70, 80, 255);
border-width:1px;
margin: 0px, 0px, 0px, 0px;
padding: 0px;
border-width:0px;
}
QTabWidget::tab{