This commit is contained in:
2022-03-24 13:28:52 +01:00
45 changed files with 3665 additions and 1568 deletions
+17 -17
View File
@@ -223,14 +223,14 @@ class LocalMagnitude(Magnitude):
in 'Z3']
# checking horizontal count and calculating power_sum accordingly
if len(power) == 1:
print ('WARNING: Only one horizontal found for station {0}.'.format(st[0].stats.station))
power_sum = power[0]
print('WARNING: Only one horizontal found for station {0}.'.format(st[0].stats.station))
power_sum = power[0]
elif len(power) == 2:
power_sum = power[0] + power[1]
else:
raise ValueError('Wood-Anderson aomplitude defintion only valid for'
' up to two horizontals: {0} given'.format(len(power)))
power_sum = power[0] + power[1]
else:
raise ValueError('Wood-Anderson aomplitude defintion only valid for'
' up to two horizontals: {0} given'.format(len(power)))
sqH = np.sqrt(power_sum)
# get time array
@@ -325,7 +325,7 @@ class LocalMagnitude(Magnitude):
if self.verbose:
print(
"Local Magnitude for station {0}: ML = {1:3.1f}".format(
station, magnitude.mag))
station, magnitude.mag))
magnitude.origin_id = self.origin_id
magnitude.waveform_id = pick.waveform_id
magnitude.amplitude_id = amplitude.resource_id
@@ -404,7 +404,7 @@ class MomentMagnitude(Magnitude):
if not wf:
continue
try:
scopy = wf.copy()
scopy = wf.copy()
except AssertionError:
print("WARNING: Something's wrong with the data,"
"station {},"
@@ -464,7 +464,7 @@ def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
# additional common parameters for calculating Mo
# average radiation pattern of P waves (Aki & Richards, 1980)
rP = 2 / np.sqrt(15)
rP = 2 / np.sqrt(15)
freesurf = 2.0 # free surface correction, assuming vertical incidence
Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf)
@@ -524,7 +524,7 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
iplot = 2
else:
iplot = 0
# get Q value
Q, A = qp
@@ -594,13 +594,13 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
# fft
fny = freq / 2
#l = len(xdat) / freq
# l = len(xdat) / freq
# number of fft bins after Bath
#n = freq * l
# n = freq * l
# find next power of 2 of data length
m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
N = min(int(np.power(m, 2)), 16384)
#N = int(np.power(m, 2))
# N = int(np.power(m, 2))
y = dt * np.fft.fft(xdat, N)
Y = abs(y[: N / 2])
L = (N - 1) / freq
@@ -643,8 +643,8 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
w0 = np.median([w01, w02])
Fc = np.median([fc1, fc2])
if verbosity:
print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (
w0, Fc))
print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (
w0, Fc))
if iplot >= 1:
f1 = plt.figure()
tLdat = np.arange(0, len(Ldat) * dt, dt)
@@ -733,7 +733,7 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
iplot = 2
else:
iplot = 0
w0 = []
stdw0 = []
fc = []
+25 -24
View File
@@ -1,18 +1,18 @@
#!/usr/bin/env pyth n
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import copy
import os
from PySide2.QtWidgets import QMessageBox
from obspy import read_events
from obspy.core import read, Stream, UTCDateTime
from obspy.core.event import Event as ObsPyEvent
from obspy.io.sac import SacIOError
from PySide.QtGui import QMessageBox
import pylot.core.loc.velest as velest
import pylot.core.loc.focmec as focmec
import pylot.core.loc.hypodd as hypodd
import pylot.core.loc.velest as velest
from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \
picksdict_from_pilot, merge_picks, PylotParameter
from pylot.core.util.errors import FormatError, OverwriteError
@@ -21,13 +21,14 @@ from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
from pylot.core.util.utils import fnConstructor, full_range, check4rotated, \
check4gapsAndMerge, trim_station_components
class Data(object):
"""
Data container with attributes wfdata holding ~obspy.core.stream.
:type parent: PySide.QtGui.QWidget object, optional
:param parent: A PySide.QtGui.QWidget object utilized when
called by a GUI to display a PySide.QtGui.QMessageBox instead of printing
:type parent: PySide2.QtWidgets.QWidget object, optional
:param parent: A PySide2.QtWidgets.QWidget object utilized when
called by a GUI to display a PySide2.QtWidgets.QMessageBox instead of printing
to standard out.
:type evtdata: ~obspy.core.event.Event object, optional
:param evtdata ~obspy.core.event.Event object containing all derived or
@@ -47,10 +48,10 @@ class Data(object):
elif isinstance(evtdata, dict):
evt = readPILOTEvent(**evtdata)
evtdata = evt
elif type(evtdata) in [str, unicode]:
elif isinstance(evtdata, str):
try:
cat = read_events(evtdata)
if len(cat) is not 1:
if len(cat) != 1:
raise ValueError('ambiguous event information for file: '
'{file}'.format(file=evtdata))
evtdata = cat[0]
@@ -99,7 +100,7 @@ class Data(object):
old_pick.phase_hint == new_pick.phase_hint,
old_pick.method_id == new_pick.method_id]
if all(comparison):
del(old_pick)
del (old_pick)
old_picks.append(new_pick)
elif not other.isNew() and self.isNew():
new = other + self
@@ -111,7 +112,7 @@ class Data(object):
return self + other
else:
raise ValueError("both Data objects have differing "
"unique Event identifiers")
"unique Event identifiers")
return self
def getPicksStr(self):
@@ -162,9 +163,8 @@ class Data(object):
def checkEvent(self, event, fcheck, forceOverwrite=False):
"""
Check information in supplied event and own event and replace own
information with supplied information if own information not exiisting
or forced by forceOverwrite
Check information in supplied event and own event and replace with own
information if no other information are given or forced by forceOverwrite
:param event: Event that supplies information for comparison
:type event: pylot.core.util.event.Event
:param fcheck: check and delete existing information
@@ -250,7 +250,7 @@ class Data(object):
for pick in self.get_evt_data().picks:
if picktype in str(pick.method_id.id):
picks.append(pick)
def exportEvent(self, fnout, fnext='.xml', fcheck='auto', upperErrors=None):
"""
Export event to file
@@ -260,7 +260,7 @@ class Data(object):
can be a str or a list of strings of ['manual', 'auto', 'origin', 'magnitude']
"""
from pylot.core.util.defaults import OUTPUTFORMATS
if not type(fcheck) == list:
fcheck = [fcheck]
@@ -293,7 +293,7 @@ class Data(object):
return
self.checkEvent(event, fcheck)
self.setEvtData(event)
self.get_evt_data().write(fnout + fnext, format=evtformat)
# try exporting event
@@ -318,7 +318,7 @@ class Data(object):
del picks_copy[k]
break
lendiff = len(picks) - len(picks_copy)
if lendiff is not 0:
if lendiff != 0:
print("Manual as well as automatic picks available. Prefered the {} manual ones!".format(lendiff))
if upperErrors:
@@ -360,13 +360,13 @@ class Data(object):
nllocfile = open(fnout + fnext)
l = nllocfile.readlines()
# Adding A0/Generic Amplitude to .obs file
#l2 = []
#for li in l:
# l2 = []
# for li in l:
# for amp in evtdata_org.amplitudes:
# if amp.waveform_id.station_code == li[0:5].strip():
# li = li[0:64] + '{:0.2e}'.format(amp.generic_amplitude) + li[73:-1] + '\n'
# l2.append(li)
#l = l2
# l = l2
nllocfile.close()
l.insert(0, header)
nllocfile = open(fnout + fnext, 'w')
@@ -503,7 +503,8 @@ class Data(object):
real_or_syn_data[synthetic] += read(fname, format='GSE2', starttime=self.tstart, endtime=self.tstop)
except Exception as e:
try:
real_or_syn_data[synthetic] += read(fname, format='SEGY', starttime=self.tstart, endtime=self.tstop)
real_or_syn_data[synthetic] += read(fname, format='SEGY', starttime=self.tstart,
endtime=self.tstop)
except Exception as e:
warnmsg += '{0}\n{1}\n'.format(fname, e)
except SacIOError as se:
@@ -794,8 +795,8 @@ class PilotDataStructure(GenericDataStructure):
def __init__(self, **fields):
if not fields:
fields = {'database': '2006.01',
'root': '/data/Egelados/EVENT_DATA/LOCAL'}
fields = {'database': '',
'root': ''}
GenericDataStructure.__init__(self, **fields)
+3 -3
View File
@@ -515,9 +515,9 @@ defaults = {'rootpath': {'type': str,
'namestring': 'TauPy model'},
'taup_phases': {'type': str,
'tooltip': 'Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.',
'value': 'ttall',
'namestring': 'TauPy phases'},
'tooltip': 'Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.',
'value': 'ttall',
'namestring': 'TauPy phases'},
}
settings_main = {
+20 -12
View File
@@ -8,14 +8,18 @@
Edited for use in PyLoT
JG, igem, 01/2022
"""
<<<<<<< HEAD
import os
import argparse
import numpy as np
import matplotlib.pyplot as plt
=======
import glob
>>>>>>> 83ba63a3fdd8ad0a212333928cfcd6f45fcb5baa
from obspy.core.event import read_events
from pyproj import Proj
import glob
"""
Creates an eventlist file summarizing all events found in a certain folder. Only called by pressing UI Button eventlis_xml_action
@@ -24,14 +28,15 @@ Creates an eventlist file summarizing all events found in a certain folder. Only
:param path: Path to root folder where single Event folder are to found
"""
def geteventlistfromxml(path, outpath):
p = Proj(proj='utm', zone=32, ellps='WGS84')
# open eventlist file and write header
evlist = outpath + '/eventlist'
evlistobj = open(evlist, 'w')
evlistobj.write('EventID Date To Lat Lon EAST NORTH Dep Ml NoP NoS RMS errH errZ Gap \n')
evlistobj.write(
'EventID Date To Lat Lon EAST NORTH Dep Ml NoP NoS RMS errH errZ Gap \n')
# data path
dp = path + "/e*/*.xml"
@@ -52,28 +57,31 @@ def geteventlistfromxml(path, outpath):
NoP = []
NoS = []
except IndexError:
print ('Insufficient data found for event (not localised): ' + names.split('/')[-1].split('_')[-1][:-4] + ' Skipping event for eventlist.' )
print('Insufficient data found for event (not localised): ' + names.split('/')[-1].split('_')[-1][
:-4] + ' Skipping event for eventlist.')
continue
for i in range(len(cat.events[0].origins[0].arrivals)):
if cat.events[0].origins[0].arrivals[i].phase == 'P':
NoP.append(cat.events[0].origins[0].arrivals[i].phase)
elif cat.events[0].origins[0].arrivals[i].phase == 'S':
NoS.append(cat.events[0].origins[0].arrivals[i].phase)
#NoP = cat.events[0].origins[0].quality.used_station_count
# NoP = cat.events[0].origins[0].quality.used_station_count
errH = cat.events[0].origins[0].origin_uncertainty.max_horizontal_uncertainty
errZ = cat.events[0].origins[0].depth_errors.uncertainty
Gap = cat.events[0].origins[0].quality.azimuthal_gap
#evID = names.split('/')[6]
# evID = names.split('/')[6]
evID = names.split('/')[-1].split('_')[-1][:-4]
Date = str(st.year) + str('%02d' % st.month) + str('%02d' % st.day)
To = str('%02d' % st.hour) + str('%02d' % st.minute) + str('%02d' % st.second) + \
'.' + str('%06d' % st.microsecond)
'.' + str('%06d' % st.microsecond)
# write into eventlist
evlistobj.write('%s %s %s %9.6f %9.6f %13.6f %13.6f %8.6f %3.1f %d %d NaN %d %d %d\n' %(evID, \
Date, To, Lat, Lon, EAST, NORTH, Dep, Ml, len(NoP), len(NoS), errH, errZ, Gap))
print ('Adding Event ' + names.split('/')[-1].split('_')[-1][:-4] + ' to eventlist')
evlistobj.write('%s %s %s %9.6f %9.6f %13.6f %13.6f %8.6f %3.1f %d %d NaN %d %d %d\n' % (evID, \
Date, To, Lat, Lon,
EAST, NORTH, Dep, Ml,
len(NoP), len(NoS),
errH, errZ, Gap))
print('Adding Event ' + names.split('/')[-1].split('_')[-1][:-4] + ' to eventlist')
print('Eventlist created and saved in: ' + outpath)
evlistobj.close()
-139
View File
@@ -1,139 +0,0 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Script to get onset uncertainties from Quakeml.xml files created by PyLoT.
Uncertainties are tranformed into quality classes and visualized via histogram if desired.
Ludger Küperkoch, BESTEC GmbH, 07/2017
rev.: Ludger Küperkoch, igem, 10/2020
Edited for usage in PyLoT: Jeldrik Gaal, igem, 01/2022
"""
import argparse
import numpy as np
import matplotlib.pyplot as plt
from obspy.core.event import read_events
import glob
def getQualitiesfromxml(path):
# uncertainties
ErrorsP = [0.02, 0.04, 0.08, 0.16]
ErrorsS = [0.04, 0.08, 0.16, 0.32]
Pw0 = []
Pw1 = []
Pw2 = []
Pw3 = []
Pw4 = []
Sw0 = []
Sw1 = []
Sw2 = []
Sw3 = []
Sw4 = []
# data path
dp = path + '/e*/*.xml'
# list of all available xml-files
xmlnames = glob.glob(dp)
# read all onset weights
for names in xmlnames:
print("Getting onset weights from {}".format(names))
cat = read_events(names)
arrivals = cat.events[0].picks
for Pick in arrivals:
if Pick.phase_hint[0] == 'P':
if Pick.time_errors.uncertainty <= ErrorsP[0]:
Pw0.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsP[0] and \
Pick.time_errors.uncertainty <= ErrorsP[1]:
Pw1.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsP[1] and \
Pick.time_errors.uncertainty <= ErrorsP[2]:
Pw2.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsP[2] and \
Pick.time_errors.uncertainty <= ErrorsP[3]:
Pw3.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsP[3]:
Pw4.append(Pick.time_errors.uncertainty)
else:
pass
elif Pick.phase_hint[0] == 'S':
if Pick.time_errors.uncertainty <= ErrorsS[0]:
Sw0.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsS[0] and \
Pick.time_errors.uncertainty <= ErrorsS[1]:
Sw1.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsS[1] and \
Pick.time_errors.uncertainty <= ErrorsS[2]:
Sw2.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsS[2] and \
Pick.time_errors.uncertainty <= ErrorsS[3]:
Sw3.append(Pick.time_errors.uncertainty)
elif Pick.time_errors.uncertainty > ErrorsS[3]:
Sw4.append(Pick.time_errors.uncertainty)
else:
pass
else:
print("Phase hint not defined for picking!")
pass
# get percentage of weights
numPweights = np.sum([len(Pw0), len(Pw1), len(Pw2), len(Pw3), len(Pw4)])
numSweights = np.sum([len(Sw0), len(Sw1), len(Sw2), len(Sw3), len(Sw4)])
try:
P0perc = 100.0 / numPweights * len(Pw0)
except:
P0perc = 0
try:
P1perc = 100.0 / numPweights * len(Pw1)
except:
P1perc = 0
try:
P2perc = 100.0 / numPweights * len(Pw2)
except:
P2perc = 0
try:
P3perc = 100.0 / numPweights * len(Pw3)
except:
P3perc = 0
try:
P4perc = 100.0 / numPweights * len(Pw4)
except:
P4perc = 0
try:
S0perc = 100.0 / numSweights * len(Sw0)
except:
Soperc = 0
try:
S1perc = 100.0 / numSweights * len(Sw1)
except:
S1perc = 0
try:
S2perc = 100.0 / numSweights * len(Sw2)
except:
S2perc = 0
try:
S3perc = 100.0 / numSweights * len(Sw3)
except:
S3perc = 0
try:
S4perc = 100.0 / numSweights * len(Sw4)
except:
S4perc = 0
weights = ('0', '1', '2', '3', '4')
y_pos = np.arange(len(weights))
width = 0.34
p1, = plt.bar(0 - width, P0perc, width, color='black')
p2, = plt.bar(0, S0perc, width, color='red')
plt.bar(y_pos - width, [P0perc, P1perc, P2perc, P3perc, P4perc], width, color='black')
plt.bar(y_pos, [S0perc, S1perc, S2perc, S3perc, S4perc], width, color='red')
plt.ylabel('%')
plt.xticks(y_pos, weights)
plt.xlim([-0.5, 4.5])
plt.xlabel('Qualities')
plt.title('{0} P-Qualities, {1} S-Qualities'.format(numPweights, numSweights))
plt.legend([p1, p2], ['P-Weights', 'S-Weights'])
plt.show()
+124 -158
View File
@@ -1,13 +1,13 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import pdb
import glob
import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import obspy.core.event as ope
import os
import scipy.io as sio
import warnings
from obspy.core import UTCDateTime
from obspy.core.event import read_events
from obspy.core.util import AttribDict
@@ -17,7 +17,7 @@ from pylot.core.io.location import create_event, \
create_magnitude
from pylot.core.pick.utils import select_for_phase, get_quality_class
from pylot.core.util.utils import getOwner, full_range, four_digits, transformFilterString4Export, \
backtransformFilterString
backtransformFilterString, loopIdentifyPhase, identifyPhase
def add_amplitudes(event, amplitudes):
@@ -232,7 +232,7 @@ def picksdict_from_picks(evt):
for pick in evt.picks:
phase = {}
station = pick.waveform_id.station_code
if pick.waveform_id.channel_code == None:
if pick.waveform_id.channel_code is None:
channel = ''
else:
channel = pick.waveform_id.channel_code
@@ -375,8 +375,7 @@ def picks_from_picksdict(picks, creation_info=None):
def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0):
import glob
# TODO: change root to datapath
db_root = os.path.join(root_dir, db_dir)
evt_list = glob.glob1(db_root, 'e????.???.??')
@@ -391,9 +390,7 @@ def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None
from pylot.core.io.inputs import PylotParameter
from pylot.core.pick.utils import earllatepicker
if fn_param is None:
fn_param = defaults.AUTOMATIC_DEFAULTS
# TODO: change root to datapath
default = PylotParameter(fn_param, verbosity)
@@ -483,7 +480,6 @@ def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None
os.makedirs(out_dir)
fnout_prefix = os.path.join(out_dir, 'PyLoT_{0}.'.format(event_id))
evt.write(fnout_prefix + 'xml', format='QUAKEML')
# evt.write(fnout_prefix + 'cnv', format='VELEST')
def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
@@ -511,14 +507,14 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
:param eventinfo: optional, needed for VELEST-cnv file
and FOCMEC- and HASH-input files
:type eventinfo: `obspy.core.event.Event` object
"""
"""
if fformat == 'NLLoc':
print("Writing phases to %s for NLLoc" % filename)
fid = open("%s" % filename, 'w')
# write header
fid.write('# EQEVENT: %s Label: EQ%s Loc: X 0.00 Y 0.00 Z 10.00 OT 0.00 \n' %
(parameter.get('database'), parameter.get('eventID')))
arrivals = chooseArrivals(arrivals)
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
for key in arrivals:
# P onsets
if arrivals[key].has_key('P'):
@@ -572,20 +568,20 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
sweight = 0 # do not use pick
except KeyError as e:
print(str(e) + '; no weight set during processing')
Ao = arrivals[key]['S']['Ao'] # peak-to-peak amplitude
Ao = arrivals[key]['S']['Ao'] # peak-to-peak amplitude
if Ao == None:
Ao = 0.0
#fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 %d \n' % (key,
# fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 0 0 0 %d \n' % (key,
fid.write('%s ? ? ? S %s %d%02d%02d %02d%02d %7.4f GAU 0 %9.2f 0 0 %d \n' % (key,
fm,
year,
month,
day,
hh,
mm,
ss_ms,
Ao,
sweight))
fm,
year,
month,
day,
hh,
mm,
ss_ms,
Ao,
sweight))
fid.close()
elif fformat == 'HYPO71':
@@ -594,7 +590,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
# write header
fid.write(' %s\n' %
parameter.get('eventID'))
arrivals = chooseArrivals(arrivals)
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
for key in arrivals:
if arrivals[key]['P']['weight'] < 4:
stat = key
@@ -671,7 +667,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
fid = open("%s" % filename, 'w')
# write header
fid.write('%s, event %s \n' % (parameter.get('database'), parameter.get('eventID')))
arrivals = chooseArrivals(arrivals)
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
for key in arrivals:
# P onsets
if arrivals[key].has_key('P') and arrivals[key]['P']['mpp'] is not None:
@@ -769,7 +765,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
arrivals = picksdict_from_picks(evt)
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrival(arrivals)
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
# P onsets
if usedarrivals[key].has_key('P'):
@@ -781,7 +777,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
Ponset = usedarrivals[key]['P']['mpp']
Pweight = usedarrivals[key]['P']['weight']
Prt = Ponset - stime # onset time relative to source time
if n % 6 is not 0:
if n % 6 != 0:
fid.write('%-4sP%d%6.2f' % (stat, Pweight, Prt))
else:
fid.write('%-4sP%d%6.2f\n' % (stat, Pweight, Prt))
@@ -795,7 +791,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
Sonset = usedarrivals[key]['S']['mpp']
Sweight = usedarrivals[key]['S']['weight']
Srt = Ponset - stime # onset time relative to source time
if n % 6 is not 0:
if n % 6 != 0:
fid.write('%-4sS%d%6.2f' % (stat, Sweight, Srt))
else:
fid.write('%-4sS%d%6.2f\n' % (stat, Sweight, Srt))
@@ -815,9 +811,9 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
event = eventinfo['pylot_id']
hddID = event.split('.')[0][1:5]
except:
print ("Error 1111111!")
hddID = "00000"
# write header
print("Error 1111111!")
hddID = "00000"
# write header
fid.write('# %d %d %d %d %d %5.2f %7.4f +%6.4f %7.4f %4.2f 0.1 0.5 %4.2f %s\n' % (
stime.year, stime.month, stime.day, stime.hour, stime.minute, stime.second,
eventsource['latitude'], eventsource['longitude'], eventsource['depth'] / 1000,
@@ -830,7 +826,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
arrivals = picksdict_from_picks(evt)
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrival(arrivals)
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
if usedarrivals[key].has_key('P'):
# P onsets
@@ -881,7 +877,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
arrivals = picksdict_from_picks(evt)
# check for automatic and manual picks
# prefer manual picks
usedarrivals = chooseArrival(arrivals)
usedarrivals = chooseArrivals(arrivals)
for key in usedarrivals:
if usedarrivals[key].has_key('P'):
if usedarrivals[key]['P']['weight'] < 4 and usedarrivals[key]['P']['fm'] is not None:
@@ -962,7 +958,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
erh, erz, eventinfo.magnitudes[0]['mag'],
hashID))
# Prefer Manual Picks over automatic ones if possible
arrivals = chooseArrivals(arrivals)
arrivals = chooseArrivals(arrivals) # MP MP what is chooseArrivals? It is not defined anywhere
# write phase lines
for key in arrivals:
if arrivals[key].has_key('P'):
@@ -1009,7 +1005,8 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
fid1.close()
fid2.close()
def chooseArrival(arrivals):
def chooseArrivals(arrivals):
"""
takes arrivals and returns the manual picks if manual and automatic ones are there
returns automatic picks if only automatic picks are there
@@ -1058,37 +1055,63 @@ def merge_picks(event, picks):
return event
def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
def getQualitiesfromxml(path, errorsP, errorsS, plotflag=1, figure=None, verbosity=0):
"""
Script to get onset uncertainties from Quakeml.xml files created by PyLoT.
Uncertainties are tranformed into quality classes and visualized via histogram if desired.
Ludger Küperkoch, BESTEC GmbH, 07/2017
:param xmlnames: list of xml obspy event files containing picks
:type xmlnames: list
:param ErrorsP: time errors of P waves for the four discrete quality classes
:type ErrorsP:
:param ErrorsS: time errors of S waves for the four discrete quality classes
:type ErrorsS:
:param path: path containing xml files
:type path: str
:param errorsP: time errors of P waves for the four discrete quality classes
:type errorsP:
:param errorsS: time errors of S waves for the four discrete quality classes
:type errorsS:
:param plotflag:
:type plotflag:
:return:
:rtype:
"""
from pylot.core.pick.utils import get_quality_class
from pylot.core.util.utils import loopIdentifyPhase, identifyPhase
def calc_perc(uncertainties, ntotal):
''' simple function that calculates percentage of number of uncertainties (list length)'''
if len(uncertainties) == 0:
return 0
else:
return 100. / ntotal * len(uncertainties)
def calc_weight_perc(psweights, weight_ids):
''' calculate percentages of different weights (pick classes!?) of total number of uncertainties of a phase'''
# count total number of list items for this phase
numWeights = np.sum([len(weight) for weight in psweights.values()])
# iterate over all available weights to return a list with percentages for plotting
plot_list = []
for weight_id in weight_ids:
plot_list.append(calc_perc(psweights[weight_id], numWeights))
return plot_list, numWeights
# get all xmlfiles in path (maybe this should be changed to one xml file for this function, selectable via GUI?)
xmlnames = glob.glob(os.path.join(path, '*.xml'))
if len(xmlnames) == 0:
print(f'No files found in path {path}.')
return False
# first define possible phases here
phases = ['P', 'S']
# define possible weights (0-4)
weight_ids = list(range(5))
# put both error lists in a dictionary with P/S key so that amount of code can be halfed by simply using P/S as key
errors = dict(P=errorsP, S=errorsS)
# create dictionaries for each phase (P/S) with a dictionary of empty list for each weight defined in weights
# tuple above
weights = {}
for phase in phases:
weights[phase] = {weight_id: [] for weight_id in weight_ids}
# read all onset weights
Pw0 = []
Pw1 = []
Pw2 = []
Pw3 = []
Pw4 = []
Sw0 = []
Sw1 = []
Sw2 = []
Sw3 = []
Sw4 = []
for names in xmlnames:
print("Getting onset weights from {}".format(names))
cat = read_events(names)
@@ -1096,117 +1119,60 @@ def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
arrivals = cat.events[0].picks
arrivals_copy = cat_copy.events[0].picks
# Prefere manual picks if qualities are sufficient!
for Pick in arrivals:
if Pick.method_id.id.split('/')[1] == 'manual':
mstation = Pick.waveform_id.station_code
for pick in arrivals:
if pick.method_id.id.split('/')[1] == 'manual':
mstation = pick.waveform_id.station_code
mstation_ext = mstation + '_'
for mpick in arrivals_copy:
phase = identifyPhase(loopIdentifyPhase(Pick.phase_hint))
if phase == 'P':
if ((mpick.waveform_id.station_code == mstation) or
(mpick.waveform_id.station_code == mstation_ext)) and \
(mpick.method_id.split('/')[1] == 'auto') and \
(mpick.time_errors['uncertainty'] <= ErrorsP[3]):
del mpick
break
elif phase == 'S':
if ((mpick.waveform_id.station_code == mstation) or
(mpick.waveform_id.station_code == mstation_ext)) and \
(mpick.method_id.split('/')[1] == 'auto') and \
(mpick.time_errors['uncertainty'] <= ErrorsS[3]):
del mpick
break
phase = identifyPhase(loopIdentifyPhase(pick.phase_hint)) # MP MP catch if this fails?
if ((mpick.waveform_id.station_code == mstation) or
(mpick.waveform_id.station_code == mstation_ext)) and \
(mpick.method_id.id.split('/')[1] == 'auto') and \
(mpick.time_errors['uncertainty'] <= errors[phase][3]):
del mpick
break
lendiff = len(arrivals) - len(arrivals_copy)
if lendiff is not 0:
if lendiff != 0:
print("Found manual as well as automatic picks, prefered the {} manual ones!".format(lendiff))
for Pick in arrivals_copy:
phase = identifyPhase(loopIdentifyPhase(Pick.phase_hint))
if phase == 'P':
Pqual = get_quality_class(Pick.time_errors.uncertainty, ErrorsP)
if Pqual == 0:
Pw0.append(Pick.time_errors.uncertainty)
elif Pqual == 1:
Pw1.append(Pick.time_errors.uncertainty)
elif Pqual == 2:
Pw2.append(Pick.time_errors.uncertainty)
elif Pqual == 3:
Pw3.append(Pick.time_errors.uncertainty)
elif Pqual == 4:
Pw4.append(Pick.time_errors.uncertainty)
elif phase == 'S':
Squal = get_quality_class(Pick.time_errors.uncertainty, ErrorsS)
if Squal == 0:
Sw0.append(Pick.time_errors.uncertainty)
elif Squal == 1:
Sw1.append(Pick.time_errors.uncertainty)
elif Squal == 2:
Sw2.append(Pick.time_errors.uncertainty)
elif Squal == 3:
Sw3.append(Pick.time_errors.uncertainty)
elif Squal == 4:
Sw4.append(Pick.time_errors.uncertainty)
else:
for pick in arrivals_copy:
phase = identifyPhase(loopIdentifyPhase(pick.phase_hint))
uncertainty = pick.time_errors.uncertainty
if not uncertainty:
if verbosity > 0:
print('No uncertainty, pick {} invalid!'.format(pick.method_id.id))
continue
# check P/S phase
if phase not in phases:
print("Phase hint not defined for picking!")
pass
continue
qual = get_quality_class(uncertainty, errors[phase])
weights[phase][qual].append(uncertainty)
if plotflag == 0:
Punc = [Pw0, Pw1, Pw2, Pw3, Pw4]
Sunc = [Sw0, Sw1, Sw2, Sw3, Sw4]
return Punc, Sunc
p_unc = [weights['P'][weight_id] for weight_id in weight_ids]
s_unc = [weights['S'][weight_id] for weight_id in weight_ids]
return p_unc, s_unc
else:
if not figure:
fig = plt.figure()
ax = fig.add_subplot(111)
# get percentage of weights
numPweights = np.sum([len(Pw0), len(Pw1), len(Pw2), len(Pw3), len(Pw4)])
numSweights = np.sum([len(Sw0), len(Sw1), len(Sw2), len(Sw3), len(Sw4)])
if len(Pw0) > 0:
P0perc = 100 / numPweights * len(Pw0)
else:
P0perc = 0
if len(Pw1) > 0:
P1perc = 100 / numPweights * len(Pw1)
else:
P1perc = 0
if len(Pw2) > 0:
P2perc = 100 / numPweights * len(Pw2)
else:
P2perc = 0
if len(Pw3) > 0:
P3perc = 100 / numPweights * len(Pw3)
else:
P3perc = 0
if len(Pw4) > 0:
P4perc = 100 / numPweights * len(Pw4)
else:
P4perc = 0
if len(Sw0) > 0:
S0perc = 100 / numSweights * len(Sw0)
else:
S0perc = 0
if len(Sw1) > 0:
S1perc = 100 / numSweights * len(Sw1)
else:
S1perc = 0
if len(Sw2) > 0:
S2perc = 100 / numSweights * len(Sw2)
else:
S2perc = 0
if len(Sw3) > 0:
S3perc = 100 / numSweights * len(Sw3)
else:
S3perc = 0
if len(Sw4) > 0:
S4perc = 100 / numSweights * len(Sw4)
else:
S4perc = 0
listP, numPweights = calc_weight_perc(weights['P'], weight_ids)
listS, numSweights = calc_weight_perc(weights['S'], weight_ids)
weights = ('0', '1', '2', '3', '4')
y_pos = np.arange(len(weights))
y_pos = np.arange(len(weight_ids))
width = 0.34
plt.bar(y_pos - width, [P0perc, P1perc, P2perc, P3perc, P4perc], width, color='black')
plt.bar(y_pos, [S0perc, S1perc, S2perc, S3perc, S4perc], width, color='red')
plt.ylabel('%')
plt.xticks(y_pos, weights)
plt.xlim([-0.5, 4.5])
plt.xlabel('Qualities')
plt.title('{0} P-Qualities, {1} S-Qualities'.format(numPweights, numSweights))
plt.show()
ax.bar(y_pos - width, listP, width, color='black')
ax.bar(y_pos, listS, width, color='red')
ax.set_ylabel('%')
ax.set_xticks(y_pos, weight_ids)
ax.set_xlim([-0.5, 4.5])
ax.set_xlabel('Qualities')
ax.set_title('{0} P-Qualities, {1} S-Qualities'.format(numPweights, numSweights))
if not figure:
fig.show()
return listP, listS
+2 -1
View File
@@ -4,11 +4,12 @@
import glob
import os
import subprocess
from obspy import read_events
from pylot.core.io.phases import writephases
from pylot.core.util.utils import getPatternLine, runProgram
from pylot.core.util.gui import which
from pylot.core.util.utils import getPatternLine, runProgram
from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString()
+84 -61
View File
@@ -9,21 +9,21 @@ function conglomerate utils.
:author: MAGS2 EP3 working group / Ludger Kueperkoch
"""
import copy
import traceback
import matplotlib.pyplot as plt
import numpy as np
import traceback
from obspy import Trace
from obspy.taup import TauPyModel
from pylot.core.pick.charfuns import CharacteristicFunction
from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
from pylot.core.pick.picker import AICPicker, PragPicker
from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
getSNR, fmpicker, checkPonsets, wadaticheck, get_pickparams, get_quality_class
from pylot.core.util.utils import getPatternLine, gen_Pool,\
getSNR, fmpicker, checkPonsets, wadaticheck, get_quality_class
from pylot.core.util.utils import getPatternLine, gen_Pool, \
get_Bool, identifyPhaseID, get_None, correct_iplot
from obspy.taup import TauPyModel
from obspy import Trace
def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None, ncores=0, metadata=None, origin=None):
"""
@@ -182,25 +182,25 @@ class PickingResults(dict):
# TODO What are those?
self.w0 = None
self.fc = None
self.Ao = None # Wood-Anderson peak-to-peak amplitude
self.Ao = None # Wood-Anderson peak-to-peak amplitude
# Station information
self.network = None
self.channel = None
# pick information
self.picker = 'auto' # type of pick
self.picker = 'auto' # type of pick
self.marked = []
# pick results
self.epp = None # earliest possible pick
self.mpp = None # most likely onset
self.lpp = None # latest possible pick
self.fm = 'N' # first motion polarity, can be set to 'U' (Up) or 'D' (Down)
self.snr = None # signal-to-noise ratio of onset
self.snrdb = None # signal-to-noise ratio of onset [dB]
self.spe = None # symmetrized picking error
self.weight = 4 # weight of onset
self.epp = None # earliest possible pick
self.mpp = None # most likely onset
self.lpp = None # latest possible pick
self.fm = 'N' # first motion polarity, can be set to 'U' (Up) or 'D' (Down)
self.snr = None # signal-to-noise ratio of onset
self.snrdb = None # signal-to-noise ratio of onset [dB]
self.spe = None # symmetrized picking error
self.weight = 4 # weight of onset
# to correctly provide dot access to dictionary attributes, all attribute access of the class is forwarded to the
# dictionary
@@ -335,9 +335,10 @@ class AutopickStation(object):
"""
waveform_data = {}
for key in self.channelorder:
waveform_data[key] = self.wfstream.select(component=key) # try ZNE first
waveform_data[key] = self.wfstream.select(component=key) # try ZNE first
if len(waveform_data[key]) == 0:
waveform_data[key] = self.wfstream.select(component=str(self.channelorder[key])) # use 123 as second option
waveform_data[key] = self.wfstream.select(
component=str(self.channelorder[key])) # use 123 as second option
return waveform_data['Z'], waveform_data['N'], waveform_data['E']
def get_traces_from_streams(self):
@@ -474,13 +475,9 @@ class AutopickStation(object):
# TODO here the pickparams is modified, instead of a copy
self.pickparams["pstart"] = 0
if self.pickparams["use_taup"] is False or not self.origin or not self.metadata:
if self.pickparams["use_taup"] is False:
# correct user mistake where a relative cuttime is selected (pstart < 0) but use of taupy is disabled/ has
# not the required parameters
if not self.origin:
print('Requested use_taup but no origin given. Exit taupy.')
if not self.metadata:
print('Requested use_taup but no metadata given. Exit taupy.')
exit_taupy()
return
@@ -528,7 +525,7 @@ class AutopickStation(object):
self.plot_pick_results()
self.finish_picking()
return [{'P': self.p_results, 'S':self.s_results}, self.ztrace.stats.station]
return [{'P': self.p_results, 'S': self.s_results}, self.ztrace.stats.station]
def finish_picking(self):
@@ -577,7 +574,7 @@ class AutopickStation(object):
self.s_results.channel = self.etrace.stats.channel
self.s_results.network = self.etrace.stats.network
self.s_results.fm = None # override default value 'N'
self.s_results.fm = None # override default value 'N'
def plot_pick_results(self):
if self.iplot > 0:
@@ -592,12 +589,14 @@ class AutopickStation(object):
plt_flag = 0
fig._tight = True
ax1 = fig.add_subplot(311)
tdata = np.linspace(start=0, stop=self.ztrace.stats.endtime-self.ztrace.stats.starttime, num=self.ztrace.stats.npts)
tdata = np.linspace(start=0, stop=self.ztrace.stats.endtime - self.ztrace.stats.starttime,
num=self.ztrace.stats.npts)
# plot tapered trace filtered with bpz2 filter settings
ax1.plot(tdata, self.tr_filt_z_bpz2.data/max(self.tr_filt_z_bpz2.data), color=linecolor, linewidth=0.7, label='Data')
ax1.plot(tdata, self.tr_filt_z_bpz2.data / max(self.tr_filt_z_bpz2.data), color=linecolor, linewidth=0.7,
label='Data')
if self.p_results.weight < 4:
# plot CF of initial onset (HOScf or ARZcf)
ax1.plot(self.cf1.getTimeArray(), self.cf1.getCF()/max(self.cf1.getCF()), 'b', label='CF1')
ax1.plot(self.cf1.getTimeArray(), self.cf1.getCF() / max(self.cf1.getCF()), 'b', label='CF1')
if self.p_data.p_aic_plot_flag == 1:
aicpick = self.p_data.aicpick
refPpick = self.p_data.refPpick
@@ -635,23 +634,28 @@ class AutopickStation(object):
if self.horizontal_traces_exist() and self.s_data.Sflag == 1:
# plot E trace
ax2 = fig.add_subplot(3, 1, 2, sharex=ax1)
th1data = np.linspace(0, self.etrace.stats.endtime-self.etrace.stats.starttime, self.etrace.stats.npts)
th1data = np.linspace(0, self.etrace.stats.endtime - self.etrace.stats.starttime,
self.etrace.stats.npts)
# plot filtered and tapered waveform
ax2.plot(th1data, self.etrace.data / max(self.etrace.data), color=linecolor, linewidth=0.7, label='Data')
ax2.plot(th1data, self.etrace.data / max(self.etrace.data), color=linecolor, linewidth=0.7,
label='Data')
if self.p_results.weight < 4:
# plot initial CF (ARHcf or AR3Ccf)
ax2.plot(self.arhcf1.getTimeArray(), self.arhcf1.getCF() / max(self.arhcf1.getCF()), 'b', label='CF1')
ax2.plot(self.arhcf1.getTimeArray(), self.arhcf1.getCF() / max(self.arhcf1.getCF()), 'b',
label='CF1')
if self.s_data.aicSflag == 1 and self.s_results.weight <= 4:
aicarhpick = self.aicarhpick
refSpick = self.refSpick
# plot second cf, used for determing precise onset (ARHcf or AR3Ccf)
ax2.plot(self.arhcf2.getTimeArray(), self.arhcf2.getCF() / max(self.arhcf2.getCF()), 'm', label='CF2')
ax2.plot(self.arhcf2.getTimeArray(), self.arhcf2.getCF() / max(self.arhcf2.getCF()), 'm',
label='CF2')
# plot preliminary onset time, calculated from CF1
ax2.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'g', label='Initial S Onset')
ax2.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [1, 1], 'g')
ax2.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [-1, -1], 'g')
# plot precise onset time, calculated from CF2
ax2.plot([refSpick.getpick(), refSpick.getpick()], [-1.3, 1.3], 'g', linewidth=2, label='Final S Pick')
ax2.plot([refSpick.getpick(), refSpick.getpick()], [-1.3, 1.3], 'g', linewidth=2,
label='Final S Pick')
ax2.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [1.3, 1.3], 'g', linewidth=2)
ax2.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [-1.3, -1.3], 'g', linewidth=2)
ax2.plot([self.s_results.lpp, self.s_results.lpp], [-1.1, 1.1], 'g--', label='lpp')
@@ -671,15 +675,19 @@ class AutopickStation(object):
# plot N trace
ax3 = fig.add_subplot(3, 1, 3, sharex=ax1)
th2data= np.linspace(0, self.ntrace.stats.endtime-self.ntrace.stats.starttime, self.ntrace.stats.npts)
th2data = np.linspace(0, self.ntrace.stats.endtime - self.ntrace.stats.starttime,
self.ntrace.stats.npts)
# plot trace
ax3.plot(th2data, self.ntrace.data / max(self.ntrace.data), color=linecolor, linewidth=0.7, label='Data')
ax3.plot(th2data, self.ntrace.data / max(self.ntrace.data), color=linecolor, linewidth=0.7,
label='Data')
if self.p_results.weight < 4:
p22, = ax3.plot(self.arhcf1.getTimeArray(), self.arhcf1.getCF() / max(self.arhcf1.getCF()), 'b', label='CF1')
p22, = ax3.plot(self.arhcf1.getTimeArray(), self.arhcf1.getCF() / max(self.arhcf1.getCF()), 'b',
label='CF1')
if self.s_data.aicSflag == 1:
aicarhpick = self.aicarhpick
refSpick = self.refSpick
ax3.plot(self.arhcf2.getTimeArray(), self.arhcf2.getCF() / max(self.arhcf2.getCF()), 'm', label='CF2')
ax3.plot(self.arhcf2.getTimeArray(), self.arhcf2.getCF() / max(self.arhcf2.getCF()), 'm',
label='CF2')
ax3.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'g', label='Initial S Onset')
ax3.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [1, 1], 'g')
ax3.plot([aicarhpick.getpick() - 0.5, aicarhpick.getpick() + 0.5], [-1, -1], 'g')
@@ -720,7 +728,8 @@ class AutopickStation(object):
if aicpick.getpick() is None:
msg = "Bad initial (AIC) P-pick, skipping this onset!\nAIC-SNR={0}, AIC-Slope={1}counts/s\n " \
"(min. AIC-SNR={2}, min. AIC-Slope={3}counts/s)"
msg = msg.format(aicpick.getSNR(), aicpick.getSlope(), self.pickparams["minAICPSNR"], self.pickparams["minAICPslope"])
msg = msg.format(aicpick.getSNR(), aicpick.getSlope(), self.pickparams["minAICPSNR"],
self.pickparams["minAICPslope"])
self.vprint(msg)
return 0
# Quality check initial pick with minimum signal length
@@ -730,14 +739,16 @@ class AutopickStation(object):
if len(self.nstream) == 0 or len(self.estream) == 0:
msg = 'One or more horizontal component(s) missing!\n' \
'Signal length only checked on vertical component!\n' \
'Decreasing minsiglengh from {0} to {1}'\
.format(minsiglength, minsiglength / 2)
'Decreasing minsiglengh from {0} to {1}' \
.format(minsiglength, minsiglength / 2)
self.vprint(msg)
minsiglength = minsiglength / 2
else:
# filter, taper other traces as well since signal length is compared on all traces
trH1_filt, _ = self.prepare_wfstream(self.estream, freqmin=self.pickparams["bph1"][0], freqmax=self.pickparams["bph1"][1])
trH2_filt, _ = self.prepare_wfstream(self.nstream, freqmin=self.pickparams["bph1"][0], freqmax=self.pickparams["bph1"][1])
trH1_filt, _ = self.prepare_wfstream(self.estream, freqmin=self.pickparams["bph1"][0],
freqmax=self.pickparams["bph1"][1])
trH2_filt, _ = self.prepare_wfstream(self.nstream, freqmin=self.pickparams["bph1"][0],
freqmax=self.pickparams["bph1"][1])
zne += trH1_filt
zne += trH2_filt
minsiglength = minsiglength
@@ -823,15 +834,18 @@ class AutopickStation(object):
# get preliminary onset time from AIC-CF
self.set_current_figure('aicFig')
aicpick = AICPicker(aiccf, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot,
Tsmooth=self.pickparams["aictsmooth"], fig=self.current_figure, linecolor=self.current_linecolor)
Tsmooth=self.pickparams["aictsmooth"], fig=self.current_figure,
linecolor=self.current_linecolor)
# save aicpick for plotting later
self.p_data.aicpick = aicpick
# add pstart and pstop to aic plot
if self.current_figure:
# TODO remove plotting from picking, make own plot function
for ax in self.current_figure.axes:
ax.vlines(self.pickparams["pstart"], ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed', label='P start')
ax.vlines(self.pickparams["pstop"], ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed', label='P stop')
ax.vlines(self.pickparams["pstart"], ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed',
label='P start')
ax.vlines(self.pickparams["pstop"], ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed',
label='P stop')
ax.legend(loc=1)
Pflag = self._pick_p_quality_control(aicpick, z_copy, tr_filt)
@@ -845,7 +859,8 @@ class AutopickStation(object):
error_msg = 'AIC P onset slope to small: got {}, min {}'.format(slope, self.pickparams["minAICPslope"])
raise PickingFailedException(error_msg)
if aicpick.getSNR() < self.pickparams["minAICPSNR"]:
error_msg = 'AIC P onset SNR to small: got {}, min {}'.format(aicpick.getSNR(), self.pickparams["minAICPSNR"])
error_msg = 'AIC P onset SNR to small: got {}, min {}'.format(aicpick.getSNR(),
self.pickparams["minAICPSNR"])
raise PickingFailedException(error_msg)
self.p_data.p_aic_plot_flag = 1
@@ -853,7 +868,8 @@ class AutopickStation(object):
'autopickstation: re-filtering vertical trace...'.format(aicpick.getSlope(), aicpick.getSNR())
self.vprint(msg)
# refilter waveform with larger bandpass
tr_filt, z_copy = self.prepare_wfstream(self.zstream, freqmin=self.pickparams["bpz2"][0], freqmax=self.pickparams["bpz2"][1])
tr_filt, z_copy = self.prepare_wfstream(self.zstream, freqmin=self.pickparams["bpz2"][0],
freqmax=self.pickparams["bpz2"][1])
# save filtered trace in instance for later plotting
self.tr_filt_z_bpz2 = tr_filt
# determine new times around initial onset
@@ -865,25 +881,29 @@ class AutopickStation(object):
else:
self.cf2 = None
assert isinstance(self.cf2, CharacteristicFunction), 'cf2 is not set correctly: maybe the algorithm name () is ' \
'corrupted'.format(self.pickparams["algoP"])
'corrupted'.format(self.pickparams["algoP"])
self.set_current_figure('refPpick')
# get refined onset time from CF2
refPpick = PragPicker(self.cf2, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot, self.pickparams["ausP"],
self.pickparams["tsmoothP"], aicpick.getpick(), self.current_figure, self.current_linecolor)
refPpick = PragPicker(self.cf2, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot,
self.pickparams["ausP"],
self.pickparams["tsmoothP"], aicpick.getpick(), self.current_figure,
self.current_linecolor)
# save PragPicker result for plotting
self.p_data.refPpick = refPpick
self.p_results.mpp = refPpick.getpick()
if self.p_results.mpp is None:
msg = 'Bad initial (AIC) P-pick, skipping this onset!\n AIC-SNR={}, AIC-Slope={}counts/s\n' \
'(min. AIC-SNR={}, min. AIC-Slope={}counts/s)'
msg.format(aicpick.getSNR(), aicpick.getSlope(), self.pickparams["minAICPSNR"], self.pickparams["minAICPslope"])
msg.format(aicpick.getSNR(), aicpick.getSlope(), self.pickparams["minAICPSNR"],
self.pickparams["minAICPslope"])
self.vprint(msg)
self.s_data.Sflag = 0
raise PickingFailedException(msg)
# quality assessment, get earliest/latest pick and symmetrized uncertainty
#todo quality assessment in own function
# todo quality assessment in own function
self.set_current_figure('el_Ppick')
elpicker_results = earllatepicker(z_copy, self.pickparams["nfacP"], self.pickparams["tsnrz"], self.p_results.mpp,
elpicker_results = earllatepicker(z_copy, self.pickparams["nfacP"], self.pickparams["tsnrz"],
self.p_results.mpp,
self.iplot, fig=self.current_figure, linecolor=self.current_linecolor)
self.p_results.epp, self.p_results.lpp, self.p_results.spe = elpicker_results
snr_results = getSNR(z_copy, self.pickparams["tsnrz"], self.p_results.mpp)
@@ -891,7 +911,8 @@ class AutopickStation(object):
# weight P-onset using symmetric error
self.p_results.weight = get_quality_class(self.p_results.spe, self.pickparams["timeerrorsP"])
if self.p_results.weight <= self.pickparams["minfmweight"] and self.p_results.snr >= self.pickparams["minFMSNR"]:
if self.p_results.weight <= self.pickparams["minfmweight"] and self.p_results.snr >= self.pickparams[
"minFMSNR"]:
# if SNR is high enough, try to determine first motion of onset
self.set_current_figure('fm_picker')
self.p_results.fm = fmpicker(self.zstream, z_copy, self.pickparams["fmpickwin"], self.p_results.mpp,
@@ -964,7 +985,7 @@ class AutopickStation(object):
trH1_filt, _ = self.prepare_wfstream(self.zstream, filter_freq_min, filter_freq_max)
trH2_filt, _ = self.prepare_wfstream(self.estream, filter_freq_min, filter_freq_max)
trH3_filt, _ = self.prepare_wfstream(self.nstream, filter_freq_min, filter_freq_max)
h_copy =self. hdat.copy()
h_copy = self.hdat.copy()
h_copy[0].data = trH1_filt.data
h_copy[1].data = trH2_filt.data
h_copy[2].data = trH3_filt.data
@@ -1119,7 +1140,8 @@ class AutopickStation(object):
# get preliminary onset time from AIC cf
self.set_current_figure('aicARHfig')
aicarhpick = AICPicker(haiccf, self.pickparams["tsnrh"], self.pickparams["pickwinS"], self.iplot,
Tsmooth=self.pickparams["aictsmoothS"], fig=self.current_figure, linecolor=self.current_linecolor)
Tsmooth=self.pickparams["aictsmoothS"], fig=self.current_figure,
linecolor=self.current_linecolor)
# save pick for later plotting
self.aicarhpick = aicarhpick
@@ -1130,8 +1152,10 @@ class AutopickStation(object):
# get refined onset time from CF2
self.set_current_figure('refSpick')
refSpick = PragPicker(arhcf2, self.pickparams["tsnrh"], self.pickparams["pickwinS"], self.iplot, self.pickparams["ausS"],
self.pickparams["tsmoothS"], aicarhpick.getpick(), self.current_figure, self.current_linecolor)
refSpick = PragPicker(arhcf2, self.pickparams["tsnrh"], self.pickparams["pickwinS"], self.iplot,
self.pickparams["ausS"],
self.pickparams["tsmoothS"], aicarhpick.getpick(), self.current_figure,
self.current_linecolor)
# save refSpick for later plotitng
self.refSpick = refSpick
self.s_results.mpp = refSpick.getpick()
@@ -1155,7 +1179,6 @@ class AutopickStation(object):
self.current_linecolor = plot_style['linecolor']['rgba_mpl']
def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None, metadata=None, origin=None):
"""
Main function to calculate picks for the station.
@@ -1243,11 +1266,11 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter, fig_dict=None):
print(
"iteratepicker: The following picking parameters have been modified for iterative picking:")
print(
"pstart: %fs => %fs" % (pstart_old, pickparameter.get('pstart')))
"pstart: %fs => %fs" % (pstart_old, pickparameter.get('pstart')))
print(
"pstop: %fs => %fs" % (pstop_old, pickparameter.get('pstop')))
"pstop: %fs => %fs" % (pstop_old, pickparameter.get('pstop')))
print(
"sstop: %fs => %fs" % (sstop_old, pickparameter.get('sstop')))
"sstop: %fs => %fs" % (sstop_old, pickparameter.get('sstop')))
print("pickwinP: %fs => %fs" % (
pickwinP_old, pickparameter.get('pickwinP')))
print("Precalcwin: %fs => %fs" % (
+13 -10
View File
@@ -18,8 +18,8 @@ autoregressive prediction: application ot local and regional distances, Geophys.
"""
import numpy as np
from scipy import signal
from obspy.core import Stream
from scipy import signal
class CharacteristicFunction(object):
@@ -150,7 +150,7 @@ class CharacteristicFunction(object):
if self.cut[0] == 0 and self.cut[1] == 0:
start = 0
stop = len(self.orig_data[0])
elif self.cut[0] == 0 and self.cut[1] is not 0:
elif self.cut[0] == 0 and self.cut[1] != 0:
start = 0
stop = self.cut[1] / self.dt
else:
@@ -159,7 +159,7 @@ class CharacteristicFunction(object):
zz = self.orig_data.copy()
z1 = zz[0].copy()
zz[0].data = z1.data[int(start):int(stop)]
if zz[0].stats.npts == 0: # cut times do not fit data length!
if zz[0].stats.npts == 0: # cut times do not fit data length!
zz[0].data = z1.data # take entire data
data = zz
return data
@@ -167,7 +167,7 @@ class CharacteristicFunction(object):
if self.cut[0] == 0 and self.cut[1] == 0:
start = 0
stop = min([len(self.orig_data[0]), len(self.orig_data[1])])
elif self.cut[0] == 0 and self.cut[1] is not 0:
elif self.cut[0] == 0 and self.cut[1] != 0:
start = 0
stop = min([self.cut[1] / self.dt, len(self.orig_data[0]),
len(self.orig_data[1])])
@@ -187,7 +187,7 @@ class CharacteristicFunction(object):
start = 0
stop = min([self.cut[1] / self.dt, len(self.orig_data[0]),
len(self.orig_data[1]), len(self.orig_data[2])])
elif self.cut[0] == 0 and self.cut[1] is not 0:
elif self.cut[0] == 0 and self.cut[1] != 0:
start = 0
stop = self.cut[1] / self.dt
else:
@@ -241,7 +241,7 @@ class AICcf(CharacteristicFunction):
ff = np.where(inf is True)
if len(ff) >= 1:
cf[ff] = 0
self.cf = cf - np.mean(cf)
self.xcf = x
@@ -305,7 +305,7 @@ class HOScf(CharacteristicFunction):
if ind.size:
first = ind[0]
LTA[:first] = LTA[first]
self.cf = LTA
self.xcf = x
@@ -313,7 +313,8 @@ class HOScf(CharacteristicFunction):
class ARZcf(CharacteristicFunction):
def __init__(self, data, cut, t1, t2, pickparams):
super(ARZcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Parorder"], fnoise=pickparams["addnoise"])
super(ARZcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Parorder"],
fnoise=pickparams["addnoise"])
def calcCF(self, data):
"""
@@ -448,7 +449,8 @@ class ARZcf(CharacteristicFunction):
class ARHcf(CharacteristicFunction):
def __init__(self, data, cut, t1, t2, pickparams):
super(ARHcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"], fnoise=pickparams["addnoise"])
super(ARHcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
fnoise=pickparams["addnoise"])
def calcCF(self, data):
"""
@@ -600,7 +602,8 @@ class ARHcf(CharacteristicFunction):
class AR3Ccf(CharacteristicFunction):
def __init__(self, data, cut, t1, t2, pickparams):
super(AR3Ccf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"], fnoise=pickparams["addnoise"])
super(AR3Ccf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
fnoise=pickparams["addnoise"])
def calcCF(self, data):
"""
+6 -3
View File
@@ -2,10 +2,11 @@
# -*- coding: utf-8 -*-
import copy
import matplotlib.pyplot as plt
import numpy as np
import operator
import os
import matplotlib.pyplot as plt
import numpy as np
from obspy.core import AttribDict
from pylot.core.util.pdf import ProbabilityDensityFunction
@@ -117,7 +118,7 @@ class Comparison(object):
pdf_a = self.get('auto').generate_pdf_data(type)
pdf_b = self.get('manu').generate_pdf_data(type)
for station, phases in pdf_a.items():
if station in pdf_b.keys():
compare_pdf = dict()
@@ -401,6 +402,8 @@ class PDFstatistics(object):
Takes a path as argument.
"""
# TODO: change root to datapath
def __init__(self, directory):
"""Initiates some values needed when dealing with pdfs later"""
self._rootdir = directory
+3 -2
View File
@@ -19,9 +19,10 @@ calculated after Diehl & Kissling (2009).
:author: MAGS2 EP3 working group / Ludger Kueperkoch
"""
import warnings
import matplotlib.pyplot as plt
import numpy as np
import warnings
from scipy.signal import argrelmax, argrelmin
from pylot.core.pick.charfuns import CharacteristicFunction
@@ -476,7 +477,7 @@ class PragPicker(AutoPicker):
cfpick_r = 0
cfpick_l = 0
lpickwindow = int(round(self.PickWindow / self.dt))
#for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])):
# for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])):
# # local minimum
# if self.cf[i + 1] > self.cf[i] <= self.cf[i - 1]:
# if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
+22 -10
View File
@@ -9,10 +9,11 @@
"""
import warnings
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import argrelmax
from obspy.core import Stream, UTCDateTime
from scipy.signal import argrelmax
from pylot.core.util.utils import get_Bool, get_None, SetChannelComponents
@@ -73,7 +74,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol
x = X[0].data
t = np.linspace(0, X[0].stats.endtime - X[0].stats.starttime,
X[0].stats.npts)
X[0].stats.npts)
inoise = getnoisewin(t, Pick1, TSNR[0], TSNR[1])
# get signal window
isignal = getsignalwin(t, Pick1, TSNR[2])
@@ -218,7 +219,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
xraw = Xraw[0].data
xfilt = Xfilt[0].data
t = np.linspace(0, Xraw[0].stats.endtime - Xraw[0].stats.starttime,
Xraw[0].stats.npts)
Xraw[0].stats.npts)
# get pick window
ipick = np.where((t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick))
if len(ipick[0]) <= 1:
@@ -536,9 +537,10 @@ def getslopewin(Tcf, Pick, tslope):
:rtype: `numpy.ndarray`
"""
# TODO: fill out docstring
slope = np.where( (Tcf <= min(Pick + tslope, Tcf[-1])) & (Tcf >= Pick) )
slope = np.where((Tcf <= min(Pick + tslope, Tcf[-1])) & (Tcf >= Pick))
return slope[0]
def getResolutionWindow(snr, extent):
"""
Produce the half of the time resolution window width from given SNR value
@@ -848,7 +850,7 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
print("Presumably picked noise peak, pick is rejected!")
print("(min. signal length required: %s s)" % minsiglength)
returnflag = 0
else:
else:
# calculate minimum adjusted signal level
minsiglevel = np.mean(rms[inoise]) * nfac
# minimum adjusted number of samples over minimum signal level
@@ -1207,7 +1209,7 @@ def checkZ4S(X, pick, pickparams, iplot, fig=None, linecolor='k'):
rms = rms_dict[key]
trace = traces_dict[key]
t = np.linspace(diff_dict[key], trace.stats.endtime - trace.stats.starttime + diff_dict[key],
trace.stats.npts)
trace.stats.npts)
if i == 0:
if get_None(fig) is None:
fig = plt.figure() # self.iplot) ### WHY? MP MP
@@ -1318,6 +1320,7 @@ def get_quality_class(uncertainty, weight_classes):
:return: quality of pick (0-4)
:rtype: int
"""
if not uncertainty: return len(weight_classes)
try:
# create generator expression containing all indices of values in weight classes that are >= than uncertainty.
# call next on it once to receive first value
@@ -1328,6 +1331,7 @@ def get_quality_class(uncertainty, weight_classes):
quality = len(weight_classes)
return quality
def set_NaNs_to(data, nan_value):
"""
Replace all NaNs in data with nan_value
@@ -1343,6 +1347,7 @@ def set_NaNs_to(data, nan_value):
data[nn] = nan_value
return data
def taper_cf(cf):
"""
Taper cf data to get rid off of side maximas
@@ -1354,6 +1359,7 @@ def taper_cf(cf):
tap = np.hanning(len(cf))
return tap * cf
def cf_positive(cf):
"""
Shifts cf so that all values are positive
@@ -1364,6 +1370,7 @@ def cf_positive(cf):
"""
return cf + max(abs(cf))
def smooth_cf(cf, t_smooth, delta):
"""
Smooth cf by taking samples over t_smooth length
@@ -1392,6 +1399,7 @@ def smooth_cf(cf, t_smooth, delta):
cf_smooth -= offset # remove offset from smoothed function
return cf_smooth
def check_counts_ms(data):
"""
check if data is in counts or m/s
@@ -1451,9 +1459,9 @@ def calcSlope(Data, datasmooth, Tcf, Pick, TSNR):
if imax == 0:
print("AICPicker: Maximum for slope determination right at the beginning of the window!")
print("Choose longer slope determination window!")
raise IndexError
raise IndexError
iislope = islope[0][0:imax + 1] # cut index so it contains only the first maximum
dataslope = Data[0].data[iislope] # slope will only be calculated to the first maximum
dataslope = Data[0].data[iislope] # slope will only be calculated to the first maximum
# calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope))
P = np.polyfit(xslope, dataslope, 1)
@@ -1474,8 +1482,10 @@ def get_pickparams(pickparam):
:rtype: (dict, dict, dict, dict)
"""
# Define names of all parameters in different groups
p_parameter_names = 'algoP pstart pstop use_taup taup_model tlta tsnrz hosorder bpz1 bpz2 pickwinP aictsmooth tsmoothP ausP nfacP tpred1z tdet1z Parorder addnoise Precalcwin minAICPslope minAICPSNR timeerrorsP checkwindowP minfactorP'.split(' ')
s_parameter_names = 'algoS sstart sstop bph1 bph2 tsnrh pickwinS tpred1h tdet1h tpred2h tdet2h Sarorder aictsmoothS tsmoothS ausS minAICSslope minAICSSNR Srecalcwin nfacS timeerrorsS zfac checkwindowS minfactorS'.split(' ')
p_parameter_names = 'algoP pstart pstop use_taup taup_model tlta tsnrz hosorder bpz1 bpz2 pickwinP aictsmooth tsmoothP ausP nfacP tpred1z tdet1z Parorder addnoise Precalcwin minAICPslope minAICPSNR timeerrorsP checkwindowP minfactorP'.split(
' ')
s_parameter_names = 'algoS sstart sstop bph1 bph2 tsnrh pickwinS tpred1h tdet1h tpred2h tdet2h Sarorder aictsmoothS tsmoothS ausS minAICSslope minAICSSNR Srecalcwin nfacS timeerrorsS zfac checkwindowS minfactorS'.split(
' ')
first_motion_names = 'minFMSNR fmpickwin minfmweight'.split(' ')
signal_length_names = 'minsiglength minpercent noisefactor'.split(' ')
# Get list of values from pickparam by name
@@ -1493,6 +1503,7 @@ def get_pickparams(pickparam):
return p_params, s_params, first_motion_params, signal_length_params
def getQualityFromUncertainty(uncertainty, Errors):
# set initial quality to 4 (worst) and change only if one condition is hit
quality = 4
@@ -1516,6 +1527,7 @@ def getQualityFromUncertainty(uncertainty, Errors):
return quality
if __name__ == '__main__':
import doctest
+341 -294
View File
@@ -1,38 +1,51 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os
import traceback
import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib
import matplotlib.patheffects as PathEffects
import matplotlib.pyplot as plt
import numpy as np
import obspy
import traceback
from PySide import QtGui
from matplotlib.figure import Figure
from PySide2 import QtWidgets
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from mpl_toolkits.basemap import Basemap
from scipy.interpolate import griddata
from pylot.core.util.widgets import PickDlg, PylotCanvas
from pylot.core.pick.utils import get_quality_class
from pylot.core.util.widgets import PickDlg
plt.interactive(False)
matplotlib.use('Qt5Agg')
class Array_map(QtGui.QWidget):
def __init__(self, parent, metadata, parameter=None, figure=None, annotate=True, pointsize=25.,
class MplCanvas(FigureCanvas):
def __init__(self, parent=None, extern_axes=None, width=5, height=4, dpi=100):
if extern_axes is None:
self.fig = plt.figure(figsize=(width, height), dpi=dpi)
self.axes = self.fig.add_subplot(111)
else:
self.fig = extern_axes.figure
self.axes = extern_axes
super(MplCanvas, self).__init__(self.fig)
class Array_map(QtWidgets.QWidget):
def __init__(self, parent, metadata, parameter=None, axes=None, annotate=True, pointsize=25.,
linewidth=1.5, width=5e6, height=2e6):
'''
Create a map of the array.
:param parent: object of PyLoT Mainwindow class
:param parameter: object of PyLoT parameter class
:param figure:
'''
QtGui.QWidget.__init__(self)
assert (parameter != None or parent != None), 'either parent or parameter has to be set'
QtWidgets.QWidget.__init__(self, parent=parent)
assert (parameter is not None or parent is not None), 'either parent or parameter has to be set'
# set properties
self._parent = parent
self.metadata = metadata
self.pointsize = pointsize
self.linewidth = linewidth
self.extern_plot_axes = axes
self.width = width
self.height = height
self.annotate = annotate
@@ -43,34 +56,228 @@ class Array_map(QtGui.QWidget):
self.hybrids_dict = None
self.eventLoc = None
self.parameter = parameter if parameter else parent._inputs
self.figure = figure
self.picks_rel = {}
self.marked_stations = []
self.highlighted_stations = []
# call functions to draw everything
self.init_graphics()
self.init_stations()
self.init_basemap(resolution='l')
self.init_crtpyMap()
self.init_map()
self._style = None if not hasattr(parent, '_style') else parent._style
# self.show()
# set original map limits to fall back on when home button is pressed
self.org_xlim = self.canvas.axes.get_xlim()
self.org_ylim = self.canvas.axes.get_ylim()
# initial map without event
self.canvas.axes.set_xlim(self.org_xlim[0], self.org_xlim[1])
self.canvas.axes.set_ylim(self.org_ylim[0], self.org_ylim[1])
def update_hybrids_dict(self):
self.hybrids_dict = self.picks_dict.copy()
for station, pick in self.autopicks_dict.items():
if not station in self.hybrids_dict.keys():
self.hybrids_dict[station] = pick
return self.hybrids_dict
self._style = None if not hasattr(parent, '_style') else parent._style
self.show()
def init_map(self):
self.init_colormap()
self.connectSignals()
self.draw_everything()
self.canvas.setZoomBorders2content()
def init_graphics(self):
"""
Initializes all GUI components and figure elements to be populated by other functions
"""
# initialize figure elements
if self.extern_plot_axes is None:
self.canvas = MplCanvas(self)
self.plotWidget = FigureCanvas(self.canvas.fig)
else:
self.canvas = MplCanvas(self, extern_axes=self.extern_plot_axes)
self.plotWidget = FigureCanvas(self.canvas.fig)
# initialize GUI elements
self.status_label = QtWidgets.QLabel()
self.map_reset_button = QtWidgets.QPushButton('Reset Map View')
self.save_map_button = QtWidgets.QPushButton('Save Map')
self.go2eq_button = QtWidgets.QPushButton('Go to Event Location')
self.main_box = QtWidgets.QVBoxLayout()
self.setLayout(self.main_box)
self.top_row = QtWidgets.QHBoxLayout()
self.main_box.addLayout(self.top_row, 1)
self.comboBox_phase = QtWidgets.QComboBox()
self.comboBox_phase.insertItem(0, 'P')
self.comboBox_phase.insertItem(1, 'S')
self.comboBox_am = QtWidgets.QComboBox()
self.comboBox_am.insertItem(0, 'hybrid (prefer manual)')
self.comboBox_am.insertItem(1, 'manual')
self.comboBox_am.insertItem(2, 'auto')
self.annotations_box = QtWidgets.QCheckBox('Annotate')
self.annotations_box.setChecked(True)
self.auto_refresh_box = QtWidgets.QCheckBox('Automatic refresh')
self.auto_refresh_box.setChecked(True)
self.refresh_button = QtWidgets.QPushButton('Refresh')
self.cmaps_box = QtWidgets.QComboBox()
self.cmaps_box.setMaxVisibleItems(20)
[self.cmaps_box.addItem(map_name) for map_name in sorted(plt.colormaps())]
# try to set to hsv as default
self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('hsv'))
self.top_row.addWidget(QtWidgets.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(QtWidgets.QLabel('Pick type: '))
self.top_row.addWidget(self.comboBox_am)
self.top_row.setStretch(3, 1) # set stretch of item 1 to 1
self.top_row.addWidget(self.cmaps_box)
self.top_row.addWidget(self.annotations_box)
self.top_row.addWidget(self.auto_refresh_box)
self.top_row.addWidget(self.refresh_button)
self.main_box.addWidget(self.plotWidget, 1)
self.bot_row = QtWidgets.QHBoxLayout()
self.main_box.addLayout(self.bot_row, 0.3)
self.bot_row.addWidget(QtWidgets.QLabel(''), 5)
self.bot_row.addWidget(self.map_reset_button, 2)
self.bot_row.addWidget(self.go2eq_button, 2)
self.bot_row.addWidget(self.save_map_button, 2)
self.bot_row.addWidget(self.status_label, 5)
def init_colormap(self):
self.init_lat_lon_dimensions()
self.init_lat_lon_grid()
self.init_x_y_dimensions()
def init_crtpyMap(self):
self.canvas.axes.cla()
self.canvas.axes = plt.axes(projection=ccrs.PlateCarree())
self.canvas.axes.add_feature(cf.LAND)
self.canvas.axes.add_feature(cf.OCEAN)
self.canvas.axes.add_feature(cf.COASTLINE, linewidth=1, edgecolor='gray')
self.canvas.axes.add_feature(cf.BORDERS, alpha=0.7)
self.canvas.axes.add_feature(cf.LAKES, alpha=0.7)
self.canvas.axes.add_feature(cf.RIVERS, linewidth=1)
# parallels and meridians
self.add_merid_paral()
self.canvas.fig.tight_layout()
def add_merid_paral(self):
self.gridlines = self.canvas.axes.gridlines(draw_labels=False, alpha=0.6, color='gray',
linewidth=self.linewidth / 2, zorder=7)
# TODO: current cartopy version does not support label removal. Devs are working on it.
# Should be fixed in coming cartopy versions
# self.gridlines.xformatter = LONGITUDE_FORMATTER
# self.gridlines.yformatter = LATITUDE_FORMATTER
def remove_merid_paral(self):
if len(self.gridlines.xline_artists):
self.gridlines.xline_artists[0].remove()
self.gridlines.yline_artists[0].remove()
def org_map_view(self):
self.canvas.axes.set_xlim(self.org_xlim[0], self.org_xlim[1])
self.canvas.axes.set_ylim(self.org_ylim[0], self.org_ylim[1])
# parallels and meridians
self.remove_merid_paral()
self.add_merid_paral()
self.canvas.axes.figure.canvas.draw_idle()
def go2eq(self):
if self.eventLoc:
lats, lons = self.eventLoc
self.canvas.axes.set_xlim(lons - 10, lons + 10)
self.canvas.axes.set_ylim(lats - 5, lats + 5)
# parallels and meridians
self.remove_merid_paral()
self.add_merid_paral()
self.canvas.axes.figure.canvas.draw_idle()
else:
self.status_label.setText('No event information available')
def connectSignals(self):
self.comboBox_phase.currentIndexChanged.connect(self._refresh_drawings)
self.comboBox_am.currentIndexChanged.connect(self._refresh_drawings)
self.cmaps_box.currentIndexChanged.connect(self._refresh_drawings)
self.annotations_box.stateChanged.connect(self.switch_annotations)
self.refresh_button.clicked.connect(self._refresh_drawings)
self.map_reset_button.clicked.connect(self.org_map_view)
self.go2eq_button.clicked.connect(self.go2eq)
self.save_map_button.clicked.connect(self.saveFigure)
self.plotWidget.mpl_connect('motion_notify_event', self.mouse_moved)
self.plotWidget.mpl_connect('scroll_event', self.mouse_scroll)
self.plotWidget.mpl_connect('button_press_event', self.mouseLeftPress)
self.plotWidget.mpl_connect('button_release_event', self.mouseLeftRelease)
# set mouse events -----------------------------------------------------
def mouse_moved(self, event):
if not event.inaxes == self.canvas.axes:
return
lat = event.ydata
lon = event.xdata
self.status_label.setText('Latitude: {:3.5f}, Longitude: {:3.5f}'.format(lat, lon))
def mouse_scroll(self, event):
if not event.inaxes == self.canvas.axes:
return
zoom = {'up': 1. / 2., 'down': 2.}
if event.button in zoom:
xlim = self.canvas.axes.get_xlim()
ylim = self.canvas.axes.get_ylim()
x, y = event.xdata, event.ydata
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
self.canvas.axes.set_xlim(xl, xr)
self.canvas.axes.set_ylim(yb, yt)
# parallels and meridians
self.remove_merid_paral()
self.add_merid_paral()
self.canvas.axes.figure.canvas.draw_idle()
def mouseLeftPress(self, event):
if not event.inaxes == self.canvas.axes:
return
self.map_x = event.xdata
self.map_y = event.ydata
self.map_xlim = self.canvas.axes.get_xlim()
self.map_ylim = self.canvas.axes.get_ylim()
def mouseLeftRelease(self, event):
if not event.inaxes == self.canvas.axes:
return
new_x = event.xdata
new_y = event.ydata
dx = new_x - self.map_x
dy = new_y - self.map_y
self.canvas.axes.set_xlim((self.map_xlim[0] - dx, self.map_xlim[1] - dx))
self.canvas.axes.set_ylim(self.map_ylim[0] - dy, self.map_ylim[1] - dy)
# parallels and meridians
self.remove_merid_paral()
self.add_merid_paral()
self.canvas.axes.figure.canvas.draw_idle()
def onpick(self, event):
ind = event.ind
@@ -84,6 +291,14 @@ class Array_map(QtGui.QWidget):
elif button == 3:
self.pickInfo(ind)
# data handling -----------------------------------------------------
def update_hybrids_dict(self):
self.hybrids_dict = self.picks_dict.copy()
for station, pick in self.autopicks_dict.items():
if not station in self.hybrids_dict.keys():
self.hybrids_dict[station] = pick
return self.hybrids_dict
def deletePick(self, ind):
self.update_hybrids_dict()
for index in ind:
@@ -117,15 +332,6 @@ class Array_map(QtGui.QWidget):
except Exception as e:
print('Could not delete pick for station {}.{}: {}'.format(network, station, e))
def highlight_station(self, network, station, color):
stat_dict = self.stations_dict['{}.{}'.format(network, station)]
lat = stat_dict['latitude']
lon = stat_dict['longitude']
self.highlighted_stations.append(self.basemap.scatter(lon, lat, s=self.pointsize, edgecolors=color,
facecolors='none', latlon=True,
zorder=12, label='deleted'))
self.canvas.draw()
def pickInfo(self, ind):
self.update_hybrids_dict()
for index in ind:
@@ -139,57 +345,6 @@ class Array_map(QtGui.QWidget):
for key, info in picks.items():
print('{}: {}'.format(key, info))
def openPickDlg(self, ind):
data = self._parent.get_data().getWFData()
for index in ind:
network, station = self._station_onpick_ids[index].split('.')[:2]
pyl_mw = self._parent
try:
data = data.select(station=station)
if not data:
self._warn('No data for station {}'.format(station))
return
pickDlg = PickDlg(self._parent, parameter=self.parameter,
data=data, network=network, station=station,
picks=self._parent.get_current_event().getPick(station),
autopicks=self._parent.get_current_event().getAutopick(station),
filteroptions=self._parent.filteroptions, metadata=self.metadata,
event=pyl_mw.get_current_event())
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)
print(traceback.format_exc())
return
try:
if pickDlg.exec_():
pyl_mw.setDirty(True)
pyl_mw.update_status('picks accepted ({0})'.format(station))
pyl_mw.addPicks(station, pickDlg.getPicks(picktype='manual'), type='manual')
pyl_mw.addPicks(station, pickDlg.getPicks(picktype='auto'), type='auto')
if self.auto_refresh_box.isChecked():
self._refresh_drawings()
else:
self.highlight_station(network, station, color='yellow')
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)
print(traceback.format_exc())
def connectSignals(self):
self.comboBox_phase.currentIndexChanged.connect(self._refresh_drawings)
self.comboBox_am.currentIndexChanged.connect(self._refresh_drawings)
self.cmaps_box.currentIndexChanged.connect(self._refresh_drawings)
self.annotations_box.stateChanged.connect(self.switch_annotations)
self.refresh_button.clicked.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]
@@ -205,14 +360,6 @@ class Array_map(QtGui.QWidget):
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().split(' ')[0]
auto_manu = {'auto': self.autopicks_dict,
@@ -220,57 +367,6 @@ class Array_map(QtGui.QWidget):
'hybrid': self.hybrids_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.main_ax.set_facecolor('0.7')
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, 'hybrid (prefer manual)')
self.comboBox_am.insertItem(1, 'manual')
self.comboBox_am.insertItem(2, 'auto')
self.annotations_box = QtGui.QCheckBox('Annotate')
self.annotations_box.setChecked(True)
self.auto_refresh_box = QtGui.QCheckBox('Automatic refresh')
self.auto_refresh_box.setChecked(True)
self.refresh_button = QtGui.QPushButton('Refresh')
self.cmaps_box = QtGui.QComboBox()
self.cmaps_box.setMaxVisibleItems(20)
[self.cmaps_box.addItem(map_name) for map_name in sorted(plt.colormaps())]
# try to set to hsv as default
self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('hsv'))
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.top_row.addWidget(self.cmaps_box)
self.top_row.addWidget(self.annotations_box)
self.top_row.addWidget(self.auto_refresh_box)
self.top_row.addWidget(self.refresh_button)
self.main_box.addWidget(self.canvas, 1)
self.main_box.addWidget(self.status_label, 0)
def init_stations(self):
self.stations_dict = self.metadata.get_all_coordinates()
self.latmin = self.get_min_from_stations('latitude')
@@ -323,41 +419,6 @@ class Array_map(QtGui.QWidget):
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)
basemap = Basemap(projection='lcc', resolution=resolution, ax=self.main_ax,
width=self.width, height=self.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)#, labels=[1, 1, 0, 0])
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)#, labels=[0, 0, 1, 1])
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)
@@ -367,8 +428,7 @@ class Array_map(QtGui.QWidget):
def init_picksgrid(self):
picks, uncertainties, lats, lons = self.get_picks_lat_lon()
try:
self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid),
method='linear')
self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid), method='linear')
except Exception as e:
self._warn('Could not init picksgrid: {}'.format(e))
@@ -382,16 +442,6 @@ class Array_map(QtGui.QWidget):
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 = []
uncertainties = []
@@ -404,65 +454,82 @@ class Array_map(QtGui.QWidget):
longitudes.append(self.stations_dict[st_id]['longitude'])
return picks, uncertainties, latitudes, longitudes
def draw_contour_filled(self, nlevel='100'):
# self.test_gradient()
# plotting -----------------------------------------------------
def highlight_station(self, network, station, color):
stat_dict = self.stations_dict['{}.{}'.format(network, station)]
lat = stat_dict['latitude']
lon = stat_dict['longitude']
self.highlighted_stations.append(self.canvas.axes.scatter(lon, lat, s=self.pointsize, edgecolors=color,
facecolors='none', zorder=12,
transform=ccrs.PlateCarree(), label='deleted'))
def openPickDlg(self, ind):
data = self._parent.get_data().getWFData()
for index in ind:
network, station = self._station_onpick_ids[index].split('.')[:2]
pyl_mw = self._parent
try:
data = data.select(station=station)
if not data:
self._warn('No data for station {}'.format(station))
return
pickDlg = PickDlg(self._parent, parameter=self.parameter,
data=data, network=network, station=station,
picks=self._parent.get_current_event().getPick(station),
autopicks=self._parent.get_current_event().getAutopick(station),
filteroptions=self._parent.filteroptions, metadata=self.metadata,
event=pyl_mw.get_current_event())
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)
print(traceback.format_exc())
return
try:
if pickDlg.exec_():
pyl_mw.setDirty(True)
pyl_mw.update_status('picks accepted ({0})'.format(station))
pyl_mw.addPicks(station, pickDlg.getPicks(picktype='manual'), type='manual')
pyl_mw.addPicks(station, pickDlg.getPicks(picktype='auto'), type='auto')
if self.auto_refresh_box.isChecked():
self._refresh_drawings()
else:
self.highlight_station(network, station, color='yellow')
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)
print(traceback.format_exc())
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.contour(self.longrid, self.latgrid, self.picksgrid_active, levels,
linewidths=self.linewidth, latlon=True, zorder=8, alpha=0.7,
cmap=self.get_colormap())
self.contourf = self.canvas.axes.contourf(self.longrid, self.latgrid, self.picksgrid_active, levels,
linewidths=self.linewidth * 5, transform=ccrs.PlateCarree(),
alpha=0.4, zorder=8, cmap=self.get_colormap())
def get_colormap(self):
return plt.get_cmap(self.cmaps_box.currentText())
def test_gradient(self):
st_ids = self.picks_rel.keys()
x, y = np.gradient(self.picksgrid_active)
gradient_modulus = np.sqrt(x ** 2 + y ** 2)
global_mean_gradient = np.nanmean(gradient_modulus)
delta_gradient = []
for st_id in st_ids:
pick_item = self.picks_rel.pop(st_id)
self.init_picksgrid()
x, y = np.gradient(self.picksgrid_active)
gradient_modulus = np.sqrt(x ** 2 + y ** 2)
mean_gradient = np.nanmean(gradient_modulus)
dgradient = global_mean_gradient - mean_gradient
# print('station: {}, mean gradient: {}'.format(st_id, dgradient))
delta_gradient.append(dgradient)
self.picks_rel[st_id] = pick_item
global_std_gradient = np.nanstd(delta_gradient)
marked_stations = []
for st_id, dg in zip(st_ids, delta_gradient):
if abs(dg) > global_std_gradient:
marked_stations.append(st_id)
self.marked_stations = marked_stations
self.init_picksgrid()
# fig = plt.figure()
# x = list(range(len(st_ids)))
# gradients = zip(x, delta_gradient)
# gradients.sort(key=lambda a: a[1])
# plt.plot(gradients[0], gradients[1])
# global_var_gradient = np.nanvar(delta_gradient)
# plt.plot(x, delta_gradient)
# plt.axhline(global_std_gradient, color='green')
# plt.axhline(2 * global_std_gradient, color='blue')
# plt.axhline(global_var_gradient, color='red')
# plt.xticks(x, st_ids)
# plt.show()
def scatter_all_stations(self):
stations, lats, lons = self.get_st_lat_lon_for_plot()
self.sc = self.basemap.scatter(lons, lats, s=self.pointsize, facecolor='none', latlon=True, marker='.',
zorder=10, picker=True, edgecolor='0.5', label='Not Picked')
self.cid = self.canvas.mpl_connect('pick_event', self.onpick)
self.sc = self.canvas.axes.scatter(lons, lats, s=self.pointsize * 3, facecolor='none', marker='.',
zorder=10, picker=True, edgecolor='0.5', label='Not Picked',
transform=ccrs.PlateCarree())
self.cid = self.plotWidget.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=2*self.pointsize, facecolor='red',
latlon=True, zorder=11, label='Event (might be outside map region)')
self.sc_event = self.canvas.axes.scatter(lons, lats, s=5 * self.pointsize, facecolor='red', zorder=11,
label='Event (might be outside map region)', marker='*',
edgecolors='black',
transform=ccrs.PlateCarree())
def scatter_picked_stations(self):
picks, uncertainties, lats, lons = self.get_picks_lat_lon()
@@ -471,24 +538,16 @@ class Array_map(QtGui.QWidget):
phase = self.comboBox_phase.currentText()
timeerrors = self.parameter['timeerrors{}'.format(phase)]
sizes = np.array([self.pointsize * ((5. - get_quality_class(uncertainty, timeerrors)))
sizes = np.array([self.pointsize * (5. - get_quality_class(uncertainty, timeerrors))
for uncertainty in uncertainties])
cmap = self.get_colormap()
# 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=sizes, edgecolors='white', cmap=cmap,
c=picks[0], latlon=True, zorder=11)
if len(lons) == 2 and len(lats) == 2:
self.sc_picked = self.basemap.scatter(lons[1], lats[1], s=sizes, edgecolors='white', cmap=cmap,
c=picks[1], latlon=True, zorder=11)
if len(lons) > 2 and len(lats) > 2:
self.sc_picked = self.basemap.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap,
c=picks, latlon=True, zorder=11, label='Picked')
self.sc_picked = self.canvas.axes.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap,
c=picks, zorder=11, label='Picked', transform=ccrs.PlateCarree())
def annotate_ax(self):
self.annotations = []
stations, xs, ys = self.get_st_x_y_for_plot()
stations, ys, xs = self.get_st_lat_lon_for_plot()
# MP MP testing station highlighting if they have high impact on mean gradient of color map
# if self.picks_rel:
# self.test_gradient()
@@ -501,16 +560,21 @@ class Array_map(QtGui.QWidget):
color = 'lightgrey'
if st in self.marked_stations:
color = 'red'
self.annotations.append(self.main_ax.annotate(' %s' % st, xy=(x, y),
fontsize=self.pointsize/4., fontweight='semibold',
color=color, zorder=14))
self.legend = self.main_ax.legend(loc=1)
self.legend.get_frame().set_facecolor((1, 1, 1, 0.75))
self.annotations.append(
self.canvas.axes.annotate(' %s' % st, xy=(x + 0.003, y + 0.003), fontsize=self.pointsize / 4.,
fontweight='semibold', color=color, alpha=0.8,
transform=ccrs.PlateCarree(), zorder=14,
path_effects=[PathEffects.withStroke(
linewidth=self.pointsize / 15., foreground='k')]))
self.legend = self.canvas.axes.legend(loc=1, framealpha=1)
self.legend.set_zorder(100)
self.legend.get_frame().set_facecolor((1, 1, 1, 0.95))
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)
self.cbax_bg = inset_axes(self.canvas.axes, width="6%", height="75%", loc=5)
cbax = inset_axes(self.canvas.axes, width='2%', height='70%', loc=5)
cbar = self.canvas.axes.figure.colorbar(self.sc_picked, cax=cbax)
cbar.set_label(label)
cbax.yaxis.tick_left()
cbax.yaxis.set_label_position('left')
@@ -521,6 +585,7 @@ class Array_map(QtGui.QWidget):
self.cbax_bg.patch.set_facecolor((1, 1, 1, 0.75))
return cbar
# handle drawings -----------------------------------------------------
def refresh_drawings(self, picks=None, autopicks=None):
self.picks_dict = picks
self.autopicks_dict = autopicks
@@ -560,7 +625,7 @@ class Array_map(QtGui.QWidget):
self.comboBox_phase.setEnabled(False)
if self.annotate:
self.annotate_ax()
self.canvas.draw()
self.plotWidget.draw_idle()
def remove_drawings(self):
self.remove_annotations()
@@ -584,7 +649,7 @@ class Array_map(QtGui.QWidget):
self.remove_contourf()
del self.contourf
if hasattr(self, 'cid'):
self.canvas.mpl_disconnect(self.cid)
self.plotWidget.mpl_disconnect(self.cid)
del self.cid
try:
self.sc.remove()
@@ -594,7 +659,7 @@ class Array_map(QtGui.QWidget):
self.legend.remove()
except Exception as e:
print('Warning: could not remove legend. Reason: {}'.format(e))
self.canvas.draw()
self.plotWidget.draw_idle()
def remove_contourf(self):
for item in self.contourf.collections:
@@ -605,34 +670,16 @@ class Array_map(QtGui.QWidget):
annotation.remove()
self.annotations = []
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 saveFigure(self):
if self.canvas.fig:
fd = QtWidgets.QFileDialog()
fname, filter = fd.getSaveFileName(self.parent(), filter='Images (*.png *.svg *.jpg)')
if not fname:
return
if not any([fname.endswith(item) for item in ['.png', '.svg', '.jpg']]):
fname += '.png'
self.canvas.fig.savefig(fname)
def _warn(self, message):
self.qmb = QtGui.QMessageBox(QtGui.QMessageBox.Icon.Warning,
'Warning', message)
self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Warning, 'Warning', message)
self.qmb.show()
+2 -1
View File
@@ -2,12 +2,13 @@
# -*- coding: utf-8 -*-
try:
# noinspection PyUnresolvedReferences
from urllib2 import urlopen
except:
from urllib.request import urlopen
def checkurl(url='https://ariadne.geophysik.ruhr-uni-bochum.de/trac/PyLoT/'):
def checkurl(url='https://git.geophysik.ruhr-uni-bochum.de/marcel/pylot/'):
"""
check if URL is available
:param url: url
+18 -15
View File
@@ -2,9 +2,10 @@
# -*- coding: utf-8 -*-
import glob
import numpy as np
import os
import sys
import numpy as np
from obspy import UTCDateTime, read_inventory, read
from obspy.io.xseed import Parser
@@ -46,7 +47,7 @@ class Metadata(object):
def __repr__(self):
return self.__str__()
def add_inventory(self, path_to_inventory, obspy_dmt_inv = False):
def add_inventory(self, path_to_inventory, obspy_dmt_inv=False):
"""
Add path to list of inventories.
:param path_to_inventory: Path to a folder
@@ -83,10 +84,10 @@ class Metadata(object):
print('Path {} not in inventories list.'.format(path_to_inventory))
return
self.inventories.remove(path_to_inventory)
for filename in self.inventory_files.keys():
for filename in list(self.inventory_files.keys()):
if filename.startswith(path_to_inventory):
del (self.inventory_files[filename])
for seed_id in self.seed_ids.keys():
for seed_id in list(self.seed_ids.keys()):
if self.seed_ids[seed_id].startswith(path_to_inventory):
del (self.seed_ids[seed_id])
# have to clean self.stations_dict as well
@@ -211,6 +212,7 @@ class Metadata(object):
self.stations_dict[st_id] = {'latitude': station[0].latitude,
'longitude': station[0].longitude,
'elevation': station[0].elevation}
read_stat = {'xml': stat_info_from_inventory,
'dless': stat_info_from_parser}
@@ -351,25 +353,25 @@ def check_time(datetime):
:type datetime: list
:return: returns True if Values are in supposed range, returns False otherwise
>>> check_time([1999, 01, 01, 23, 59, 59, 999000])
>>> check_time([1999, 1, 1, 23, 59, 59, 999000])
True
>>> check_time([1999, 01, 01, 23, 59, 60, 999000])
>>> check_time([1999, 1, 1, 23, 59, 60, 999000])
False
>>> check_time([1999, 01, 01, 23, 59, 59, 1000000])
>>> check_time([1999, 1, 1, 23, 59, 59, 1000000])
False
>>> check_time([1999, 01, 01, 23, 60, 59, 999000])
>>> check_time([1999, 1, 1, 23, 60, 59, 999000])
False
>>> check_time([1999, 01, 01, 23, 60, 59, 999000])
>>> check_time([1999, 1, 1, 23, 60, 59, 999000])
False
>>> check_time([1999, 01, 01, 24, 59, 59, 999000])
>>> check_time([1999, 1, 1, 24, 59, 59, 999000])
False
>>> check_time([1999, 01, 31, 23, 59, 59, 999000])
>>> check_time([1999, 1, 31, 23, 59, 59, 999000])
True
>>> check_time([1999, 02, 30, 23, 59, 59, 999000])
>>> check_time([1999, 2, 30, 23, 59, 59, 999000])
False
>>> check_time([1999, 02, 29, 23, 59, 59, 999000])
>>> check_time([1999, 2, 29, 23, 59, 59, 999000])
False
>>> check_time([2000, 02, 29, 23, 59, 59, 999000])
>>> check_time([2000, 2, 29, 23, 59, 59, 999000])
True
>>> check_time([2000, 13, 29, 23, 59, 59, 999000])
False
@@ -381,6 +383,7 @@ def check_time(datetime):
return False
# TODO: change root to datapath
def get_file_list(root_dir):
"""
Function uses a directorie to get all the *.gse files from it.
@@ -444,7 +447,7 @@ def evt_head_check(root_dir, out_dir=None):
"""
if not out_dir:
print('WARNING files are going to be overwritten!')
inp = str(raw_input('Continue? [y/N]'))
inp = str(input('Continue? [y/N]'))
if not inp == 'y':
sys.exit()
filelist = get_file_list(root_dir)
+10 -8
View File
@@ -2,6 +2,7 @@
# -*- coding: utf-8 -*-
import os
from obspy import UTCDateTime
from obspy.core.event import Event as ObsPyEvent
from obspy.core.event import Origin, ResourceIdentifier
@@ -21,12 +22,13 @@ class Event(ObsPyEvent):
:param path: path to event directory
:type path: str
"""
# TODO: remove rootpath and database
self.pylot_id = path.split('/')[-1]
# initialize super class
super(Event, self).__init__(resource_id=ResourceIdentifier('smi:local/' + self.pylot_id))
self.path = path
self.database = path.split('/')[-2]
self.datapath = path.split('/')[-3]
self.datapath = os.path.split(path)[0] # path.split('/')[-3]
self.rootpath = '/' + os.path.join(*path.split('/')[:-3])
self.pylot_autopicks = {}
self.pylot_picks = {}
@@ -73,7 +75,7 @@ class Event(ObsPyEvent):
text = lines[0]
self.addNotes(text)
try:
datetime = UTCDateTime(path.split('/')[-1])
datetime = UTCDateTime(self.path.split('/')[-1])
origin = Origin(resource_id=self.resource_id, time=datetime, latitude=0, longitude=0, depth=0)
self.origins.append(origin)
except:
@@ -296,13 +298,13 @@ class Event(ObsPyEvent):
:rtype: None
"""
try:
import cPickle
import pickle
except ImportError:
import _pickle as cPickle
import _pickle as pickle
try:
outfile = open(filename, 'wb')
cPickle.dump(self, outfile, -1)
pickle.dump(self, outfile, -1)
self.dirty = False
except Exception as e:
print('Could not pickle PyLoT event. Reason: {}'.format(e))
@@ -317,11 +319,11 @@ class Event(ObsPyEvent):
:rtype: Event
"""
try:
import cPickle
import pickle
except ImportError:
import _pickle as cPickle
import _pickle as pickle
infile = open(filename, 'rb')
event = cPickle.load(infile)
event = pickle.load(infile)
event.dirty = False
print('Loaded %s' % filename)
return event
+38 -16
View File
@@ -3,24 +3,37 @@
# small script that creates array maps for each event within a previously generated PyLoT project
import os
num_thread = "16"
os.environ["OMP_NUM_THREADS"] = num_thread
os.environ["OPENBLAS_NUM_THREADS"] = num_thread
os.environ["MKL_NUM_THREADS"] = num_thread
os.environ["VECLIB_MAXIMUM_THREADS"] = num_thread
os.environ["NUMEXPR_NUM_THREADS"] = num_thread
os.environ["NUMEXPR_MAX_THREADS"] = num_thread
import multiprocessing
import sys
import glob
import matplotlib
matplotlib.use('Qt5Agg')
sys.path.append(os.path.join('/'.join(sys.argv[0].split('/')[:-1]), '../../..'))
from PyLoT import Project
from pylot.core.util.dataprocessing import Metadata
from pylot.core.util.array_map import Array_map
import matplotlib.pyplot as plt
import argparse
def main(project_file_path, manual=False, auto=True, file_format='png', f_ext='', ncores=None):
project = Project.load(project_file_path)
nEvents = len(project.eventlist)
input_list = []
print('\n')
for index, event in enumerate(project.eventlist):
# MP MP TESTING +++
#if not eventdir.endswith('20170908_044946.a'):
# continue
# MP MP ----
kwargs = dict(project=project, event=event, nEvents=nEvents, index=index, manual=manual, auto=auto,
file_format=file_format, f_ext=f_ext)
input_list.append(kwargs)
@@ -30,14 +43,20 @@ def main(project_file_path, manual=False, auto=True, file_format='png', f_ext=''
array_map_worker(item)
else:
pool = multiprocessing.Pool(ncores)
result = pool.map(array_map_worker, input_list)
pool.map(array_map_worker, input_list)
pool.close()
pool.join()
def array_map_worker(input_dict):
event = input_dict['event']
eventdir = event.path
print('Working on event: {} ({}/{})'.format(eventdir, input_dict['index'] + 1, input_dict['nEvents']))
xml_picks = glob.glob(os.path.join(eventdir, f'*{input_dict["f_ext"]}.xml'))
if not len(xml_picks):
print('Event {} does not have any picks associated with event file extension {}'.format(eventdir,
input_dict['f_ext']))
return
# check for picks
manualpicks = event.getPicks()
autopicks = event.getAutopicks()
@@ -51,11 +70,12 @@ def array_map_worker(input_dict):
continue
if not metadata:
metadata = Metadata(inventory=metadata_path, verbosity=0)
# create figure to plot on
fig = plt.figure(figsize=(16,9))
fig, ax = plt.subplots(figsize=(15, 9))
# create array map object
map = Array_map(None, metadata, parameter=input_dict['project'].parameter, figure=fig,
width=2.13e6, height=1.2e6, pointsize=15., linewidth=1.0)
map = Array_map(None, metadata, parameter=input_dict['project'].parameter, axes=ax,
width=2.13e6, height=1.2e6, pointsize=25., linewidth=1.0)
# set combobox to auto/manual to plot correct pick type
map.comboBox_am.setCurrentIndex(map.comboBox_am.findText(pick_type))
# add picks to map and save file
@@ -65,11 +85,13 @@ def array_map_worker(input_dict):
fig.savefig(fpath_out, dpi=300.)
print('Wrote file: {}'.format(fpath_out))
if __name__ == '__main__':
dataroot = '/home/marcel'
infiles=['alparray_all_events_0.03-0.1_mantle_correlated_v3.plp']
for infile in infiles:
main(os.path.join(dataroot, infile), f_ext='_correlated_0.1Hz', ncores=10)
#main('E:\Shared\AlpArray\\test_aa.plp', f_ext='_correlated_0.5Hz', ncores=1)
#main('/home/marcel/alparray_m6.5-6.9_mantle_correlated_v3.plp', f_ext='_correlated_0.5Hz')
if __name__ == '__main__':
cl = argparse.ArgumentParser()
cl.add_argument('--dataroot', help='Directory containing the PyLoT .plp file', type=str)
cl.add_argument('--infiles', help='.plp files to use', nargs='+')
cl.add_argument('--ncores', hepl='Specify number of parallel processes', type=int, default=1)
args = cl.parse_args()
for infile in args.infiles:
main(os.path.join(args.dataroot, infile), f_ext='_correlated_0.03-0.1', ncores=args.ncores)
+4 -6
View File
@@ -7,11 +7,10 @@ try:
except Exception as e:
print('Warning: Could not import module pyqtgraph.')
try:
from PySide import QtCore
from PySide2 import QtCore
except Exception as e:
print('Warning: Could not import module QtCore.')
from pylot.core.util.utils import pick_color
@@ -49,11 +48,11 @@ def which(program, parameter):
:rtype: str
"""
try:
from PySide.QtCore import QSettings
from PySide2.QtCore import QSettings
settings = QSettings()
for key in settings.allKeys():
if 'binPath' in key:
os.environ['PATH'] += ':{0}'.format(settings.value(key))
os.environ['PATH'] += ':{0}'.format(settings.value(key))
nllocpath = ":" + parameter.get('nllocbin')
os.environ['PATH'] += nllocpath
except Exception as e:
@@ -73,7 +72,7 @@ def which(program, parameter):
return program
else:
for path in os.environ["PATH"].split(os.pathsep):
exe_file = os.path.join(path, program)
exe_file = os.path.join(path, program)
for candidate in ext_candidates(exe_file):
if is_exe(candidate):
return candidate
@@ -101,4 +100,3 @@ def make_pen(picktype, phase, key, quality):
linestyle, width = pick_linestyle_pg(picktype, key)
pen = pg.mkPen(rgba, width=width, style=linestyle)
return pen
+15 -6
View File
@@ -2,6 +2,7 @@
# -*- coding: utf-8 -*-
import os
from obspy import UTCDateTime
@@ -34,17 +35,25 @@ def qml_from_obspyDMT(path):
if not os.path.exists(path):
return IOError('Could not find Event at {}'.format(path))
infile = open(path, 'rb')
event_dmt = pickle.load(infile)#, fix_imports=True)
with open(path, 'rb') as infile:
event_dmt = pickle.load(infile) # , fix_imports=True)
event_dmt['origin_id'].id = str(event_dmt['origin_id'].id)
ev = Event(resource_id=event_dmt['event_id'])
#small bugfix "unhashable type: 'newstr' "
# small bugfix "unhashable type: 'newstr' "
event_dmt['origin_id'].id = str(event_dmt['origin_id'].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 = 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
+2 -1
View File
@@ -1,8 +1,9 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import warnings
import numpy as np
from obspy import UTCDateTime
from pylot.core.util.utils import fit_curve, clims
+28 -6
View File
@@ -4,8 +4,8 @@ import os
import sys
import traceback
from PySide.QtCore import QThread, Signal, Qt, Slot, QRunnable, QObject
from PySide.QtGui import QDialog, QProgressBar, QLabel, QHBoxLayout
from PySide2.QtCore import QThread, Signal, Qt, Slot, QRunnable, QObject
from PySide2.QtWidgets import QDialog, QProgressBar, QLabel, QHBoxLayout, QPushButton
class Thread(QThread):
@@ -22,9 +22,11 @@ class Thread(QThread):
self.abortButton = abortButton
self.finished.connect(self.hideProgressbar)
self.showProgressbar()
self.old_stdout = None
def run(self):
if self.redirect_stdout:
self.old_stdout = sys.stdout
sys.stdout = self
try:
if self.arg is not None:
@@ -39,7 +41,7 @@ class Thread(QThread):
exctype, value = sys.exc_info()[:2]
self._executedErrorInfo = '{} {} {}'. \
format(exctype, value, traceback.format_exc())
sys.stdout = sys.__stdout__
sys.stdout = self.old_stdout
def showProgressbar(self):
if self.progressText:
@@ -49,7 +51,23 @@ class Thread(QThread):
# self.pb_widget.setWindowFlags(Qt.SplashScreen)
# self.pb_widget.setModal(True)
self.pb_widget.label.setText(self.progressText)
# generate widget if not given in init
if not self.pb_widget:
self.pb_widget = QDialog(self.parent())
self.pb_widget.setWindowFlags(Qt.SplashScreen)
self.pb_widget.setModal(True)
# add button
delete_button = QPushButton('X')
delete_button.clicked.connect(self.exit)
hl = QHBoxLayout()
pb = QProgressBar()
pb.setRange(0, 0)
hl.addWidget(pb)
hl.addWidget(QLabel(self.progressText))
if self.abortButton:
hl.addWidget(delete_button)
self.pb_widget.setLayout(hl)
self.pb_widget.show()
def hideProgressbar(self):
@@ -80,10 +98,12 @@ class Worker(QRunnable):
self.progressText = progressText
self.pb_widget = pb_widget
self.redirect_stdout = redirect_stdout
self.old_stdout = None
@Slot()
def run(self):
if self.redirect_stdout:
self.old_stdout = sys.stdout
sys.stdout = self
try:
@@ -96,7 +116,7 @@ class Worker(QRunnable):
self.signals.result.emit(result)
finally:
self.signals.finished.emit('Done')
sys.stdout = sys.__stdout__
sys.stdout = self.old_stdout
def write(self, text):
self.signals.message.emit(text)
@@ -128,11 +148,13 @@ class MultiThread(QThread):
self.progressText = progressText
self.pb_widget = pb_widget
self.redirect_stdout = redirect_stdout
self.old_stdout = None
self.finished.connect(self.hideProgressbar)
self.showProgressbar()
def run(self):
if self.redirect_stdout:
self.old_stdout = sys.stdout
sys.stdout = self
try:
if not self.ncores:
@@ -148,7 +170,7 @@ class MultiThread(QThread):
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 = self.old_stdout
def showProgressbar(self):
if self.progressText:
+11 -7
View File
@@ -2,12 +2,13 @@
# -*- coding: utf-8 -*-
import hashlib
import numpy as np
import os
import platform
import re
import subprocess
import warnings
import numpy as np
from obspy import UTCDateTime, read
from obspy.core import AttribDict
from obspy.signal.rotate import rotate2zne
@@ -227,7 +228,7 @@ def findComboBoxIndex(combo_box, val):
:type val: basestring
:return: index value of item with name val or 0
"""
return combo_box.findText(val) if combo_box.findText(val) is not -1 else 0
return combo_box.findText(val) if combo_box.findText(val) != -1 else 0
def find_in_list(list, str):
@@ -950,10 +951,13 @@ def check4rotated(data, metadata=None, verbosity=1):
if any(rotation_required):
t_start = full_range(wfstream)
try:
azimuts = [metadata.get_coordinates(tr_id, t_start)['azimuth'] for tr_id in trace_ids]
dips = [metadata.get_coordinates(tr_id, t_start)['dip'] for tr_id in trace_ids]
azimuts = []
dips = []
for tr_id in trace_ids:
azimuts.append(metadata.get_coordinates(tr_id, t_start)['azimuth'])
dips.append(metadata.get_coordinates(tr_id, t_start)['dip'])
except (KeyError, TypeError) as e:
print('Failed to rotate trace {}, no azimuth or dip available in metadata'.format(trace_id))
print('Failed to rotate trace {}, no azimuth or dip available in metadata'.format(tr_id))
return wfstream
if len(wfstream) < 3:
print('Failed to rotate Stream {}, not enough components available.'.format(wfstream))
@@ -964,7 +968,7 @@ def check4rotated(data, metadata=None, verbosity=1):
z, n, e = rotate2zne(wfstream[0], azimuts[0], dips[0],
wfstream[1], azimuts[1], dips[1],
wfstream[2], azimuts[2], dips[2])
print('check4rotated: rotated trace {} to ZNE'.format(trace_id))
print('check4rotated: rotated trace {} to ZNE'.format(trace_ids))
# replace old data with rotated data, change the channel code to ZNE
z_index = dips.index(min(
dips)) # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90 being vertical)
@@ -1013,7 +1017,7 @@ def scaleWFData(data, factor=None, components='all'):
:return: scaled waveform data
:rtype: `~obspy.core.stream.Stream` object
"""
if components is not 'all':
if components != 'all':
for comp in components:
if factor is None:
max_val = np.max(np.abs(data.select(component=comp)[0].data))
+1 -1
View File
@@ -35,9 +35,9 @@ from __future__ import print_function
__all__ = "get_git_version"
import inspect
# NO IMPORTS FROM PYLOT IN THIS FILE! (file gets used at installation time)
import os
import inspect
from subprocess import Popen, PIPE
# NO IMPORTS FROM PYLOT IN THIS FILE! (file gets used at installation time)
+280 -232
View File
File diff suppressed because it is too large Load Diff