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

Merge branch 'develop' into feature/documentation

Browse Source
This commit is contained in:
Darius Arnold 2018-09-14 21:46:29 +02:00
commit 2d233a7099
46 changed files with 3190 additions and 1683 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -1,3 +1,5 @@
*.pyc
*~
pylot/RELEASE-VERSION
*.idea
autopylot.sh*

964
PyLoT.py
View File

File diff suppressed because it is too large Load Diff

107
autoPyLoT.py
View File

@ -7,6 +7,7 @@ import argparse
import datetime
import glob
import os
import traceback
import pylot.core.loc.focmec as focmec
import pylot.core.loc.hash as hash
@ -22,7 +23,7 @@ from pylot.core.analysis.magnitude import MomentMagnitude, LocalMagnitude
from pylot.core.io.data import Data
from pylot.core.io.inputs import PylotParameter
from pylot.core.pick.autopick import autopickevent, iteratepicker
from pylot.core.util.dataprocessing import restitute_data, read_metadata
from pylot.core.util.dataprocessing import restitute_data, read_metadata, Metadata
from pylot.core.util.defaults import SEPARATOR
from pylot.core.util.event import Event
from pylot.core.util.structure import DATASTRUCTURE
@ -34,10 +35,11 @@ __version__ = _getVersionString()
def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, eventid=None, savepath=None,
savexml=True, station='all', iplot=0, ncores=0):
savexml=True, station='all', iplot=0, ncores=0, obspyDMT_wfpath=False):
"""
Determine phase onsets automatically utilizing the automatic picking
algorithms by Kueperkoch et al. 2010/2012.
:param obspyDMT_wfpath: if obspyDMT is used, name of data directory ("raw" or "processed")
:param input_dict:
:type input_dict:
:param parameter: PylotParameter object containing parameters used for automatic picking
@ -53,7 +55,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
:type savepath: str
:param savexml: export results in XML file if True
:type savexml: bool
:param station: list of station names or 'all' to pick all stations
:param station: choose specific station name or 'all' to pick all stations
:type station: str
:param iplot: logical variable for plotting: 0=none, 1=partial, 2=all
:type iplot: int
@ -94,7 +96,6 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
fig_dict = None
fig_dict_wadatijack = None
locflag = 1
if input_dict and isinstance(input_dict, dict):
if 'parameter' in input_dict:
parameter = input_dict['parameter']
@ -110,15 +111,15 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
eventid = input_dict['eventid']
if 'iplot' in input_dict:
iplot = input_dict['iplot']
if 'locflag' in input_dict:
locflag = input_dict['locflag']
if 'savexml' in input_dict:
savexml = input_dict['savexml']
if 'obspyDMT_wfpath' in input_dict:
obspyDMT_wfpath = input_dict['obspyDMT_wfpath']
if not parameter:
if inputfile:
parameter = PylotParameter(inputfile)
#iplot = parameter['iplot']
# iplot = parameter['iplot']
else:
infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
print('Using default input file {}'.format(infile))
@ -150,8 +151,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
datastructure.modifyFields(**dsfields)
datastructure.setExpandFields(exf)
# check if default location routine NLLoc is available
if real_None(parameter['nllocbin']) and locflag:
# check if default location routine NLLoc is available and all stations are used
if real_None(parameter['nllocbin']) and station == 'all':
locflag = 1
# get NLLoc-root path
nllocroot = parameter.get('nllocroot')
# get path to NLLoc executable
@ -167,7 +169,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
ttpat = parameter.get('ttpatter')
# pattern of NLLoc-output file
nllocoutpatter = parameter.get('outpatter')
maxnumit = 3 # maximum number of iterations for re-picking
maxnumit = 2 # maximum number of iterations for re-picking
else:
locflag = 0
print(" !!! ")
@ -175,6 +177,14 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
print("!!No source parameter estimation possible!!")
print(" !!! ")
wfpath_extension = ''
if obspyDMT_wfpath not in [None, False, 'False', '']:
wfpath_extension = obspyDMT_wfpath
print('Using obspyDMT structure. There will be no restitution, as pre-processed data are expected.')
if wfpath_extension != 'processed':
print('WARNING: Expecting wfpath_extension to be "processed" for'
' pre-processed data but received "{}" instead!!!'.format(wfpath_extension))
if not input_dict:
# started in production mode
datapath = datastructure.expandDataPath()
@ -192,8 +202,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
events.append(os.path.join(datapath, eventID))
else:
# autoPyLoT was initialized from GUI
events = []
events.append(eventid)
events = [eventid]
evID = os.path.split(eventid)[-1]
locflag = 2
else:
@ -219,6 +228,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
glocflag = locflag
for eventpath in events:
evID = os.path.split(eventpath)[-1]
event_datapath = os.path.join(eventpath, wfpath_extension)
fext = '.xml'
filename = os.path.join(eventpath, 'PyLoT_' + evID + fext)
try:
@ -231,7 +241,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
pylot_event = Event(eventpath) # event should be path to event directory
data.setEvtData(pylot_event)
if fnames == 'None':
data.setWFData(glob.glob(os.path.join(datapath, eventpath, '*')))
data.setWFData(glob.glob(os.path.join(datapath, event_datapath, '*')))
# the following is necessary because within
# multiple event processing no event ID is provided
# in autopylot.in
@ -262,18 +272,26 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
if not wfdat:
print('Could not find station {}. STOP!'.format(station))
return
wfdat = remove_underscores(wfdat)
#wfdat = remove_underscores(wfdat)
# trim components for each station to avoid problems with different trace starttimes for one station
wfdat = check4gaps(wfdat)
wfdat = check4doubled(wfdat)
wfdat = trim_station_components(wfdat, trim_start=True, trim_end=False)
metadata = read_metadata(parameter.get('invdir'))
# rotate stations to ZNE
wfdat = check4rotated(wfdat, metadata)
if not wfpath_extension:
metadata = Metadata(parameter.get('invdir'))
else:
metadata = Metadata(os.path.join(eventpath, 'resp'))
corr_dat = None
if locflag:
print("Restitute data ...")
corr_dat = restitute_data(wfdat.copy(), *metadata, ncores=ncores)
if metadata:
# rotate stations to ZNE
try:
wfdat = check4rotated(wfdat, metadata)
except Exception as e:
print('Could not rotate station {} to ZNE:\n{}'.format(wfdat[0].stats.station,
traceback.format_exc()))
if locflag:
print("Restitute data ...")
corr_dat = restitute_data(wfdat.copy(), metadata, ncores=ncores)
if not corr_dat and locflag:
locflag = 2
print('Working on event %s. Stations: %s' % (eventpath, station))
@ -299,7 +317,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
ttpat)
# locate the event
nll.locate(ctrfile, inputfile)
nll.locate(ctrfile, parameter)
# !iterative picking if traces remained unpicked or occupied with bad picks!
# get theoretical onset times for picks with weights >= 4
@ -331,7 +349,8 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
picks[stats]['P'].update(props)
evt = moment_mag.updated_event()
net_mw = moment_mag.net_magnitude()
print("Network moment magnitude: %4.1f" % net_mw.mag)
if net_mw is not None:
print("Network moment magnitude: %4.1f" % net_mw.mag)
# calculate local (Richter) magntiude
WAscaling = parameter.get('WAscaling')
magscaling = parameter.get('magscaling')
@ -347,9 +366,17 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
WAscaling[2]))
evt = local_mag.updated_event(magscaling)
net_ml = local_mag.net_magnitude(magscaling)
print("Network local magnitude: %4.1f" % net_ml.mag)
print("Network local magnitude scaled with:")
print("%f * Ml + %f" % (magscaling[0], magscaling[1]))
if net_ml:
print("Network local magnitude: %4.1f" % net_ml.mag)
if magscaling is None:
scaling = False
elif magscaling[0] != 0 and magscaling[1] != 0:
scaling = False
else:
scaling = True
if scaling:
print("Network local magnitude scaled with:")
print("%f * Ml + %f" % (magscaling[0], magscaling[1]))
else:
print("autoPyLoT: No NLLoc-location file available!")
print("No source parameter estimation possible!")
@ -379,7 +406,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
# remove actual NLLoc-location file to keep only the last
os.remove(nllocfile)
# locate the event
nll.locate(ctrfile, inputfile)
nll.locate(ctrfile, parameter)
print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
# get updated NLLoc-location file
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
@ -406,7 +433,8 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
picks[stats]['P'].update(props)
evt = moment_mag.updated_event()
net_mw = moment_mag.net_magnitude()
print("Network moment magnitude: %4.1f" % net_mw.mag)
if net_mw is not None:
print("Network moment magnitude: %4.1f" % net_mw.mag)
# calculate local (Richter) magntiude
WAscaling = parameter.get('WAscaling')
magscaling = parameter.get('magscaling')
@ -423,9 +451,17 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
WAscaling[2]))
evt = local_mag.updated_event(magscaling)
net_ml = local_mag.net_magnitude(magscaling)
print("Network local magnitude: %4.1f" % net_ml.mag)
print("Network local magnitude scaled with:")
print("%f * Ml + %f" % (magscaling[0], magscaling[1]))
if net_ml:
print("Network local magnitude: %4.1f" % net_ml.mag)
if magscaling is None:
scaling = False
elif magscaling[0] != 0 and magscaling[1] != 0:
scaling = False
else:
scaling = True
if scaling:
print("Network local magnitude scaled with:")
print("%f * Ml + %f" % (magscaling[0], magscaling[1]))
else:
print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
locflag = 9
@ -468,7 +504,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
endsplash = '''------------------------------------------\n'
-----Finished event %s!-----\n'
------------------------------------------'''.format \
(version=_getVersionString()) % evID
(version=_getVersionString()) % evID
print(endsplash)
locflag = glocflag
if locflag == 0:
@ -497,9 +533,9 @@ if __name__ == "__main__":
action='store',
help='''full path to the file containing the input
parameters for autoPyLoT''')
parser.add_argument('-p', '-P', '--iplot', type=int,
parser.add_argument('-p', '-P', '--iplot', type=int,
action='store', default=0,
help='''optional, logical variable for plotting: 0=none, 1=partial, 2=all''')
help='''optional, logical variable for plotting: 0=none, 1=partial, 2=all''')
parser.add_argument('-f', '-F', '--fnames', type=str,
action='store',
help='''optional, list of data file names''')
@ -512,9 +548,12 @@ if __name__ == "__main__":
parser.add_argument('-c', '-C', '--ncores', type=int,
action='store', default=0,
help='''optional, number of CPU cores used for parallel processing (default: all available(=0))''')
parser.add_argument('-dmt', '-DMT', '--obspy_dmt_wfpath', type=str,
action='store', default=False,
help='''optional, wftype (raw, processed) used for obspyDMT database structure''')
cla = parser.parse_args()
picks = autoPyLoT(inputfile=str(cla.inputfile), fnames=str(cla.fnames),
eventid=str(cla.eventid), savepath=str(cla.spath),
ncores=cla.ncores, iplot=int(cla.iplot))
ncores=cla.ncores, iplot=int(cla.iplot), obspyDMT_wfpath=str(cla.obspy_dmt_wfpath))

1
autopylot.sh Normal file
View File

@ -0,0 +1 @@
python ./autoPyLoT.py -i /home/marcel/.pylot/pylot_global.in -dmt processed -c 80

2
icons_rc_2.py
View File

@ -7,7 +7,7 @@
#
# WARNING! All changes made in this file will be lost!
from PyQt4 import QtCore
from PySide import QtCore
qt_resource_data = "\
\x00\x00\x9e\x04\

2
icons_rc_3.py
View File

@ -7,7 +7,7 @@
#
# WARNING! All changes made in this file will be lost!
from PyQt4 import QtCore
from PySide import QtCore
qt_resource_data = "\
\x00\x00\x9e\x04\

66
inputs/pylot_local.in
View File

@ -4,19 +4,19 @@
%Parameters are optimized for %extent data sets!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#main settings#
#rootpath# %project path
#datapath# %data path
#database# %name of data base
#eventID# %event ID for single event processing (* for all events found in database)
#invdir# %full path to inventory or dataless-seed file
/DATA/Insheim #rootpath# %project path
EVENT_DATA/LOCAL #datapath# %data path
2018.02_Insheim #database# %name of data base
e0006.038.18 #eventID# %event ID for single event processing (* for all events found in database)
/DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure# %choose data structure
True #apverbose# %choose 'True' or 'False' for terminal output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#NLLoc settings#
None #nllocbin# %path to NLLoc executable
None #nllocroot# %root of NLLoc-processing directory
None #phasefile# %name of autoPyLoT-output phase file for NLLoc
None #ctrfile# %name of autoPyLoT-output control file for NLLoc
/home/ludger/NLLOC #nllocbin# %path to NLLoc executable
/home/ludger/NLLOC/Insheim #nllocroot# %root of NLLoc-processing directory
AUTOPHASES.obs #phasefile# %name of autoPyLoT-output phase file for NLLoc
Insheim_min1d032016_auto.in #ctrfile# %name of autoPyLoT-output control file for NLLoc
ttime #ttpatter# %pattern of NLLoc ttimes from grid
AUTOLOC_nlloc #outpatter# %pattern of NLLoc-output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -27,31 +27,31 @@ AUTOLOC_nlloc #outpatter# %pattern of
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#settings local magnitude#
1.11 0.0009 -2.0 #WAscaling# %Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm] If zeros are set, original Richter magnitude is calculated!
1.0382 -0.447 #magscaling# %Scaling relation for derived local magnitude [a*Ml+b]. If zeros are set, no scaling of network magnitude is applied!
0.0 0.0 #magscaling# %Scaling relation for derived local magnitude [a*Ml+b]. If zeros are set, no scaling of network magnitude is applied!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#filter settings#
1.0 1.0 #minfreq# %Lower filter frequency [P, S]
10.0 10.0 #maxfreq# %Upper filter frequency [P, S]
2 2 #filter_order# %filter order [P, S]
2.0 2.0 #minfreq# %Lower filter frequency [P, S]
30.0 15.0 #maxfreq# %Upper filter frequency [P, S]
3 3 #filter_order# %filter order [P, S]
bandpass bandpass #filter_type# %filter type (bandpass, bandstop, lowpass, highpass) [P, S]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#common settings picker#
local #extent# %extent of array ("local", "regional" or "global")
15.0 #pstart# %start time [s] for calculating CF for P-picking
60.0 #pstop# %end time [s] for calculating CF for P-picking
-1.0 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking
7.0 #pstart# %start time [s] for calculating CF for P-picking
16.0 #pstop# %end time [s] for calculating CF for P-picking
-0.5 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking
10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking
True #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
False #use_taup# %use estimated traveltimes from TauPy for calculating windows for CF
iasp91 #taup_model# %define TauPy model for traveltime estimation
2.0 10.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
2.0 12.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
2.0 8.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
2.0 10.0 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz]
2.0 20.0 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
2.0 30.0 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
2.0 10.0 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz]
2.0 15.0 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz]
#special settings for calculating CF#
%!!Edit the following only if you know what you are doing!!%
#Z-component#
HOS #algoP# %choose algorithm for P-onset determination (HOS, ARZ, or AR3)
7.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s]
4.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s]
4 #hosorder# %for HOS-picker, order of Higher Order Statistics
2 #Parorder# %for AR-picker, order of AR process of Z-component
1.2 #tdet1z# %for AR-picker, length of AR determination window [s] for Z-component, 1st pick
@ -59,12 +59,12 @@ HOS #algoP# %choose algo
0.6 #tdet2z# %for AR-picker, length of AR determination window [s] for Z-component, 2nd pick
0.2 #tpred2z# %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick
0.001 #addnoise# %add noise to seismogram for stable AR prediction
3.0 0.1 0.5 1.0 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]
3.0 0.0 1.0 0.5 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]
3.0 #pickwinP# %for initial AIC pick, length of P-pick window [s]
6.0 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)
0.2 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s]
0.4 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s]
0.1 #tsmoothP# %for HOS/AR, take average of samples for smoothing CF [s]
0.001 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P)
0.4 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P)
1.3 #nfacP# %for HOS/AR, noise factor for noise level determination (P)
#H-components#
ARH #algoS# %choose algorithm for S-onset determination (ARH or AR3)
@ -75,7 +75,7 @@ ARH #algoS# %choose algo
4 #Sarorder# %for AR-picker, order of AR process of H-components
5.0 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
4.0 #pickwinS# %for initial AIC pick, length of S-pick window [s]
2.0 0.3 1.5 1.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]
2.0 0.2 1.5 1.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]
1.0 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s]
0.7 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S)
0.9 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S)
@ -85,16 +85,16 @@ ARH #algoS# %choose algo
2.0 #minFMSNR# %miniumum required SNR for first-motion determination
0.2 #fmpickwin# %pick window around P onset for calculating zero crossings
#quality assessment#
0.02 0.04 0.08 0.16 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P
0.04 0.08 0.16 0.32 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S
0.04 0.08 0.16 0.32 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P
0.05 0.10 0.20 0.40 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S
0.8 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected
1.1 #minAICPSNR# %below this SNR the initial P pick is rejected
1.0 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected
1.5 #minAICSSNR# %below this SNR the initial S pick is rejected
1.0 #minsiglength# %length of signal part for which amplitudes must exceed noiselevel [s]
1.0 #noisefactor# %noiselevel*noisefactor=threshold
10.0 #minpercent# %required percentage of amplitudes exceeding threshold
1.5 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude
6.0 #mdttolerance# %maximum allowed deviation of P picks from median [s]
1.1 #noisefactor# %noiselevel*noisefactor=threshold
50.0 #minpercent# %required percentage of amplitudes exceeding threshold
1.1 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude
5.0 #mdttolerance# %maximum allowed deviation of P picks from median [s]
1.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram
5.0 #jackfactor# %pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor
2.0 #jackfactor# %pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor

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

@ -6,6 +6,7 @@ Revised/extended summer 2017.
:author: Ludger Küperkoch / MAGS2 EP3 working group
"""
import matplotlib.pyplot as plt
import numpy as np
import obspy.core.event as ope
@ -16,6 +17,7 @@ from pylot.core.util.utils import common_range, fit_curve
from scipy import integrate, signal
from scipy.optimize import curve_fit
def richter_magnitude_scaling(delta):
distance = np.array([0, 10, 20, 25, 30, 35, 40, 45, 50, 60, 70, 75, 85, 90, 100, 110,
120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 230, 240, 250,
@ -121,7 +123,13 @@ class Magnitude(object):
def net_magnitude(self, magscaling=None):
if self:
if magscaling is not None and str(magscaling) is not '[0.0, 0.0]':
if magscaling is None:
scaling = False
elif magscaling[0] != 0 and magscaling[1] != 0:
scaling = False
else:
scaling = True
if scaling:
# scaling necessary
print("Scaling network magnitude ...")
mag = ope.Magnitude(
@ -140,7 +148,6 @@ class Magnitude(object):
station_count=len(self.magnitudes),
azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap)
return mag
return None
class LocalMagnitude(Magnitude):
@ -219,7 +226,7 @@ class LocalMagnitude(Magnitude):
sqH = np.sqrt(power_sum)
# get time array
th = np.arange(0, len(sqH) * dt, dt)
th = np.arange(0, st[0].stats.npts / st[0].stats.sampling_rate, st[0].stats.delta)
# get maximum peak within pick window
iwin = getsignalwin(th, t0 - stime, self.calc_win)
ii = min([iwin[len(iwin) - 1], len(th)])
@ -239,9 +246,9 @@ class LocalMagnitude(Magnitude):
ax.plot(th[iwin], sqH[iwin], 'g')
ax.plot([t0 - stime, t0 - stime], [0, max(sqH)], 'r', linewidth=2)
ax.set_title('Station %s, Channel %s, RMS Horizontal Trace, '
'WA-peak-to-peak=%6.3f mm' % (st[0].stats.station,
'WA-peak-to-peak=%6.3f mm' % (st[0].stats.station,
st[0].stats.channel,
wapp))
wapp))
ax.set_xlabel('Time [s]')
ax.set_ylabel('Displacement [mm]')
ax = fig.add_subplot(212)
@ -251,15 +258,16 @@ class LocalMagnitude(Magnitude):
ax.plot([t0 - stime, t0 - stime], [0, max(sqH)], 'r', linewidth=2)
ax.set_title('Channel %s, RMS Horizontal Trace, '
'WA-peak-to-peak=%6.3f mm' % (st[1].stats.channel,
wapp))
wapp))
ax.set_xlabel('Time [s]')
ax.set_ylabel('Displacement [mm]')
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return wapp, fig
def calc(self):
@ -303,7 +311,7 @@ class LocalMagnitude(Magnitude):
a0 = a0 * 1e03 # mm to nm (see Havskov & Ottemöller, 2010)
magnitude = ope.StationMagnitude(mag=np.log10(a0) \
+ self.wascaling[0] * np.log10(delta) + self.wascaling[1]
* delta + self.wascaling[
* delta + self.wascaling[
2])
magnitude.origin_id = self.origin_id
magnitude.waveform_id = pick.waveform_id
@ -366,7 +374,7 @@ class MomentMagnitude(Magnitude):
def calc(self):
for a in self.arrivals:
if a.phase not in 'pP':
if a.phase not in 'pP':
continue
# make sure calculating Mo only from reliable onsets
# NLLoc: time_weight = 0 => do not use onset!
@ -504,6 +512,9 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
zdat = select_for_phase(wfstream, "P")
if len(zdat) == 0:
raise IOError('No vertical component found in stream:\n{}'.format(wfstream))
dt = zdat[0].stats.delta
freq = zdat[0].stats.sampling_rate

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

@ -15,6 +15,7 @@ from pylot.core.util.event import Event
from pylot.core.util.utils import fnConstructor, full_range, remove_underscores, check4gaps, check4doubled, \
check4rotated, trim_station_components
import pylot.core.loc.velest as velest
from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
class Data(object):
@ -43,7 +44,7 @@ class Data(object):
elif isinstance(evtdata, dict):
evt = readPILOTEvent(**evtdata)
evtdata = evt
elif isinstance(evtdata, str):
elif type(evtdata) in [str, unicode]:
try:
cat = read_events(evtdata)
if len(cat) is not 1:
@ -60,6 +61,8 @@ class Data(object):
raise NotImplementedError('PILOT location information '
'read support not yet '
'implemeted.')
elif 'event.pkl' in evtdata:
evtdata = qml_from_obspyDMT(evtdata)
else:
raise e
else:
@ -72,6 +75,7 @@ class Data(object):
self.wforiginal = None
self.cuttimes = None
self.dirty = False
self.processed = None
def __str__(self):
return str(self.wfdata)
@ -281,7 +285,7 @@ class Data(object):
mstation_ext = mstation + '_'
for k in range(len(picks_copy)):
if ((picks_copy[k].waveform_id.station_code == mstation) or
(picks_copy[k].waveform_id.station_code == mstation_ext)) and \
(picks_copy[k].waveform_id.station_code == mstation_ext)) and \
(picks_copy[k].method_id == 'auto'):
del picks_copy[k]
break
@ -296,7 +300,7 @@ class Data(object):
for i in range(len(picks_copy)):
if picks_copy[i].phase_hint[0] == 'P':
if (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[0]) or \
(picks_copy[i].time_errors['uncertainty'] == None):
(picks_copy[i].time_errors['uncertainty'] is None):
print("Uncertainty exceeds or equal adjusted upper time error!")
print("Adjusted uncertainty: {}".format(upperErrors[0]))
print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty']))
@ -308,7 +312,7 @@ class Data(object):
break
if picks_copy[i].phase_hint[0] == 'S':
if (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[1]) or \
(picks_copy[i].time_errors['uncertainty'] == None):
(picks_copy[i].time_errors['uncertainty'] is None):
print("Uncertainty exceeds or equal adjusted upper time error!")
print("Adjusted uncertainty: {}".format(upperErrors[1]))
print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty']))
@ -368,7 +372,7 @@ class Data(object):
data.filter(**kwargs)
self.dirty = True
def setWFData(self, fnames, checkRotated=False, metadata=None, obspy_dmt=False):
def setWFData(self, fnames, fnames_syn=None, checkRotated=False, metadata=None):
"""
Clear current waveform data and set given waveform data
:param fnames: waveform data names to append
@ -377,30 +381,31 @@ class Data(object):
self.wfdata = Stream()
self.wforiginal = None
self.wfsyn = Stream()
wffnames = None
wffnames_syn = None
wfdir = 'processed' if 'processed' in [fname.split('/')[-1] for fname in fnames] else 'raw'
if obspy_dmt:
for fpath in fnames:
if fpath.endswith(wfdir):
wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
if 'syngine' in fpath.split('/')[-1]:
wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
else:
wffnames = fnames
if wffnames is not None:
self.appendWFData(wffnames)
if wffnames_syn is not None:
self.appendWFData(wffnames_syn, synthetic=True)
# if obspy_dmt:
# wfdir = 'raw'
# self.processed = False
# for fname in fnames:
# if fname.endswith('processed'):
# wfdir = 'processed'
# self.processed = True
# break
# for fpath in fnames:
# if fpath.endswith(wfdir):
# wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
# if 'syngine' in fpath.split('/')[-1]:
# wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)]
# else:
# wffnames = fnames
if fnames is not None:
self.appendWFData(fnames)
if fnames_syn is not None:
self.appendWFData(fnames_syn, synthetic=True)
else:
return False
# various pre-processing steps:
# remove possible underscores in station names
self.wfdata = remove_underscores(self.wfdata)
# check for gaps and doubled channels
check4gaps(self.wfdata)
check4doubled(self.wfdata)
#self.wfdata = remove_underscores(self.wfdata)
# check for stations with rotated components
if checkRotated and metadata is not None:
self.wfdata = check4rotated(self.wfdata, metadata, verbosity=0)
@ -412,7 +417,6 @@ class Data(object):
self.dirty = False
return True
def appendWFData(self, fnames, synthetic=False):
"""
Read waveform data from fnames and append it to current wf data
@ -430,7 +434,7 @@ class Data(object):
False: self.wfdata}
warnmsg = ''
for fname in fnames:
for fname in set(fnames):
try:
real_or_syn_data[synthetic] += read(fname)
except TypeError:
@ -441,7 +445,7 @@ class Data(object):
except SacIOError as se:
warnmsg += '{0}\n{1}\n'.format(fname, se)
if warnmsg:
warnmsg = 'WARNING: unable to read\n' + warnmsg
warnmsg = 'WARNING in appendWFData: unable to read waveform data\n' + warnmsg
print(warnmsg)
def getWFData(self):
@ -502,8 +506,10 @@ class Data(object):
# check for automatic picks
print("Writing phases to ObsPy-quakeml file")
for key in picks:
if not picks[key].get('P'):
continue
if picks[key]['P']['picker'] == 'auto':
print("Existing picks will be overwritten!")
print("Existing auto-picks will be overwritten in pick-dictionary!")
picks = picks_from_picksdict(picks)
break
else:
@ -729,6 +735,22 @@ class PilotDataStructure(GenericDataStructure):
self.setExpandFields(['root', 'database'])
class ObspyDMTdataStructure(GenericDataStructure):
"""
Object containing the data access information for the old PILOT data
structure.
"""
def __init__(self, **fields):
if not fields:
fields = {'database': '',
'root': ''}
GenericDataStructure.__init__(self, **fields)
self.setExpandFields(['root', 'database'])
class SeiscompDataStructure(GenericDataStructure):
"""
Dictionary containing the data access information for an SDS data archive:

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

@ -63,7 +63,7 @@ defaults = {'rootpath': {'type': str,
'ctrfile': {'type': str,
'tooltip': 'name of autoPyLoT-output control file for NLLoc',
'value': 'Insheim_min1d2015_auto.in',
'value': '',
'namestring': 'Control filename'},
'ttpatter': {'type': str,
@ -456,11 +456,11 @@ defaults = {'rootpath': {'type': str,
'namestring': 'Wadati tolerance'},
'jackfactor': {'type': float,
'tooltip': 'pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor',
'value': 5.0,
'min': 0.,
'max': np.inf,
'namestring': 'Jackknife safety factor'},
'tooltip': 'pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor',
'value': 5.0,
'min': 0.,
'max': np.inf,
'namestring': 'Jackknife safety factor'},
'WAscaling': {'type': (float, float, float),
'tooltip': 'Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm] \

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

@ -48,6 +48,7 @@ class PylotParameter(object):
self.__init_default_paras()
self.__init_subsettings()
self.__filename = fnin
self.__parameter = {}
self._verbosity = verbosity
self._parFileCont = {}
# io from parsed arguments alternatively
@ -273,8 +274,8 @@ class PylotParameter(object):
:rtype: None
"""
defaults = self.get_defaults()
for param in defaults:
self.setParamKV(param, defaults[param]['value'])
for param_name, param in defaults.items():
self.setParamKV(param_name, param['value'])
def from_file(self, fnin=None):
"""

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

@ -245,7 +245,7 @@ def picksdict_from_picks(evt):
if picker.startswith('smi:local/'):
picker = picker.split('smi:local/')[1]
except IndexError:
picker = 'manual' # MP MP TODO maybe improve statement
picker = 'manual' # MP MP TODO maybe improve statement
try:
onsets = picksdict[picker][station]
except KeyError as e:
@ -346,6 +346,7 @@ def picks_from_picksdict(picks, creation_info=None):
picks_list.append(pick)
return picks_list
def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0):
import glob
@ -499,7 +500,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
except KeyError as e:
print(e)
fm = None
if fm == None:
if fm is None:
fm = '?'
onset = arrivals[key]['P']['mpp']
year = onset.year
@ -916,9 +917,9 @@ def merge_picks(event, picks):
network = pick.waveform_id.network_code
method = pick.method_id
for p in event.picks:
if p.waveform_id.station_code == station\
and p.waveform_id.network_code == network\
and p.phase_hint == phase\
if p.waveform_id.station_code == station \
and p.waveform_id.network_code == network \
and p.phase_hint == phase \
and (str(p.method_id) in str(method)
or str(method) in str(p.method_id)):
p.time, p.time_errors, p.waveform_id.network_code, p.method_id = time, err, network, method
@ -965,22 +966,22 @@ def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
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':
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.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.waveform_id.station_code == mstation_ext)) and \
(mpick.method_id.split('/')[1] == 'auto') and \
(mpick.time_errors['uncertainty'] <= ErrorsS[3]):
del mpick
break
@ -1032,19 +1033,19 @@ def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
P0perc = 0
if len(Pw1) > 0:
P1perc = 100 / numPweights * len(Pw1)
else:
else:
P1perc = 0
if len(Pw2) > 0:
P2perc = 100 / numPweights * len(Pw2)
else:
else:
P2perc = 0
if len(Pw3) > 0:
P3perc = 100 / numPweights * len(Pw3)
else:
else:
P3perc = 0
if len(Pw4) > 0:
P4perc = 100 / numPweights * len(Pw4)
else:
else:
P4perc = 0
if len(Sw0) > 0:
S0perc = 100 / numSweights * len(Sw0)
@ -1052,19 +1053,19 @@ def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
S0perc = 0
if len(Sw1) > 0:
S1perc = 100 / numSweights * len(Sw1)
else:
else:
S1perc = 0
if len(Sw2) > 0:
S2perc = 100 / numSweights * len(Sw2)
else:
else:
S2perc = 0
if len(Sw3) > 0:
S3perc = 100 / numSweights * len(Sw3)
else:
else:
S3perc = 0
if len(Sw4) > 0:
S4perc = 100 / numSweights * len(Sw4)
else:
else:
S4perc = 0
weights = ('0', '1', '2', '3', '4')

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

@ -73,17 +73,15 @@ def modify_inputs(ctrfn, root, nllocoutn, phasefn, tttn):
nllfile.close()
def locate(fnin, infile=None):
def locate(fnin, parameter=None):
"""
takes an external program name and tries to run it
:param parameter: PyLoT Parameter object
:param fnin: external program name
:return: None
"""
if infile is None:
exe_path = which('NLLoc')
else:
exe_path = which('NLLoc', infile)
exe_path = which('NLLoc', parameter)
if exe_path is None:
raise NLLocError('NonLinLoc executable not found; check your '
'environment variables')

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

@ -8,6 +8,7 @@ function conglomerate utils.
:author: MAGS2 EP3 working group / Ludger Kueperkoch
"""
import traceback
import matplotlib.pyplot as plt
import numpy as np
@ -16,7 +17,7 @@ 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
from pylot.core.util.utils import getPatternLine, gen_Pool,\
from pylot.core.util.utils import getPatternLine, gen_Pool, \
real_Bool, identifyPhaseID
from obspy.taup import TauPyModel
@ -70,42 +71,56 @@ def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None,
for station in stations:
topick = data.select(station=station)
if iplot is None or iplot == 'None' or iplot == 0:
input_tuples.append((topick, param, apverbose, metadata, origin))
if iplot > 0:
all_onsets[station] = autopickstation(topick, param, verbose=apverbose,
iplot=iplot, fig_dict=fig_dict,
metadata=metadata, origin=origin)
input_tuples.append((topick, param, apverbose, iplot, fig_dict, metadata, origin))
if iplot > 0:
print('iPlot Flag active: NO MULTIPROCESSING possible.')
return all_onsets
ncores = 1
# rename str for ncores in case ncores == 0 (use all cores)
# rename ncores for string representation in case ncores == 0 (use all cores)
ncores_str = ncores if ncores != 0 else 'all available'
print('Autopickstation: Distribute autopicking for {} '
'stations on {} cores.'.format(len(input_tuples), ncores_str))
pool = gen_Pool(ncores)
result = pool.map(call_autopickstation, input_tuples)
pool.close()
if ncores == 1:
results = serial_picking(input_tuples)
else:
results = parallel_picking(input_tuples, ncores)
for pick in result:
if pick:
station = pick['station']
pick.pop('station')
all_onsets[station] = pick
for result, station in results:
if type(result) == dict:
all_onsets[station] = result
else:
if result is None:
result = 'Picker exited unexpectedly.'
print('Could not pick a station: {}\nReason: {}'.format(station, result))
# no Wadati/JK for single station (also valid for tuning mode)
if len(stations) == 1:
return all_onsets
# quality control
# median check and jackknife on P-onset times
jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, 1, fig_dict_wadatijack)
jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, jackfactor, iplot, fig_dict_wadatijack)
# check S-P times (Wadati)
wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, 1, fig_dict_wadatijack)
wadationsets = wadaticheck(jk_checked_onsets, wdttolerance, iplot, fig_dict_wadatijack)
return wadationsets
def serial_picking(input_tuples):
result = []
for input_tuple in input_tuples:
result.append(call_autopickstation(input_tuple))
return result
def parallel_picking(input_tuples, ncores):
pool = gen_Pool(ncores)
result = pool.imap_unordered(call_autopickstation, input_tuples)
pool.close()
return result
def call_autopickstation(input_tuple):
"""
helper function used for multiprocessing
@ -114,27 +129,16 @@ def call_autopickstation(input_tuple):
:return: dictionary containing P pick, S pick and station name
:rtype: dict
"""
wfstream, pickparam, verbose, metadata, origin = input_tuple
wfstream, pickparam, verbose, iplot, fig_dict, metadata, origin = input_tuple
if fig_dict:
print('Running in interactive mode')
# multiprocessing not possible with interactive plotting
return autopickstation(wfstream, pickparam, verbose, iplot=0, metadata=metadata, origin=origin)
def get_source_coords(parser, station_id):
"""
retrieves station coordinates from metadata
:param parser: Parser object containing metadata read from inventory file
:type parser: ~obspy.io.xseed.parser.Parser
:param station_id: station id of which the coordinates should be retrieved
:type station_id: str
:return: dictionary containing 'latitude', 'longitude', 'elevation' and 'local_depth' of station
:rtype: dict
"""
station_coords = None
try:
station_coords = parser.get_coordinates(station_id)
return autopickstation(wfstream, pickparam, verbose, fig_dict=fig_dict, iplot=iplot, metadata=metadata,
origin=origin)
except Exception as e:
print('Could not get source coordinates for station {}: {}'.format(station_id, e))
return station_coords
tbe = traceback.format_exc()
return tbe, wfstream[0].stats.station
def autopickstation(wfstream, pickparam, verbose=False,
@ -179,6 +183,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
nfacP = pickparam.get('nfacP')
tpred1z = pickparam.get('tpred1z')
tdet1z = pickparam.get('tdet1z')
tpred2z = pickparam.get('tpred2z')
tdet2z = pickparam.get('tdet2z')
Parorder = pickparam.get('Parorder')
addnoise = pickparam.get('addnoise')
Precalcwin = pickparam.get('Precalcwin')
@ -248,17 +254,23 @@ def autopickstation(wfstream, pickparam, verbose=False,
# split components
zdat = wfstream.select(component="Z")
if len(zdat) == 0: # check for other components
print('HIT: 3')
zdat = wfstream.select(component="3")
edat = wfstream.select(component="E")
if len(edat) == 0: # check for other components
edat = wfstream.select(component="2")
print('HIT: 2')
ndat = wfstream.select(component="N")
if len(ndat) == 0: # check for other components
ndat = wfstream.select(component="1")
print('HIT: 1')
picks = {}
station = wfstream[0].stats.station
if not zdat:
print('No z-component found for station {}. STOP'.format(wfstream[0].stats.station))
return
print('No z-component found for station {}. STOP'.format(station))
return picks, station
if algoP == 'HOS' or algoP == 'ARZ' and zdat is not None:
msg = '##################################################\nautopickstation:' \
@ -283,12 +295,11 @@ def autopickstation(wfstream, pickparam, verbose=False,
if use_taup is True:
Lc = np.inf
print('autopickstation: use_taup flag active.')
if not metadata[1]:
if not metadata:
print('Warning: Could not use TauPy to estimate onsets as there are no metadata given.')
else:
station_id = wfstream[0].get_id()
parser = metadata[1]
station_coords = get_source_coords(parser, station_id)
station_coords = metadata.get_coordinates(station_id, time=wfstream[0].stats.starttime)
if station_coords and origin:
source_origin = origin[0]
model = TauPyModel(taup_model)
@ -320,7 +331,7 @@ def autopickstation(wfstream, pickparam, verbose=False,
# make sure pstart and pstop are inside zdat[0]
pstart = max(pstart, 0)
pstop = min(pstop, len(zdat[0])*zdat[0].stats.delta)
pstop = min(pstop, len(zdat[0]) * zdat[0].stats.delta)
if not use_taup is True or origin:
Lc = pstop - pstart
@ -328,10 +339,10 @@ def autopickstation(wfstream, pickparam, verbose=False,
Lwf = zdat[0].stats.endtime - zdat[0].stats.starttime
if not Lwf > 0:
print('autopickstation: empty trace! Return!')
return
return picks, station
Ldiff = Lwf - abs(Lc)
if Ldiff < 0 or pstop <= pstart:
if Ldiff <= 0 or pstop <= pstart or pstop - pstart <= thosmw:
msg = 'autopickstation: Cutting times are too large for actual ' \
'waveform!\nUsing entire waveform instead!'
if verbose: print(msg)
@ -487,7 +498,7 @@ def autopickstation(wfstream, pickparam, verbose=False,
elif algoP == 'ARZ':
# calculate ARZ-CF using subclass ARZcf of class
# CharcteristicFunction
cf2 = ARZcf(z_copy, cuttimes2, tpred1z, Parorder, tdet1z,
cf2 = ARZcf(z_copy, cuttimes2, tpred2z, Parorder, tdet2z,
addnoise) # instance of ARZcf
##############################################################
# get refined onset time from CF2 using class Picker
@ -562,7 +573,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
SNRPdB,
FM)
print(msg)
msg = 'autopickstation: Refined P-Pick: {} s | P-Error: {} s'.format(mpickP, Perror)
msg = 'autopickstation: Refined P-Pick: {} s | P-Error: {} s'.format(zdat[0].stats.starttime \
+ mpickP, Perror)
print(msg)
Sflag = 1
@ -596,7 +608,7 @@ def autopickstation(wfstream, pickparam, verbose=False,
ndat = edat
pickSonset = (edat is not None and ndat is not None and len(edat) > 0 and len(
ndat) > 0 and Pweight < 4)
ndat) > 0 and Pweight < 4)
if pickSonset:
# determine time window for calculating CF after P onset
@ -604,8 +616,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
round(max([mpickP + sstart, 0])), # MP MP relative time axis
round(min([
mpickP + sstop,
edat[0].stats.endtime-edat[0].stats.starttime,
ndat[0].stats.endtime-ndat[0].stats.starttime
edat[0].stats.endtime - edat[0].stats.starttime,
ndat[0].stats.endtime - ndat[0].stats.starttime
]))
]
@ -700,7 +712,7 @@ def autopickstation(wfstream, pickparam, verbose=False,
if not slope:
slope = 0
if (slope >= minAICSslope and
aicarhpick.getSNR() >= minAICSSNR and aicarhpick.getpick() is not None):
aicarhpick.getSNR() >= minAICSSNR and aicarhpick.getpick() is not None):
aicSflag = 1
msg = 'AIC S-pick passes quality control: Slope: {0} counts/s, ' \
'SNR: {1}\nGo on with refined picking ...\n' \
@ -841,7 +853,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
lpickS = lpick[ipick]
Serror = pickerr[ipick]
msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(mpickS, Serror)
msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(hdat[0].stats.starttime \
+ mpickS, Serror)
print(msg)
# get SNR
@ -852,7 +865,7 @@ def autopickstation(wfstream, pickparam, verbose=False,
Sweight = 0
elif timeerrorsS[0] < Serror <= timeerrorsS[1]:
Sweight = 1
elif Perror > timeerrorsS[1] and Serror <= timeerrorsS[2]:
elif timeerrorsS[1] < Serror <= timeerrorsS[2]:
Sweight = 2
elif timeerrorsS[2] < Serror <= timeerrorsS[3]:
Sweight = 3
@ -887,7 +900,7 @@ def autopickstation(wfstream, pickparam, verbose=False,
# re-create stream object including both horizontal components
hdat = edat.copy()
hdat += ndat
else:
print('autopickstation: No horizontal component data available or '
'bad P onset, skipping S picking!')
@ -958,115 +971,114 @@ def autopickstation(wfstream, pickparam, verbose=False,
ax1.set_ylim([-1.5, 1.5])
ax1.set_ylabel('Normalized Counts')
# fig.suptitle(tr_filt.stats.starttime)
try:
len(edat[0])
except:
edat = ndat
try:
len(ndat[0])
except:
ndat = edat
if len(edat[0]) > 1 and len(ndat[0]) > 1 and Sflag == 1:
# plot horizontal traces
ax2 = fig.add_subplot(3, 1, 2, sharex=ax1)
th1data = np.arange(0,
trH1_filt.stats.npts /
trH1_filt.stats.sampling_rate,
trH1_filt.stats.delta)
# check equal length of arrays, sometimes they are different!?
wfldiff = len(trH1_filt.data) - len(th1data)
if wfldiff < 0:
th1data = th1data[0:len(th1data) - abs(wfldiff)]
ax2.plot(th1data, trH1_filt.data / max(trH1_filt.data), color=linecolor, linewidth=0.7, label='Data')
if Pweight < 4:
ax2.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b', label='CF1')
if aicSflag == 1 and Sweight < 4:
ax2.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(arhcf2.getCF()), 'm', label='CF2')
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')
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([lpickS, lpickS], [-1.1, 1.1], 'g--', label='lpp')
ax2.plot([epickS, epickS], [-1.1, 1.1], 'g--', label='epp')
ax2.set_title('%s, S Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f' % (
trH1_filt.stats.channel,
Sweight, SNRS, SNRSdB))
else:
ax2.set_title('%s, S Weight=%d, SNR=None, SNRdB=None' % (
trH1_filt.stats.channel, Sweight))
ax2.legend(loc=1)
ax2.set_yticks([])
ax2.set_ylim([-1.5, 1.5])
ax2.set_ylabel('Normalized Counts')
# fig.suptitle(trH1_filt.stats.starttime)
# only continue if one horizontal stream exists
if (ndat or edat) and Sflag == 1:
# mirror components in case one does not exist
if not edat:
edat = ndat
if not ndat:
ndat = edat
if len(edat[0]) > 1 and len(ndat[0]) > 1:
# plot horizontal traces
ax2 = fig.add_subplot(3, 1, 2, sharex=ax1)
th1data = np.arange(0,
trH1_filt.stats.npts /
trH1_filt.stats.sampling_rate,
trH1_filt.stats.delta)
# check equal length of arrays, sometimes they are different!?
wfldiff = len(trH1_filt.data) - len(th1data)
if wfldiff < 0:
th1data = th1data[0:len(th1data) - abs(wfldiff)]
ax2.plot(th1data, trH1_filt.data / max(trH1_filt.data), color=linecolor, linewidth=0.7, label='Data')
if Pweight < 4:
ax2.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b', label='CF1')
if aicSflag == 1 and Sweight < 4:
ax2.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(arhcf2.getCF()), 'm', label='CF2')
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')
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([lpickS, lpickS], [-1.1, 1.1], 'g--', label='lpp')
ax2.plot([epickS, epickS], [-1.1, 1.1], 'g--', label='epp')
ax2.set_title('%s, S Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f' % (
trH1_filt.stats.channel,
Sweight, SNRS, SNRSdB))
else:
ax2.set_title('%s, S Weight=%d, SNR=None, SNRdB=None' % (
trH1_filt.stats.channel, Sweight))
ax2.legend(loc=1)
ax2.set_yticks([])
ax2.set_ylim([-1.5, 1.5])
ax2.set_ylabel('Normalized Counts')
# fig.suptitle(trH1_filt.stats.starttime)
ax3 = fig.add_subplot(3, 1, 3, sharex=ax1)
th2data = np.arange(0,
trH2_filt.stats.npts /
trH2_filt.stats.sampling_rate,
trH2_filt.stats.delta)
# check equal length of arrays, sometimes they are different!?
wfldiff = len(trH2_filt.data) - len(th2data)
if wfldiff < 0:
th2data = th2data[0:len(th2data) - abs(wfldiff)]
ax3.plot(th2data, trH2_filt.data / max(trH2_filt.data), color=linecolor, linewidth=0.7, label='Data')
if Pweight < 4:
p22, = ax3.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b', label='CF1')
if aicSflag == 1:
ax3.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(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')
ax3.plot([refSpick.getpick(), refSpick.getpick()],
[-1.3, 1.3], 'g', linewidth=2, label='Final S Pick')
ax3.plot(
[refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[1.3, 1.3], 'g', linewidth=2)
ax3.plot(
[refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[-1.3, -1.3], 'g', linewidth=2)
ax3.plot([lpickS, lpickS], [-1.1, 1.1], 'g--', label='lpp')
ax3.plot([epickS, epickS], [-1.1, 1.1], 'g--', label='epp')
ax3.legend(loc=1)
ax3.set_yticks([])
ax3.set_ylim([-1.5, 1.5])
ax3.set_xlabel('Time [s] after %s' % tr_filt.stats.starttime)
ax3.set_ylabel('Normalized Counts')
ax3.set_title(trH2_filt.stats.channel)
if plt_flag == 1:
fig.show()
try:
input()
except SyntaxError:
pass
plt.close(fig)
ax3 = fig.add_subplot(3, 1, 3, sharex=ax1)
th2data = np.arange(0,
trH2_filt.stats.npts /
trH2_filt.stats.sampling_rate,
trH2_filt.stats.delta)
# check equal length of arrays, sometimes they are different!?
wfldiff = len(trH2_filt.data) - len(th2data)
if wfldiff < 0:
th2data = th2data[0:len(th2data) - abs(wfldiff)]
ax3.plot(th2data, trH2_filt.data / max(trH2_filt.data), color=linecolor, linewidth=0.7, label='Data')
if Pweight < 4:
p22, = ax3.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b', label='CF1')
if aicSflag == 1:
ax3.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(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')
ax3.plot([refSpick.getpick(), refSpick.getpick()],
[-1.3, 1.3], 'g', linewidth=2, label='Final S Pick')
ax3.plot(
[refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[1.3, 1.3], 'g', linewidth=2)
ax3.plot(
[refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[-1.3, -1.3], 'g', linewidth=2)
ax3.plot([lpickS, lpickS], [-1.1, 1.1], 'g--', label='lpp')
ax3.plot([epickS, epickS], [-1.1, 1.1], 'g--', label='epp')
ax3.legend(loc=1)
ax3.set_yticks([])
ax3.set_ylim([-1.5, 1.5])
ax3.set_xlabel('Time [s] after %s' % tr_filt.stats.starttime)
ax3.set_ylabel('Normalized Counts')
ax3.set_title(trH2_filt.stats.channel)
if plt_flag == 1:
fig.show()
try:
input()
except SyntaxError:
pass
plt.close(fig)
##########################################################################
# calculate "real" onset times
if lpickP is not None and lpickP == mpickP:
@ -1084,12 +1096,22 @@ def autopickstation(wfstream, pickparam, verbose=False,
epickP = zdat[0].stats.starttime - timeerrorsP[3]
mpickP = zdat[0].stats.starttime
# create dictionary
# for P phase
ccode = zdat[0].stats.channel
ncode = zdat[0].stats.network
ppick = dict(channel=ccode, network=ncode, lpp=lpickP, epp=epickP, mpp=mpickP, spe=Perror, snr=SNRP,
snrdb=SNRPdB, weight=Pweight, fm=FM, w0=None, fc=None, Mo=None,
Mw=None, picker=picker, marked=Pmarker)
if edat:
hdat = edat[0]
elif ndat:
hdat = ndat[0]
else:
return
# no horizontal components given
picks = dict(P=ppick)
return picks, station
if lpickS is not None and lpickS == mpickS:
lpickS += hdat.stats.delta
@ -1106,21 +1128,14 @@ def autopickstation(wfstream, pickparam, verbose=False,
epickS = hdat.stats.starttime - timeerrorsS[3]
mpickS = hdat.stats.starttime
# create dictionary
# for P phase
ccode = zdat[0].stats.channel
ncode = zdat[0].stats.network
ppick = dict(channel=ccode, network=ncode, lpp=lpickP, epp=epickP, mpp=mpickP, spe=Perror, snr=SNRP,
snrdb=SNRPdB, weight=Pweight, fm=FM, w0=None, fc=None, Mo=None,
Mw=None, picker=picker, marked=Pmarker)
# add S phase
ccode = hdat.stats.channel
ncode = hdat.stats.network
spick = dict(channel=ccode, network=ncode, lpp=lpickS, epp=epickS, mpp=mpickS, spe=Serror, snr=SNRS,
snrdb=SNRSdB, weight=Sweight, fm=None, picker=picker, Ao=Ao)
# merge picks into returning dictionary
picks = dict(P=ppick, S=spick, station=zdat[0].stats.station)
return picks
picks = dict(P=ppick, S=spick)
return picks, station
def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter, fig_dict=None):
@ -1191,11 +1206,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" % (
@ -1205,7 +1220,7 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter, fig_dict=None):
print("zfac: %f => %f" % (zfac_old, pickparameter.get('zfac')))
# repick station
newpicks = autopickstation(wf2pick, pickparameter, fig_dict=fig_dict)
newpicks, _ = autopickstation(wf2pick, pickparameter, fig_dict=fig_dict)
# replace old dictionary with new one
picks[badpicks[i][0]] = newpicks

2
pylot/core/pick/charfuns.py
View File

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

3
pylot/core/pick/compare.py
View File

@ -506,7 +506,8 @@ class PDFstatistics(object):
return rlist
def writeThetaToFile(self, array, out_dir):
@staticmethod
def writeThetaToFile(array, out_dir):
"""
Method to write array like data to file. Useful since acquiring can take
serious amount of time when dealing with large databases.

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

@ -23,7 +23,7 @@ import warnings
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import argrelmax
from scipy.signal import argrelmax, argrelmin
from pylot.core.pick.charfuns import CharacteristicFunction
from pylot.core.pick.utils import getnoisewin, getsignalwin
@ -164,9 +164,9 @@ class AICPicker(AutoPicker):
iplot = int(self.iplot)
except:
if self.iplot == True or self.iplot == 'True':
iplot = 2
iplot = 2
else:
iplot = 0
iplot = 0
# find NaN's
nn = np.isnan(self.cf)
@ -191,16 +191,36 @@ class AICPicker(AutoPicker):
# remove offset in AIC function
offset = abs(min(aic) - min(aicsmooth))
aicsmooth = aicsmooth - offset
cf = self.Data[0].data
# get maximum of HOS/AR-CF as startimg point for searching
# minimum in AIC function
icfmax = np.argmax(self.Data[0].data)
icfmax = np.argmax(cf)
# MP MP testing threshold
thresh_hit = False
thresh_factor = 0.7
thresh = thresh_factor * cf[icfmax]
for index, sample in enumerate(cf):
if sample >= thresh:
thresh_hit = True
# go on searching for the following maximum
if index > 0 and thresh_hit:
if sample <= cf[index - 1]:
icfmax = index - 1
break
# MP MP ---
# find minimum in AIC-CF front of maximum of HOS/AR-CF
lpickwindow = int(round(self.PickWindow / self.dt))
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1):
if aicsmooth[i - 1] >= aicsmooth[i]:
self.Pick = self.Tcf[i]
break
tsafety = self.TSNR[1] # safety gap, AIC is usually a little bit too late
left_corner_ind = max([icfmax - lpickwindow, 2])
right_corner_ind = icfmax + int(tsafety / self.dt)
aic_snip = aicsmooth[left_corner_ind: right_corner_ind]
minima = argrelmin(aic_snip)[0] # 0th entry of tuples for axes
if len(minima) > 0:
pickindex = minima[-1] + left_corner_ind
self.Pick = self.Tcf[pickindex]
# if no minimum could be found:
# search in 1st derivative of AIC-CF
if self.Pick is None:
@ -215,17 +235,12 @@ class AICPicker(AutoPicker):
for i in range(icfmax - 1, max([icfmax - lpickwindow, 2]), -1):
if diffcf[i - 1] >= diffcf[i]:
self.Pick = self.Tcf[i]
pickindex = i
break
if self.Pick is not None:
# get noise window
inoise = getnoisewin(self.Tcf, self.Pick, self.TSNR[0], self.TSNR[1])
# check, if these are counts or m/s, important for slope estimation!
# this is quick and dirty, better solution?
if max(self.Data[0].data < 1e-3) and max(self.Data[0].data >= 1e-6):
self.Data[0].data = self.Data[0].data * 1000000.
elif max(self.Data[0].data < 1e-6):
self.Data[0].data = self.Data[0].data * 1e13
# get signal window
isignal = getsignalwin(self.Tcf, self.Pick, self.TSNR[2])
if len(isignal) == 0:
@ -233,33 +248,39 @@ class AICPicker(AutoPicker):
ii = min([isignal[len(isignal) - 1], len(self.Tcf)])
isignal = isignal[0:ii]
try:
self.Data[0].data[isignal]
cf[isignal]
except IndexError as e:
msg = "Time series out of bounds! {}".format(e)
print(msg)
return
# calculate SNR from CF
self.SNR = max(abs(self.Data[0].data[isignal] - np.mean(self.Data[0].data[isignal]))) / \
max(abs(self.Data[0].data[inoise] - np.mean(self.Data[0].data[inoise])))
self.SNR = max(abs(cf[isignal])) / \
abs(np.mean(cf[inoise]))
# calculate slope from CF after initial pick
# get slope window
tslope = self.TSNR[3] # slope determination window
islope = np.where((self.Tcf <= min([self.Pick + tslope, self.Tcf[-1]])) \
& (self.Tcf >= self.Pick)) # TODO: put this in a seperate function like getsignalwin
if tsafety >= 0:
islope = np.where((self.Tcf <= min([self.Pick + tslope + tsafety, self.Tcf[-1]])) \
& (self.Tcf >= self.Pick)) # TODO: put this in a seperate function like getsignalwin
else:
islope = np.where((self.Tcf <= min([self.Pick + tslope, self.Tcf[-1]])) \
& (
self.Tcf >= self.Pick + tsafety)) # TODO: put this in a seperate function like getsignalwin
# find maximum within slope determination window
# 'cause slope should be calculated up to first local minimum only!
try:
dataslope = self.Data[0].data[islope[0][0:-1]]
dataslope = cf[islope[0][0:-1]]
except IndexError:
print("Slope Calculation: empty array islope, check signal window")
return
if len(dataslope) <= 1:
if len(dataslope) < 2:
print('No or not enough data in slope window found!')
return
imaxs, = argrelmax(dataslope)
if imaxs.size:
try:
imaxs, = argrelmax(dataslope)
imax = imaxs[0]
else:
except (ValueError, IndexError) as e:
print(e, 'picker: argrelmax not working!')
imax = np.argmax(dataslope)
iislope = islope[0][0:imax + 1]
if len(iislope) < 2:
@ -271,52 +292,64 @@ class AICPicker(AutoPicker):
print("AICPicker: Maximum for slope determination right at the beginning of the window!")
print("Choose longer slope determination window!")
if self.iplot > 1:
if self.fig == None or self.fig == 'None':
if self.fig is None or self.fig == 'None':
fig = plt.figure()
plt_flag = 1
plt_flag = iplot
else:
fig = self.fig
ax = fig.add_subplot(111)
x = self.Data[0].data
ax.plot(self.Tcf, x / max(x), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
cf = cf
ax.plot(self.Tcf, cf / max(cf), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
ax.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
ax.legend(loc=1)
ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax.set_yticks([])
ax.set_title(self.Data[0].stats.station)
if plt_flag == 1:
if plt_flag in [1, 2]:
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return
iislope = islope[0][0:imax+1]
dataslope = self.Data[0].data[iislope]
iislope = islope[0][0:imax + 1]
# MP MP change slope calculation
# get all maxima of aicsmooth
iaicmaxima = argrelmax(aicsmooth)[0]
# get first index of maximum after pickindex (indices saved in iaicmaxima)
aicmax = iaicmaxima[np.where(iaicmaxima > pickindex)[0]]
if len(aicmax) > 0:
iaicmax = aicmax[0]
else:
iaicmax = -1
dataslope = aicsmooth[pickindex: iaicmax]
# calculate slope as polynomal fit of order 1
xslope = np.arange(0, len(dataslope), 1)
P = np.polyfit(xslope, dataslope, 1)
datafit = np.polyval(P, xslope)
if datafit[0] >= datafit[-1]:
print('AICPicker: Negative slope, bad onset skipped!')
return
self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
else:
self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
# normalize slope to maximum of cf to make it unit independent
self.slope /= aicsmooth[iaicmax]
else:
self.SNR = None
self.slope = None
if iplot > 1:
if self.fig == None or self.fig == 'None':
if self.fig is None or self.fig == 'None':
fig = plt.figure() # self.iplot)
plt_flag = 1
plt_flag = iplot
else:
fig = self.fig
fig._tight = True
ax1 = fig.add_subplot(211)
x = self.Data[0].data
if len(self.Tcf) > len(self.Data[0].data): # why? LK
self.Tcf = self.Tcf[0:len(self.Tcf)-1]
ax1.plot(self.Tcf, x / max(x), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
if len(self.Tcf) > len(cf): # why? LK
self.Tcf = self.Tcf[0:len(self.Tcf) - 1]
ax1.plot(self.Tcf, cf / max(cf), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
ax1.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
if self.Pick is not None:
ax1.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b', linewidth=2, label='AIC-Pick')
@ -326,7 +359,7 @@ class AICPicker(AutoPicker):
if self.Pick is not None:
ax2 = fig.add_subplot(2, 1, 2, sharex=ax1)
ax2.plot(self.Tcf, x, color=self._linecolor, linewidth=0.7, label='Data')
ax2.plot(self.Tcf, aicsmooth, color='r', linewidth=0.7, label='Data')
ax1.axvspan(self.Tcf[inoise[0]], self.Tcf[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
ax1.axvspan(self.Tcf[isignal[0]], self.Tcf[isignal[-1]], color='b', alpha=0.2, lw=0,
label='Signal Window')
@ -338,28 +371,34 @@ class AICPicker(AutoPicker):
label='Signal Window')
ax2.axvspan(self.Tcf[iislope[0]], self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0,
label='Slope Window')
ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2, label='Slope')
ax2.plot(self.Tcf[pickindex: iaicmax], datafit, 'g', linewidth=2,
label='Slope') # MP MP changed temporarily!
ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope))
if self.slope is not None:
ax1.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope))
else:
ax1.set_title('Station %s, SNR=%7.2f' % (self.Data[0].stats.station, self.SNR))
ax2.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax2.set_ylabel('Counts')
ax2.set_yticks([])
ax2.legend(loc=1)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
plt.close(fig)
else:
ax1.set_title(self.Data[0].stats.station)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
plt.close(fig)
if self.Pick == None:
if plt_flag in [1, 2]:
fig.show()
try:
input()
except SyntaxError:
pass
plt.close(fig)
if plt_flag == 3:
stats = self.Data[0].stats
netstlc = '{}.{}.{}'.format(stats.network, stats.station, stats.location)
fig.savefig('aicfig_{}_{}.png'.format(netstlc, stats.channel))
if self.Pick is None:
print('AICPicker: Could not find minimum, picking window too short?')
return
@ -376,9 +415,9 @@ class PragPicker(AutoPicker):
iplot = int(self.getiplot())
except:
if self.getiplot() == True or self.getiplot() == 'True':
iplot = 2
iplot = 2
else:
iplot = 0
iplot = 0
if self.getpick1() is not None:
print('PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...')
@ -417,11 +456,11 @@ class PragPicker(AutoPicker):
# prominent trend: decrease aus
# flat: use given aus
cfdiff = np.diff(cfipick)
if len(cfdiff)<20:
if len(cfdiff) < 20:
print('PragPicker: Very few samples for CF. Check LTA window dimensions!')
i0diff = np.where(cfdiff > 0)
cfdiff = cfdiff[i0diff]
if len(cfdiff)<1:
if len(cfdiff) < 1:
print('PragPicker: Negative slope for CF. Check LTA window dimensions! STOP')
self.Pick = None
return
@ -445,7 +484,7 @@ class PragPicker(AutoPicker):
break
# now we look to the left
if len(self.cf) > ipick1 +1:
if len(self.cf) > ipick1 + 1:
for i in range(ipick1, max([ipick1 - lpickwindow + 1, 2]), -1):
if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]:
if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
@ -456,7 +495,7 @@ class PragPicker(AutoPicker):
cfpick_l = self.cf[i]
break
else:
msg ='PragPicker: Initial onset too close to start of CF! \
msg = 'PragPicker: Initial onset too close to start of CF! \
Stop finalizing pick to the left.'
print(msg)
@ -476,7 +515,7 @@ class PragPicker(AutoPicker):
pickflag = 0
if iplot > 1:
if self.fig == None or self.fig == 'None':
if self.fig is None or self.fig == 'None':
fig = plt.figure() # self.getiplot())
plt_flag = 1
else:
@ -486,15 +525,18 @@ class PragPicker(AutoPicker):
ax.plot(Tcfpick, cfipick, color=self._linecolor, linewidth=0.7, label='CF')
ax.plot(Tcfpick, cfsmoothipick, 'r', label='Smoothed CF')
if pickflag > 0:
ax.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], self._pickcolor_p, linewidth=2, label='Pick')
ax.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], self._pickcolor_p, linewidth=2,
label='Pick')
ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax.set_yticks([])
ax.set_title(self.Data[0].stats.station)
ax.legend(loc=1)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return

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

@ -13,8 +13,7 @@ import warnings
import matplotlib.pyplot as plt
import numpy as np
from obspy.core import Stream, UTCDateTime
from pylot.core.util.utils import real_Bool, real_None
from pylot.core.util.utils import real_Bool, real_None, SetChannelComponents
def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecolor='k'):
@ -143,13 +142,16 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol
ax.plot(t, x, color=linecolor, linewidth=0.7, label='Data')
ax.axvspan(t[inoise[0]], t[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
ax.axvspan(t[isignal[0]], t[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window')
ax.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], color=linecolor, linewidth=0.7, linestyle='dashed', label='Noise Level')
ax.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], color=linecolor, linewidth=0.7, linestyle='dashed',
label='Noise Level')
ax.plot(t[pis[zc]], np.zeros(len(zc)), '*g',
markersize=14, label='Zero Crossings')
ax.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], color=linecolor, linewidth=0.7, linestyle='dashed')
ax.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2, label='mpp')
ax.plot([LPick, LPick], [max(x) / 2, -max(x) / 2], color=linecolor, linewidth=0.7, linestyle='dashed', label='lpp')
ax.plot([EPick, EPick], [max(x) / 2, -max(x) / 2], color=linecolor, linewidth=0.7, linestyle='dashed', label='epp')
ax.plot([LPick, LPick], [max(x) / 2, -max(x) / 2], color=linecolor, linewidth=0.7, linestyle='dashed',
label='lpp')
ax.plot([EPick, EPick], [max(x) / 2, -max(x) / 2], color=linecolor, linewidth=0.7, linestyle='dashed',
label='epp')
ax.plot([Pick1 + PickError, Pick1 + PickError],
[max(x) / 2, -max(x) / 2], 'r--', label='spe')
ax.plot([Pick1 - PickError, Pick1 - PickError],
@ -162,8 +164,10 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol
ax.legend(loc=1)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return EPick, LPick, PickError
@ -197,9 +201,9 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
iplot = int(iplot)
except:
if iplot == True or iplot == 'True':
iplot = 2
iplot = 2
else:
iplot = 0
iplot = 0
warnings.simplefilter('ignore', np.RankWarning)
@ -226,8 +230,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
# get zero crossings after most likely pick
# initial onset is assumed to be the first zero crossing
# first from unfiltered trace
zc1 = []
zc1.append(Pick)
zc1 = [Pick]
index1 = []
i = 0
for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]):
@ -272,8 +275,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
# now using filterd trace
# next zero crossings after most likely pick
zc2 = []
zc2.append(Pick)
zc2 = [Pick]
index2 = []
i = 0
for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]):
@ -361,8 +363,10 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
ax2.set_yticks([])
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return FM
@ -552,8 +556,6 @@ def select_for_phase(st, phase):
:rtype: `~obspy.core.stream.Stream`
"""
from pylot.core.util.defaults import SetChannelComponents
sel_st = Stream()
compclass = SetChannelComponents()
if phase.upper() == 'P':
@ -602,14 +604,18 @@ def wadaticheck(pickdic, dttolerance, iplot=0, fig_dict=None):
ibad = 0
for key in list(pickdic.keys()):
if pickdic[key]['P']['weight'] < 4 and pickdic[key]['S']['weight'] < 4:
ppick = pickdic[key].get('P')
spick = pickdic[key].get('S')
if not ppick or not spick:
continue
if ppick['weight'] < 4 and spick['weight'] < 4:
# calculate S-P time
spt = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp']
spt = spick['mpp'] - ppick['mpp']
# add S-P time to dictionary
pickdic[key]['SPt'] = spt
# add P onsets and corresponding S-P times to list
UTCPpick = UTCDateTime(pickdic[key]['P']['mpp'])
UTCSpick = UTCDateTime(pickdic[key]['S']['mpp'])
UTCPpick = UTCDateTime(ppick['mpp'])
UTCSpick = UTCDateTime(spick['mpp'])
Ppicks.append(UTCPpick.timestamp)
Spicks.append(UTCSpick.timestamp)
SPtimes.append(spt)
@ -639,11 +645,11 @@ def wadaticheck(pickdic, dttolerance, iplot=0, fig_dict=None):
# check, if deviation is larger than adjusted
if wddiff > dttolerance:
# remove pick from dictionary
pickdic.pop(key)
# # mark onset and downgrade S-weight to 9
# # mark onset and downgrade S-weight to 9, also set SPE to None (disregarded in GUI)
# # (not used anymore)
# marker = 'badWadatiCheck'
# pickdic[key]['S']['weight'] = 9
marker = 'badWadatiCheck'
pickdic[key]['S']['weight'] = 9
pickdic[key]['S']['spe'] = None
badstations.append(key)
ibad += 1
else:
@ -655,8 +661,7 @@ def wadaticheck(pickdic, dttolerance, iplot=0, fig_dict=None):
checkedSPtime = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp']
checkedSPtimes.append(checkedSPtime)
pickdic[key]['S']['marked'] = marker
#pickdic[key]['S']['marked'] = marker
pickdic[key]['S']['marked'] = marker
print("wadaticheck: the following stations failed the check:")
print(badstations)
@ -684,7 +689,7 @@ def wadaticheck(pickdic, dttolerance, iplot=0, fig_dict=None):
wfitflag = 1
# plot results
if iplot > 0:
if iplot > 0 or fig_dict:
if fig_dict:
fig = fig_dict['wadati']
linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl']
@ -698,8 +703,8 @@ def wadaticheck(pickdic, dttolerance, iplot=0, fig_dict=None):
ax.plot(Ppicks, SPtimes, 'ro', label='Skipped S-Picks')
if wfitflag == 0:
ax.plot(Ppicks, wdfit, color=linecolor, linewidth=0.7, label='Wadati 1')
ax.plot(Ppicks, wdfit+dttolerance, color='0.9', linewidth=0.5, label='Wadati 1 Tolerance')
ax.plot(Ppicks, wdfit-dttolerance, color='0.9', linewidth=0.5)
ax.plot(Ppicks, wdfit + dttolerance, color='0.9', linewidth=0.5, label='Wadati 1 Tolerance')
ax.plot(Ppicks, wdfit - dttolerance, color='0.9', linewidth=0.5)
ax.plot(checkedPpicks, wdfit2, 'g', label='Wadati 2')
ax.plot(checkedPpicks, checkedSPtimes, color=linecolor,
linewidth=0, marker='o', label='Reliable S-Picks')
@ -766,9 +771,9 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot=0, fi
iplot = int(iplot)
except:
if real_Bool(iplot):
iplot = 2
iplot = 2
else:
iplot = 0
iplot = 0
assert isinstance(X, Stream), "%s is not a stream object" % str(X)
@ -830,8 +835,10 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot=0, fi
ax.set_yticks([])
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return returnflag
@ -865,9 +872,12 @@ def checkPonsets(pickdic, dttolerance, jackfactor=5, iplot=0, fig_dict=None):
Ppicks = []
stations = []
for station in pickdic:
if pickdic[station]['P']['weight'] < 4:
pick = pickdic[station].get('P')
if not pick:
continue
if pick['weight'] < 4:
# add P onsets to list
UTCPpick = UTCDateTime(pickdic[station]['P']['mpp'])
UTCPpick = UTCDateTime(pick['mpp'])
Ppicks.append(UTCPpick.timestamp)
stations.append(station)
@ -926,7 +936,7 @@ def checkPonsets(pickdic, dttolerance, jackfactor=5, iplot=0, fig_dict=None):
checkedonsets = pickdic
if iplot > 0:
if iplot > 0 or fig_dict:
if fig_dict:
fig = fig_dict['jackknife']
plt_flag = 0
@ -1056,16 +1066,15 @@ def checkZ4S(X, pick, zfac, checkwin, iplot, fig=None, linecolor='k'):
:return: returnflag; 0 if onset failed test, 1 if onset passed test
:rtype: int
"""
plt_flag = 0
try:
iplot = int(iplot)
except:
if real_Bool(iplot):
iplot = 2
iplot = 2
else:
iplot = 0
iplot = 0
assert isinstance(X, Stream), "%s is not a stream object" % str(X)
@ -1167,12 +1176,81 @@ def checkZ4S(X, pick, zfac, checkwin, iplot, fig=None, linecolor='k'):
ax.set_xlabel('Time [s] since %s' % zdat[0].stats.starttime)
if plt_flag == 1:
fig.show()
try: input()
except SyntaxError: pass
try:
input()
except SyntaxError:
pass
plt.close(fig)
return returnflag
def getPickQuality(wfdata, picks, inputs, phase, compclass=None):
quality = 4
components4phases = {'P': ['Z'],
'S': ['N', 'E']}
timeErrors4phases = {'P': 'timeerrorsP',
'S': 'timeerrorsS'}
tsnr4phases = {'P': 'tsnrz',
'S': 'tsnrh'}
if not phase in components4phases.keys():
raise IOError('getPickQuality: Could not understand phase: {}'.format(phase))
if not compclass:
print('Warning: No settings for channel components found. Using default')
compclass = SetChannelComponents()
picks = picks[phase]
mpp = picks.get('mpp')
uncertainty = picks.get('spe')
if not mpp:
print('getPickQuality: No pick found!')
return quality
if not uncertainty:
print('getPickQuality: No pick uncertainty (spe) found!')
return quality
tsnr = inputs[tsnr4phases[phase]]
timeErrors = inputs[timeErrors4phases[phase]]
snrdb_final = 0
for component in components4phases[phase]:
alter_comp = compclass.getCompPosition(component)
st_select = wfdata.select(component=component)
st_select += wfdata.select(component=alter_comp)
if st_select:
trace = st_select[0]
_, snrdb, _ = getSNR(st_select, tsnr,
mpp - trace.stats.starttime)
if snrdb > snrdb_final:
snrdb_final = snrdb
quality = getQualityFromUncertainty(uncertainty, timeErrors)
quality += getQualityFromSNR(snrdb_final)
return quality
def getQualityFromSNR(snrdb):
quality_modifier = 4
if not snrdb:
print('getQualityFromSNR: No snrdb!')
return quality_modifier
# MP MP ++++ experimental,
# raise pick quality by x classes if snrdb is lower than corresponding key
quality4snrdb = {3: 4,
5: 3,
7: 2,
9: 1,
11: 0}
# MP MP ---
# iterate over all thresholds and check whether snrdb is larger, if so, set new quality_modifier
for snrdb_threshold in sorted(list(quality4snrdb.keys())):
if snrdb > snrdb_threshold:
quality_modifier = quality4snrdb[snrdb_threshold]
return quality_modifier
def getQualityFromUncertainty(uncertainty, Errors):
"""
Script to transform uncertainty into quality classes 0-4 regarding adjusted time errors
@ -1193,19 +1271,20 @@ def getQualityFromUncertainty(uncertainty, Errors):
if uncertainty <= Errors[0]:
quality = 0
elif (uncertainty > Errors[0]) and \
(uncertainty < Errors[1]):
(uncertainty <= Errors[1]):
quality = 1
elif (uncertainty > Errors[1]) and \
(uncertainty < Errors[2]):
(uncertainty <= Errors[2]):
quality = 2
elif (uncertainty > Errors[2]) and \
(uncertainty < Errors[3]):
(uncertainty <= Errors[3]):
quality = 3
elif uncertainty > Errors[3]:
quality = 4
return quality
if __name__ == '__main__':
import doctest

448
pylot/core/util/array_map.py Normal file
View File

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

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

@ -12,6 +12,303 @@ from pylot.core.util.utils import key_for_set_value, find_in_list, \
remove_underscores, gen_Pool
class Metadata(object):
def __init__(self, inventory=None, verbosity=1):
self.inventories = []
# saves read metadata objects (Parser/inventory) for a filename
self.inventory_files = {}
# saves filenames holding metadata for a seed_id
# seed id as key, path to file as value
self.seed_ids = {}
self.stations_dict = {}
if inventory:
if os.path.isdir(inventory):
self.add_inventory(inventory)
if os.path.isfile(inventory):
self.add_inventory_file(inventory)
self.verbosity = verbosity
def __str__(self):
repr = 'PyLoT Metadata object including the following inventories:\n\n'
ntotal = len(self.inventories)
for index, inventory in enumerate(self.inventories):
if index < 2 or (ntotal - index) < 3:
repr += '{}\n'.format(inventory)
if ntotal > 4 and int(ntotal / 2) == index:
repr += '...\n'
if ntotal > 4:
repr += '\nTotal of {} inventories. Use Metadata.inventories to see all.'.format(ntotal)
return repr
def __repr__(self):
return self.__str__()
def add_inventory(self, path_to_inventory):
"""
Add path to list of inventories.
:param path_to_inventory: Path to a folder
:type path_to_inventory: str
:return: None
"""
assert (os.path.isdir(path_to_inventory)), '{} is no directory'.format(path_to_inventory)
if not path_to_inventory in self.inventories:
self.inventories.append(path_to_inventory)
def add_inventory_file(self, path_to_inventory_file):
"""
Add the folder in which the file exists to the list of inventories.
:param path_to_inventory_file: full path including filename
:type path_to_inventory_file: str
:return: None
"""
assert (os.path.isfile(path_to_inventory_file)), '{} is no file'.format(path_to_inventory_file)
self.add_inventory(os.path.split(path_to_inventory_file)[0])
if not path_to_inventory_file in self.inventory_files.keys():
self.read_single_file(path_to_inventory_file)
def remove_all_inventories(self):
self.__init__()
def remove_inventory(self, path_to_inventory):
"""
Remove a path from inventories list. If path is not in inventories list, do nothing.
:param path_to_inventory: Path to a folder
"""
if not path_to_inventory in self.inventories:
print('Path {} not in inventories list.'.format(path_to_inventory))
return
self.inventories.remove(path_to_inventory)
for filename in self.inventory_files.keys():
if filename.startswith(path_to_inventory):
del (self.inventory_files[filename])
for seed_id in self.seed_ids.keys():
if self.seed_ids[seed_id].startswith(path_to_inventory):
del (self.seed_ids[seed_id])
def get_metadata(self, seed_id, time=None):
"""
Get metadata for seed id at time. When time is not specified, metadata for current time is fetched.
:param seed_id: Seed id such as BW.WETR..HHZ (Network.Station.Location.Channel)
:type seed_id: str
:param time: Time for which the metadata should be returned
:type time: UTCDateTime
:return: Dictionary with keys data and invtype.
data is a obspy.io.xseed.parser.Parser or an obspy.core.inventory.inventory.Inventory depending on the metadata
file.
invtype is a string denoting of which type the value of the data key is. It can take the values 'dless',
'dseed', 'xml', 'resp', according to the filetype of the metadata.
:rtype: dict
"""
# try most recent data if no time is specified
if not time:
time = UTCDateTime()
# get metadata for a specific seed_id, if not already read, try to read from inventories
if not seed_id in self.seed_ids.keys():
self._read_inventory_data(seed_id)
# if seed id is not found read all inventories and try to find it there
if not seed_id in self.seed_ids.keys():
if self.verbosity:
print('No data found for seed id {}. Trying to find it in all known inventories...'.format(seed_id))
self.read_all()
for inv_fname, metadata_dict in self.inventory_files.items():
# use get_coordinates to check for seed_id
try:
metadata_dict['data'].get_coordinates(seed_id, time)
self.seed_ids[seed_id] = inv_fname
if self.verbosity:
print('Found metadata for station {}!'.format(seed_id))
return metadata_dict
except Exception as e:
continue
print('Could not find metadata for station {}'.format(seed_id))
return None
fname = self.seed_ids[seed_id]
return self.inventory_files[fname]
def read_all(self):
"""
Read all metadata files found in all inventories
"""
for inventory in self.inventories:
for inv_fname in os.listdir(inventory):
inv_fname = os.path.join(inventory, inv_fname)
if not self.read_single_file(inv_fname):
continue
def read_single_file(self, inv_fname):
"""
Try to read a single file as Parser/Inventory and add its dictionary to inventory files if reading sudceeded.
:param inv_fname: path/filename of inventory file
:type inv_fname: str
:rtype: None
"""
# return if it was read already
if self.inventory_files.get(inv_fname, None):
return
try:
invtype, robj = self._read_metadata_file(inv_fname)
if robj is None:
return
except Exception as e:
print('Could not read file {}'.format(inv_fname))
return
self.inventory_files[inv_fname] = {'invtype': invtype,
'data': robj}
return True
def get_coordinates(self, seed_id, time=None):
"""
Get coordinates of given seed id.
:param seed_id: Seed id such as BW.WETR..HHZ (Network.Station.Location.Channel)
:type seed_id: str
:param time: Used when a station has data available at multiple time intervals
:type time: UTCDateTime
:return: dict containing position information of the station
:rtype: dict
"""
# try most recent data if no time is specified
if not time:
time = UTCDateTime()
metadata = self.get_metadata(seed_id, time)
if not metadata:
return
return metadata['data'].get_coordinates(seed_id, time)
def get_all_coordinates(self):
def stat_info_from_parser(parser):
for station in parser.stations:
station_name = station[0].station_call_letters
network_name = station[0].network_code
if not station_name in self.stations_dict.keys():
st_id = '{}.{}'.format(network_name, station_name)
self.stations_dict[st_id] = {'latitude': station[0].latitude,
'longitude': station[0].longitude}
def stat_info_from_inventory(inventory):
for network in inventory.networks:
for station in network.stations:
station_name = station.code
network_name = network_name.code
if not station_name in self.stations_dict.keys():
st_id = '{}.{}'.format(network_name, station_name)
self.stations_dict[st_id] = {'latitude': station[0].latitude,
'longitude': station[0].longitude}
read_stat = {'xml': stat_info_from_inventory,
'dless': stat_info_from_parser}
self.read_all()
for item in self.inventory_files.values():
inventory = item['data']
invtype = item['invtype']
read_stat[invtype](inventory)
return self.stations_dict
def get_paz(self, seed_id, time):
"""
:param seed_id: Seed id such as BW.WETR..HHZ (Network.Station.Location.Channel)
:type seed_id: str
:param time: Used when a station has data available at multiple time intervals
:type time: UTCDateTime
:rtype: dict
"""
metadata = self.get_metadata(seed_id)
if not metadata:
return
if metadata['invtype'] in ['dless', 'dseed']:
return metadata['data'].get_paz(seed_id, time)
elif metadata['invtype'] in ['resp', 'xml']:
resp = metadata['data'].get_response(seed_id, time)
return resp.get_paz(seed_id)
def _read_inventory_data(self, seed_id=None):
for inventory in self.inventories:
if self._read_metadata_iterator(path_to_inventory=inventory, station_seed_id=seed_id):
return
def _read_metadata_iterator(self, path_to_inventory, station_seed_id):
"""
Search for metadata for a specific station iteratively.
"""
station, network, location, channel = station_seed_id.split('.')
# seach for station seed id in filenames in invetory
fnames = glob.glob(os.path.join(path_to_inventory, '*' + station_seed_id + '*'))
if not fnames:
# search for station name in filename
fnames = glob.glob(os.path.join(path_to_inventory, '*' + station + '*'))
if not fnames:
# search for network name in filename
fnames = glob.glob(os.path.join(path_to_inventory, '*' + network + '*'))
if not fnames:
if self.verbosity:
print('Could not find filenames matching station name, network name or seed id')
return
for fname in fnames:
if fname in self.inventory_files.keys():
if self.inventory_files[fname]:
# file already read
continue
invtype, robj = self._read_metadata_file(os.path.join(path_to_inventory, fname))
try:
robj.get_coordinates(station_seed_id)
self.inventory_files[fname] = {'invtype': invtype,
'data': robj}
if station_seed_id in self.seed_ids.keys():
print('WARNING: Overwriting metadata for station {}'.format(station_seed_id))
self.seed_ids[station_seed_id] = fname
return True
except Exception as e:
continue
print('Could not find metadata for station_seed_id {} in path {}'.format(station_seed_id, path_to_inventory))
def _read_metadata_file(self, path_to_inventory_filename):
"""
function reading metadata files (either dataless seed, xml or resp)
:param path_to_inventory_filename:
:return: file type/ending, inventory object (Parser or Inventory)
:rtype: (str, obspy.io.xseed.Parser or obspy.core.inventory.inventory.Inventory)
"""
# functions used to read metadata for different file endings (or file types)
read_functions = {'dless': self._read_dless,
'dseed': self._read_dless,
'xml': self._read_inventory_file,
'resp': self._read_inventory_file}
file_ending = path_to_inventory_filename.split('.')[-1]
if file_ending in read_functions.keys():
robj, exc = read_functions[file_ending](path_to_inventory_filename)
if exc is not None:
raise exc
return file_ending, robj
# in case file endings did not match the above keys, try and error
for file_type in ['dless', 'xml']:
robj, exc = read_functions[file_type](path_to_inventory_filename)
if exc is None:
return file_type, robj
return None, None
@staticmethod
def _read_dless(path_to_inventory):
exc = None
try:
parser = Parser(path_to_inventory)
except Exception as exc:
parser = None
return parser, exc
@staticmethod
def _read_inventory_file(path_to_inventory):
exc = None
try:
inv = read_inventory(path_to_inventory)
except Exception as exc:
inv = None
return inv, exc
def time_from_header(header):
"""
Function takes in the second line from a .gse file and takes out the date and time from that line.
@ -196,12 +493,51 @@ def read_metadata(path_to_inventory):
return invtype, robj
# idea to optimize read_metadata
# def read_metadata_new(path_to_inventory):
# metadata_objects = []
# # read multiple files from directory
# if os.path.isdir(path_to_inventory):
# fnames = os.listdir(path_to_inventory)
# # read single file
# elif os.path.isfile(path_to_inventory):
# fnames = [path_to_inventory]
# else:
# print("Neither dataless-SEED file, inventory-xml file nor "
# "RESP-file found!")
# print("!!WRONG CALCULATION OF SOURCE PARAMETERS!!")
# fnames = []
#
# for fname in fnames:
# path_to_inventory_filename = os.path.join(path_to_inventory, fname)
# try:
# ftype, robj = read_metadata_file(path_to_inventory_filename)
# metadata_objects.append((ftype, robj))
# except Exception as e:
# print('Could not read metadata file {} '
# 'because of the following Exception: {}'.format(path_to_inventory_filename, e))
# return metadata_objects
def restitute_trace(input_tuple):
tr, invtype, inobj, unit, force = input_tuple
def no_metadata(tr, seed_id):
print('no metadata file found '
'for trace {0}'.format(seed_id))
return tr, True
tr, metadata, unit, force = input_tuple
remove_trace = False
seed_id = tr.get_id()
mdata = metadata.get_metadata(seed_id, time=tr.stats.starttime)
if not mdata:
return no_metadata(tr, seed_id)
invtype = mdata['invtype']
inobj = mdata['data']
# check, whether this trace has already been corrected
if 'processing' in tr.stats.keys() \
and np.any(['remove' in p for p in tr.stats.processing]) \
@ -213,8 +549,7 @@ def restitute_trace(input_tuple):
if invtype == 'resp':
fresp = find_in_list(inobj, seed_id)
if not fresp:
raise IOError('no response file found '
'for trace {0}'.format(seed_id))
return no_metadata(tr, seed_id)
fname = fresp
seedresp = dict(filename=fname,
date=stime,
@ -236,9 +571,8 @@ def restitute_trace(input_tuple):
else:
finv = invlist[0]
inventory = read_inventory(finv, format='STATIONXML')
elif invtype == None:
print("No restitution possible, as there are no station-meta data available!")
return tr, True
elif invtype is None:
return no_metadata(tr, seed_id)
else:
remove_trace = True
# apply restitution to data
@ -269,14 +603,11 @@ def restitute_trace(input_tuple):
return tr, remove_trace
def restitute_data(data, invtype, inobj, unit='VEL', force=False, ncores=0):
def restitute_data(data, metadata, unit='VEL', force=False, ncores=0):
"""
takes a data stream and a path_to_inventory and returns the corrected
waveform data stream
:param data: seismic data stream
:param invtype: type of found metadata
:param inobj: either list of metadata files or `obspy.io.xseed.Parser`
object
:param unit: unit to correct for (default: 'VEL')
:param force: force restitution for already corrected traces (default:
False)
@ -285,16 +616,16 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False, ncores=0):
restflag = list()
data = remove_underscores(data)
#data = remove_underscores(data)
# loop over traces
input_tuples = []
for tr in data:
input_tuples.append((tr, invtype, inobj, unit, force))
input_tuples.append((tr, metadata, unit, force))
data.remove(tr)
pool = gen_Pool(ncores)
result = pool.map(restitute_trace, input_tuples)
result = pool.imap_unordered(restitute_trace, input_tuples)
pool.close()
for tr, remove_trace in result:
@ -344,7 +675,7 @@ def get_prefilt(trace, tlow=(0.5, 0.9), thi=(5., 2.), verbosity=0):
fny = trace.stats.sampling_rate / 2
fc21 = fny - (fny * thi[0] / 100.)
fc22 = fny - (fny * thi[1] / 100.)
return (tlow[0], tlow[1], fc21, fc22)
return tlow[0], tlow[1], fc21, fc22
if __name__ == "__main__":

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

@ -16,7 +16,6 @@ from pylot.core.loc import hyposat
from pylot.core.loc import nll
from pylot.core.loc import velest
# determine system dependent path separator
system_name = platform.system()
if system_name in ["Linux", "Darwin"]:
@ -42,54 +41,3 @@ OUTPUTFORMATS = {'.xml': 'QUAKEML',
LOCTOOLS = dict(nll=nll, hyposat=hyposat, velest=velest, hypo71=hypo71, hypodd=hypodd)
class SetChannelComponents(object):
def __init__(self):
self.setDefaultCompPosition()
def setDefaultCompPosition(self):
# default component order
self.compPosition_Map = dict(Z=2, N=1, E=0)
self.compName_Map = {'3': 'Z',
'1': 'N',
'2': 'E'}
def _getCurrentPosition(self, component):
for key, value in self.compName_Map.items():
if value == component:
return key, value
errMsg = 'getCurrentPosition: Could not find former position of component {}.'.format(component)
raise ValueError(errMsg)
def _switch(self, component, component_alter):
# Without switching, multiple definitions of the same alter_comp are possible
old_alter_comp, _ = self._getCurrentPosition(component)
old_comp = self.compName_Map[component_alter]
if not old_alter_comp == component_alter and not old_comp == component:
self.compName_Map[old_alter_comp] = old_comp
print('switch: Automatically switched component {} to {}'.format(old_alter_comp, old_comp))
def setCompPosition(self, component_alter, component, switch=True):
component_alter = str(component_alter)
if not component_alter in self.compName_Map.keys():
errMsg = 'setCompPosition: Unrecognized alternative component {}. Expecting one of {}.'
raise ValueError(errMsg.format(component_alter, self.compName_Map.keys()))
if not component in self.compPosition_Map.keys():
errMsg = 'setCompPosition: Unrecognized target component {}. Expecting one of {}.'
raise ValueError(errMsg.format(component, self.compPosition_Map.keys()))
print('setCompPosition: set component {} to {}'.format(component_alter, component))
if switch:
self._switch(component, component_alter)
self.compName_Map[component_alter] = component
def getCompPosition(self, component):
return self._getCurrentPosition(component)[0]
def getPlotPosition(self, component):
component = str(component)
if component in self.compPosition_Map.keys():
return self.compPosition_Map[component]
elif component in self.compName_Map.keys():
return self.compPosition_Map[self.compName_Map[component]]
else:
errMsg = 'getCompPosition: Unrecognized component {}. Expecting one of {} or {}.'
raise ValueError(errMsg.format(component, self.compPosition_Map.keys(), self.compName_Map.keys()))

31
pylot/core/util/event.py
View File

@ -7,6 +7,7 @@ from obspy import UTCDateTime
from obspy.core.event import Event as ObsPyEvent
from obspy.core.event import Origin, ResourceIdentifier
from pylot.core.io.phases import picks_from_picksdict
from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
class Event(ObsPyEvent):
@ -33,6 +34,8 @@ class Event(ObsPyEvent):
self._testEvent = False
self._refEvent = False
self.get_notes()
self.get_obspy_event_info()
self.dirty = False
def get_notes_path(self):
"""
@ -43,6 +46,18 @@ class Event(ObsPyEvent):
notesfile = os.path.join(self.path, 'notes.txt')
return notesfile
def get_obspy_event_info(self):
infile_pickle = os.path.join(self.path, 'info/event.pkl')
if not os.path.isfile(infile_pickle):
return
try:
event_dmt = qml_from_obspyDMT(infile_pickle)
except Exception as e:
print('Could not get obspy event info: {}'.format(e))
return
self.magnitudes = event_dmt.magnitudes
self.origins = event_dmt.origins
def get_notes(self):
"""
set self.note attribute to content of notes file
@ -129,6 +144,7 @@ class Event(ObsPyEvent):
for index, pick in reversed(list(enumerate(self.picks))):
if picktype in str(pick.method_id):
self.picks.pop(index)
self.dirty = True
def addPicks(self, picks):
"""
@ -143,12 +159,12 @@ class Event(ObsPyEvent):
# add ObsPy picks (clear old manual and copy all new manual from pylot)
self.clearObsPyPicks('manual')
self.picks += picks_from_picksdict(self.pylot_picks)
self.dirty = True
def addAutopicks(self, autopicks):
"""
Add automatic picks to event
:param autopicks: automatic picks to add to event
:type autopicks dict:
:return:
:rtype: None
"""
@ -157,6 +173,7 @@ class Event(ObsPyEvent):
# add ObsPy picks (clear old auto and copy all new auto from pylot)
self.clearObsPyPicks('auto')
self.picks += picks_from_picksdict(self.pylot_autopicks)
self.dirty = True
def setPick(self, station, pick):
"""
@ -172,11 +189,13 @@ class Event(ObsPyEvent):
self.pylot_picks[station] = pick
else:
try:
self.pylot_picks.pop(station)
if station in self.pylot_picks:
self.pylot_picks.pop(station)
except Exception as e:
print('Could not remove pick {} from station {}: {}'.format(pick, station, e))
self.clearObsPyPicks('manual')
self.picks += picks_from_picksdict(self.pylot_picks)
self.dirty = True
def setPicks(self, picks):
"""
@ -189,6 +208,7 @@ class Event(ObsPyEvent):
self.pylot_picks = picks
self.clearObsPyPicks('manual')
self.picks += picks_from_picksdict(self.pylot_picks)
self.dirty = True
def getPick(self, station):
"""
@ -223,11 +243,13 @@ class Event(ObsPyEvent):
self.pylot_autopicks[station] = pick
else:
try:
self.pylot_autopicks.pop(station)
if station in self.pylot_autopicks:
self.pylot_autopicks.pop(station)
except Exception as e:
print('Could not remove pick {} from station {}: {}'.format(pick, station, e))
self.clearObsPyPicks('auto')
self.picks += picks_from_picksdict(self.pylot_autopicks)
self.dirty = True
def setAutopicks(self, picks):
"""
@ -240,6 +262,7 @@ class Event(ObsPyEvent):
self.pylot_autopicks = picks
self.clearObsPyPicks('auto')
self.picks += picks_from_picksdict(self.pylot_autopicks)
self.dirty = True
def getAutopick(self, station):
"""
@ -278,6 +301,7 @@ class Event(ObsPyEvent):
try:
outfile = open(filename, 'wb')
cPickle.dump(self, outfile, -1)
self.dirty = False
except Exception as e:
print('Could not pickle PyLoT event. Reason: {}'.format(e))
@ -296,5 +320,6 @@ class Event(ObsPyEvent):
import _pickle as cPickle
infile = open(filename, 'rb')
event = cPickle.load(infile)
event.dirty = False
print('Loaded %s' % filename)
return event

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

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

20
pylot/core/util/obspyDMT_interface.py
View File

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

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

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

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

@ -33,31 +33,21 @@ class Thread(QThread):
self._executed = False
self._executedError = e
traceback.print_exc()
exctype, value = sys.exc_info ()[:2]
self._executedErrorInfo = '{} {} {}'.\
exctype, value = sys.exc_info()[:2]
self._executedErrorInfo = '{} {} {}'. \
format(exctype, value, traceback.format_exc())
sys.stdout = sys.__stdout__
def showProgressbar(self):
if 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)
# # generate widget if not given in init
# if not self.pb_widget:
# self.pb_widget = ProgressBarWidget(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.label.setText(self.progressText)
self.pb_widget.show()
def hideProgressbar(self):
@ -75,6 +65,7 @@ class Worker(QRunnable):
'''
Worker class to be run by MultiThread(QThread).
'''
def __init__(self, fun, args,
progressText=None,
pb_widget=None,
@ -82,7 +73,7 @@ class Worker(QRunnable):
super(Worker, self).__init__()
self.fun = fun
self.args = args
#self.kwargs = kwargs
# self.kwargs = kwargs
self.signals = WorkerSignals()
self.progressText = progressText
self.pb_widget = pb_widget
@ -96,9 +87,9 @@ class Worker(QRunnable):
try:
result = self.fun(self.args)
except:
exctype, value = sys.exc_info ()[:2]
exctype, value = sys.exc_info()[:2]
print(exctype, value, traceback.format_exc())
self.signals.error.emit ((exctype, value, traceback.format_exc ()))
self.signals.error.emit((exctype, value, traceback.format_exc()))
else:
self.signals.result.emit(result)
finally:
@ -140,13 +131,13 @@ class MultiThread(QThread):
def run(self):
if self.redirect_stdout:
sys.stdout = self
sys.stdout = self
try:
if not self.ncores:
self.ncores = multiprocessing.cpu_count()
pool = multiprocessing.Pool(self.ncores)
self.data = pool.map_async(self.func, self.args, callback=self.emitDone)
#self.data = pool.apply_async(self.func, self.shotlist, callback=self.emitDone) #emit each time returned
# self.data = pool.apply_async(self.func, self.shotlist, callback=self.emitDone) #emit each time returned
pool.close()
self._executed = True
except Exception as e:

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

@ -22,11 +22,7 @@ from pylot.styles import style_settings
from scipy.interpolate import splrep, splev
from PySide import QtCore, QtGui
try:
import pyqtgraph as pg
except Exception as e:
print('PyLoT: Could not import pyqtgraph. {}'.format(e))
pg = None
import pyqtgraph as pg
def _pickle_method(m):
if m.im_self is None:
@ -34,6 +30,7 @@ def _pickle_method(m):
else:
return getattr, (m.im_self, m.im_func.func_name)
def getAutoFilteroptions(phase, parameter):
filtername = {'P': 'bpz2',
'S': 'bph2'}
@ -41,9 +38,10 @@ def getAutoFilteroptions(phase, parameter):
print('autoPickParameter: No filter options for phase {}.'.format(phase))
return
freqmin, freqmax = parameter.get(filtername[phase])
filteroptions = FilterOptions(type='bandpass', freq=[freqmin, freqmax], order=4) # order=4 default from obspy
filteroptions = FilterOptions(type='bandpass', freq=[freqmin, freqmax], order=4) # order=4 default from obspy
return filteroptions
def readDefaultFilterInformation(fname):
"""
Read default filter information from pylot.in file
@ -118,8 +116,12 @@ def gen_Pool(ncores=0):
"""
import multiprocessing
if ncores == 0:
ncores = multiprocessing.cpu_count()
ncores_max = multiprocessing.cpu_count()
if ncores == 0 or ncores > ncores_max:
ncores = ncores_max
if ncores > ncores_max:
print('Reduced number of requested CPU slots to available number: {}'.format(ncores))
print('gen_Pool: Generated multiprocessing Pool with {} cores\n'.format(ncores))
@ -397,6 +399,10 @@ def full_range(stream):
:return: minimum start time and maximum end time
:rtype: (`~maximum start time and minimum end time`, maximum start time and minimum end time)
"""
if not stream:
print('full_range: Empty Stream!')
return None, None
min_start = min([trace.stats.starttime for trace in stream])
max_end = max([trace.stats.endtime for trace in stream])
@ -454,7 +460,8 @@ def getLogin():
:return: login ID
:rtype: str
"""
return os.getlogin()
import getpass
return getpass.getuser()
def getOwner(fn):
@ -536,7 +543,7 @@ def isSorted(iterable):
False
"""
assert isIterable(iterable), 'object is not iterable; object: {' \
'0}'.format(iterable)
'}'.format(iterable)
if type(iterable) is str:
iterable = [s for s in iterable]
return sorted(iterable) == iterable
@ -674,7 +681,7 @@ def pick_color(picktype, phase, quality=0):
bpc = base_phase_colors(picktype, phase) # returns dict like {'modifier': 'g', 'rgba': (0, 0, 255, 255)}
rgba = bpc['rgba']
modifier = bpc['modifier']
intensity = 255.*quality/min_quality
intensity = 255. * quality / min_quality
rgba = modify_rgba(rgba, modifier, intensity)
return rgba
@ -786,6 +793,7 @@ def base_phase_colors(picktype, phase):
phasecolors = style_settings.phasecolors
return phasecolors[picktype][phase]
def transform_colors_mpl_str(colors, no_alpha=False):
"""
Transforms rgba color values to a matplotlib string of color values with a range of [0, 1]
@ -804,6 +812,7 @@ def transform_colors_mpl_str(colors, no_alpha=False):
colors_mpl = '({}, {}, {}, {})'.format(*colors_mpl)
return colors_mpl
def transform_colors_mpl(colors):
"""
Transform rgba colors from [0, 255] to [0, 1]
@ -816,6 +825,7 @@ def transform_colors_mpl(colors):
colors_mpl = tuple([color / 255. for color in colors])
return colors_mpl
def remove_underscores(data):
"""
takes a `obspy.core.stream.Stream` object and removes all underscores
@ -825,9 +835,9 @@ def remove_underscores(data):
:return: data stream
:rtype: `~obspy.core.stream.Stream`
"""
for tr in data:
# remove underscores
tr.stats.station = tr.stats.station.strip('_')
#for tr in data:
# # remove underscores
# tr.stats.station = tr.stats.station.strip('_')
return data
@ -859,6 +869,19 @@ def trim_station_components(data, trim_start=True, trim_end=True):
return data
def merge_stream(stream):
gaps = stream.get_gaps()
if gaps:
# list of merged stations (seed_ids)
merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
stream.merge()
print('Merged the following stations because of gaps:')
for merged_station in merged:
print(merged_station)
return stream, gaps
def check4gaps(data):
"""
check for gaps in Stream and remove them
@ -925,7 +948,10 @@ def get_stations(data):
def check4rotated(data, metadata=None, verbosity=1):
"""
Check all traces in data. If a trace is not in ZNE rotation (last symbol of channel code is numeric) and the trace
is in the metadata with azimuth and dip, rotate it to classical ZNE orientation.
Rotating the traces requires them to be of the same length, so, all traces will be trimmed to a common length as a
side effect.
:param data: stream object containing seismic traces
:type data: `~obspy.core.stream.Stream`
:param metadata: tuple containing metadata type string and metadata parser object
@ -942,100 +968,59 @@ def check4rotated(data, metadata=None, verbosity=1):
Azimut and dip are fetched from metadata. To be rotated, traces of a station have to be cut to the same length.
Returns unrotated traces of no metadata is provided
:param wfstream: stream containing seismic traces
:param wfstream: stream containing seismic traces of a station
:type wfstream: `~obspy.core.stream.Stream`
:param metadata: tuple containing metadata type string and metadata parser object
:type metadata: (str, `~obspy.io.xseed.parser.Parser`)
:return: stream object with traditionally oriented traces (ZNE)
:rtype: `~obspy.core.stream.Stream`
"""
try:
# indexing fails if metadata is None
metadata[0]
except TypeError:
if verbosity:
msg = 'Warning: could not rotate traces since no metadata was given\nset Inventory file!'
print(msg)
return wfstream
if metadata[0] is None:
# sometimes metadata is (None, (None,))
if verbosity:
msg = 'Warning: could not rotate traces since no metadata was given\nCheck inventory directory!'
print(msg)
return wfstream
else:
parser = metadata[1]
def get_dip_azimut(parser, trace_id):
"""
Gets azimuth and dip by trace id out of the metadata parser
:param parser: metadata parser object
:type parser: `~obspy.io.xseed.parser.Parser`
:param trace_id: eg. 'BW.RJOB..EHZ',
:type trace_id: str
:return: tuple containing dip and azimuth of the trace corresponding to trace_id
:rtype: (float, float)
"""
dip = None
azimut = None
try:
blockettes = parser._select(trace_id)
except SEEDParserException as e:
print(e)
raise ValueError
for blockette_ in blockettes:
if blockette_.id != 52:
continue
dip = blockette_.dip
azimut = blockette_.azimuth
break
if (dip is None or azimut is None) or (dip == 0 and azimut == 0):
error_msg = 'Dip and azimuth not available for trace_id {}'.format(trace_id)
raise ValueError(error_msg)
return dip, azimut
# check if any traces in this station need to be rotated
trace_ids = [trace.id for trace in wfstream]
for trace_id in trace_ids:
orientation = trace_id[-1] # last letter if trace id is orientation code, ZNE or 123
if orientation.isnumeric():
# misaligned channels have a number as orientation
azimuts = []
dips = []
for trace_id in trace_ids:
try:
dip, azimut = get_dip_azimut(parser, trace_id)
except ValueError as e:
print(e)
print('Failed to rotate station {}, no azimuth or dip available in metadata'.format(trace_id))
return wfstream
azimuts.append(azimut)
dips.append(dip)
# to rotate all traces must have same length
wfstream = trim_station_components(wfstream, trim_start=True, trim_end=True)
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 station {} to ZNE'.format(trace_id))
z_index = dips.index(min(dips)) # get z-trace index (dip is measured from 0 to -90)
wfstream[z_index].data = z
wfstream[z_index].stats.channel = wfstream[z_index].stats.channel[0:-1] + 'Z'
del trace_ids[z_index]
for trace_id in trace_ids:
dip, az = get_dip_azimut(parser, trace_id)
trace = wfstream.select(id=trace_id)[0]
if az > 315 and az <= 45 or az > 135 and az <= 225:
trace.data = n
trace.stats.channel = trace.stats.channel[0:-1] + 'N'
elif az > 45 and az <= 135 or az > 225 and az <= 315:
trace.data = e
trace.stats.channel = trace.stats.channel[0:-1] + 'E'
break
else:
continue
orientations = [trace_id[-1] for trace_id in trace_ids]
rotation_required = [orientation.isnumeric() for orientation in orientations]
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]
except (KeyError, TypeError) as e:
print('Failed to rotate trace {}, no azimuth or dip available in metadata'.format(trace_id))
return wfstream
if len(wfstream) < 3:
print('Failed to rotate Stream {}, not enough components available.'.format(wfstream))
return wfstream
# to rotate all traces must have same length, so trim them
wfstream = trim_station_components(wfstream, trim_start=True, trim_end=True)
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))
# 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)
wfstream[z_index].data = z
wfstream[z_index].stats.channel = wfstream[z_index].stats.channel[0:-1] + 'Z'
del trace_ids[z_index]
for trace_id in trace_ids:
coordinates = metadata.get_coordinates(trace_id, t_start)
dip, az = coordinates['dip'], coordinates['azimuth']
trace = wfstream.select(id=trace_id)[0]
if az > 315 or az <= 45 or az > 135 and az <= 225:
trace.data = n
trace.stats.channel = trace.stats.channel[0:-1] + 'N'
elif az > 45 and az <= 135 or az > 225 and az <= 315:
trace.data = e
trace.stats.channel = trace.stats.channel[0:-1] + 'E'
return wfstream
if metadata is None:
if verbosity:
msg = 'Warning: could not rotate traces since no metadata was given\nset Inventory file!'
print(msg)
return data
stations = get_stations(data)
for station in stations: # loop through all stations and rotate data if neccessary
wf_station = data.select(station=station)
rotate_components(wf_station, metadata)
@ -1093,7 +1078,7 @@ def runProgram(cmd, parameter=None):
subprocess.check_output('{} | tee /dev/stderr'.format(cmd), shell=True)
def which(program, infile=None):
def which(program, parameter):
"""
takes a program name and returns the full path to the executable or None
modified after: http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python
@ -1108,16 +1093,9 @@ def which(program, infile=None):
for key in settings.allKeys():
if 'binPath' in key:
os.environ['PATH'] += ':{0}'.format(settings.value(key))
if infile is None:
# use default parameter-file name
bpath = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
else:
bpath = os.path.join(os.path.expanduser('~'), '.pylot', infile)
if os.path.exists(bpath):
nllocpath = ":" + PylotParameter(bpath).get('nllocbin')
os.environ['PATH'] += nllocpath
except ImportError as e:
nllocpath = ":" + parameter.get('nllocbin')
os.environ['PATH'] += nllocpath
except Exception as e:
print(e.message)
def is_exe(fpath):
@ -1155,7 +1133,7 @@ def loopIdentifyPhase(phase):
"""
from pylot.core.util.defaults import ALTSUFFIX
if phase == None:
if phase is None:
raise NameError('Can not identify phase that is None')
phase_copy = phase
@ -1205,20 +1183,6 @@ def identifyPhaseID(phase):
return identifyPhase(loopIdentifyPhase(phase))
def has_spe(pick):
"""
Check for 'spe' key (symmetric picking error) in dict and return its value if found, else return None
:param pick: pick dictionary
:type pick: dict
:return: value of 'spe' key
:rtype: float or None
"""
if not 'spe' in pick.keys():
return None
else:
return pick['spe']
def check_all_obspy(eventlist):
ev_type = 'obspydmt'
return check_event_folders(eventlist, ev_type)
@ -1251,8 +1215,8 @@ def check_event_folder(path):
folder = path.split('/')[-1]
# for pylot: select only folders that start with 'e', containin two dots and have length 12
if (folder.startswith('e')
and len(folder.split('.')) == 3
and len(folder) == 12):
and len(folder.split('.')) == 3
and len(folder) == 12):
ev_type = 'pylot'
elif check_obspydmt_eventfolder(folder)[0]:
ev_type = 'obspydmt'
@ -1263,3 +1227,56 @@ if __name__ == "__main__":
import doctest
doctest.testmod()
class SetChannelComponents(object):
def __init__(self):
self.setDefaultCompPosition()
def setDefaultCompPosition(self):
# default component order
self.compPosition_Map = dict(Z=2, N=1, E=0)
self.compName_Map = {'3': 'Z',
'1': 'N',
'2': 'E'}
def _getCurrentPosition(self, component):
for key, value in self.compName_Map.items():
if value == component:
return key, value
errMsg = 'getCurrentPosition: Could not find former position of component {}.'.format(component)
raise ValueError(errMsg)
def _switch(self, component, component_alter):
# Without switching, multiple definitions of the same alter_comp are possible
old_alter_comp, _ = self._getCurrentPosition(component)
old_comp = self.compName_Map[component_alter]
if not old_alter_comp == component_alter and not old_comp == component:
self.compName_Map[old_alter_comp] = old_comp
print('switch: Automatically switched component {} to {}'.format(old_alter_comp, old_comp))
def setCompPosition(self, component_alter, component, switch=True):
component_alter = str(component_alter)
if not component_alter in self.compName_Map.keys():
errMsg = 'setCompPosition: Unrecognized alternative component {}. Expecting one of {}.'
raise ValueError(errMsg.format(component_alter, self.compName_Map.keys()))
if not component in self.compPosition_Map.keys():
errMsg = 'setCompPosition: Unrecognized target component {}. Expecting one of {}.'
raise ValueError(errMsg.format(component, self.compPosition_Map.keys()))
print('setCompPosition: set component {} to {}'.format(component_alter, component))
if switch:
self._switch(component, component_alter)
self.compName_Map[component_alter] = component
def getCompPosition(self, component):
return self._getCurrentPosition(component)[0]
def getPlotPosition(self, component):
component = str(component)
if component in self.compPosition_Map.keys():
return self.compPosition_Map[component]
elif component in self.compName_Map.keys():
return self.compPosition_Map[self.compName_Map[component]]
else:
errMsg = 'getCompPosition: Unrecognized component {}. Expecting one of {} or {}.'
raise ValueError(errMsg.format(component, self.compPosition_Map.keys(), self.compName_Map.keys()))

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

File diff suppressed because it is too large Load Diff

6
pylot/styles/bright.qss
View File

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

6
pylot/styles/dark.qss
View File

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

15
pylot/styles/style_settings.py
View File

@ -5,18 +5,18 @@
# the base color
phasecolors = {
'manual': {
'P':{
'P': {
'rgba': (0, 0, 255, 255),
'modifier': 'g'},
'S':{
'S': {
'rgba': (255, 0, 0, 255),
'modifier': 'b'}
},
'auto':{
'P':{
'auto': {
'P': {
'rgba': (140, 0, 255, 255),
'modifier': 'g'},
'S':{
'S': {
'rgba': (255, 140, 0, 255),
'modifier': 'b'}
}
@ -24,8 +24,8 @@ phasecolors = {
# Set plot colors and stylesheet for each style
stylecolors = {
'default':{
'linecolor':{
'default': {
'linecolor': {
'rgba': (0, 0, 0, 255)},
'background': {
'rgba': (255, 255, 255, 255)},
@ -67,4 +67,3 @@ stylecolors = {
'filename': 'bright.qss'}
}
}

2
setup.py
View File

@ -8,7 +8,7 @@ setup(
packages=['pylot', 'pylot.core', 'pylot.core.loc', 'pylot.core.pick',
'pylot.core.io', 'pylot.core.util', 'pylot.core.active',
'pylot.core.analysis', 'pylot.testing'],
requires=['obspy', 'PySide', 'matplotlib', 'numpy'],
requires=['obspy', 'PySide', 'matplotlib', 'numpy', 'scipy', 'pyqtgraph'],
url='dummy',
license='LGPLv3',
author='Sebastian Wehling-Benatelli',

27
tests/__init__.py Normal file
View File

@ -0,0 +1,27 @@
# -*- coding: utf-8 -*-
# --------------------------------------------------------
# Purpose: Convience imports for PyLoT
#
'''
================================================
PyLoT - the Python picking and Localization Tool
================================================
This python library contains a graphical user interfaces for picking
seismic phases. This software needs ObsPy (http://github.com/obspy/obspy/wiki)
and the Qt4 libraries to be installed first.
PILOT has been developed in Mathworks' MatLab. In order to distribute
PILOT without facing portability problems, it has been decided to re-
develop the software package in Python. The great work of the ObsPy
group allows easy handling of a bunch of seismic data and PyLoT will
benefit a lot compared to the former MatLab version.
The development of PyLoT is part of the joint research project MAGS2.
:copyright:
The PyLoT Development Team
:license:
GNU Lesser General Public License, Version 3
(http://www.gnu.org/copyleft/lesser.html)
'''

334
tests/test_Metadata/test_Metadata.py Normal file
View File

@ -0,0 +1,334 @@
import unittest
import os
from obspy import UTCDateTime
from obspy.io.xseed.utils import SEEDParserException
from obspy.io.xseed import Parser
from pylot.core.util.dataprocessing import Metadata
from tests.utils import HidePrints
class TestMetadata(unittest.TestCase):
def setUp(self):
self.station_id = 'BW.WETR..HH'
self.time = UTCDateTime('2012-08-01')
metadata_folder = os.path.join('test_data', 'dless_multiple_files', 'metadata1')
self.m = Metadata(metadata_folder)
def test_get_coordinates_sucess(self):
expected = {'Z': {u'elevation': 607.0, u'longitude': 12.87571, u'local_depth': 0.0, u'azimuth': 0.0,
u'latitude': 49.14502, u'dip': -90.0},
'E': {u'azimuth': 90.0, u'dip': 0.0, u'elevation': 607.0, u'latitude': 49.14502,
u'local_depth': 0.0, u'longitude': 12.87571},
'N': {u'azimuth': 0.0, u'dip': 0.0, u'elevation': 607.0, u'latitude': 49.14502, u'local_depth': 0.0,
u'longitude': 12.87571}
}
result = {}
for channel in ('Z', 'N', 'E'):
with HidePrints():
coords = self.m.get_coordinates(self.station_id+channel, time=self.time)
result[channel] = coords
self.assertDictEqual(result[channel], expected[channel])
def test_get_coordinates_sucess_no_time(self):
expected = {'Z': {u'elevation': 607.0, u'longitude': 12.87571, u'local_depth': 0.0, u'azimuth': 0.0,
u'latitude': 49.14502, u'dip': -90.0},
'E': {u'azimuth': 90.0, u'dip': 0.0, u'elevation': 607.0, u'latitude': 49.14502,
u'local_depth': 0.0, u'longitude': 12.87571},
'N': {u'azimuth': 0.0, u'dip': 0.0, u'elevation': 607.0, u'latitude': 49.14502, u'local_depth': 0.0,
u'longitude': 12.87571}
}
result = {}
for channel in ('Z', 'N', 'E'):
with HidePrints():
coords = self.m.get_coordinates(self.station_id+channel)
result[channel] = coords
self.assertDictEqual(result[channel], expected[channel])
class TestMetadataAdding(unittest.TestCase):
"""Tests if adding files and directories to a metadata object works."""
def setUp(self):
self.station_id = 'BW.WETR..HH'
self.metadata_folders = (os.path.join('test_data', 'dless_multiple_files', 'metadata1'),
os.path.join('test_data', 'dless_multiple_files', 'metadata2'))
self.m = Metadata()
def test_add_inventory_folder(self):
"""Test if add_inventory adds the folder to the list of inventories"""
self.m.add_inventory(self.metadata_folders[0])
# adding an inventory folder should append it to the list of inventories
self.assertDictEqual({}, self.m.inventory_files)
self.assertDictEqual({}, self.m.seed_ids)
self.assertEqual([self.metadata_folders[0]], self.m.inventories)
def test_add_inventory_file(self):
"""Test if add_inventory_file adds the folder containing the file to the list of inventories and
if the files is added to inventory_files"""
fpath = os.path.join(self.metadata_folders[0], 'DATALESS.BW.WETR..HHZ')
self.m.add_inventory_file(fpath)
# adding an inventory file should append its folder to the list of inventories and the file to the
self.assertEqual([os.path.join(self.metadata_folders[0], 'DATALESS.BW.WETR..HHZ')],
self.m.inventory_files.keys()) # does the filename exist in inventory files?
self.assertEqual(['data', 'invtype'], self.m.inventory_files[os.path.join(self.metadata_folders[0],
'DATALESS.BW.WETR..HHZ')].keys()) # is the required information attacht to the filename?
self.assertDictEqual({}, self.m.seed_ids)
self.assertEqual([self.metadata_folders[0]], self.m.inventories)
def test_add_inventory_invalid_path(self):
"""Test if adding an inventory that is not an existing directory fails with an exception"""
with self.assertRaises(Exception):
self.m.add_inventory('InvalidDirName')
self.assertEqual([], self.m.inventories) # inventory list should still be empty
def test_add_inventory_file_invalid_path(self):
"""Test if adding a inventory file with an invalid path fails with an exception"""
with self.assertRaises(Exception):
self.m.add_inventory_file('/invalid/file/name')
self.assertEqual([], self.m.inventories) # inventory list should still be empty
class TestMetadataRemoval(unittest.TestCase):
"""Tests if removing files and directories to a metadata object works."""
def setUp(self):
self.station_id = 'BW.WETR..HH'
self.metadata_folders = (os.path.join('test_data', 'dless_multiple_files', 'metadata1'),
os.path.join('test_data', 'dless_multiple_files', 'metadata2'))
self.m = Metadata()
def test_remove_all_inventories(self):
"""Test if function remove_inventory cleans the Metadata object """
# add multiple inventories
for folder in self.metadata_folders:
self.m.add_inventory(folder)
self.m.remove_all_inventories()
self.isEmpty(self.m)
def test_remove_inventory(self):
"""Test if remove_inventory removes single inventories"""
# add multiple inventories
for folder in self.metadata_folders:
self.m.add_inventory(folder)
self.m.remove_inventory(self.metadata_folders[0])
self.assertNotIn(self.metadata_folders[0], self.m.inventories)
self.m.remove_inventory(self.metadata_folders[1])
self.assertNotIn(self.metadata_folders[1], self.m.inventories)
self.isEmpty(self.m)
def test_remove_inventory_not_in_inventory_list(self):
"""Test if remove_inventory does not modify the metadata instance if the given inventory to remove does not
exist in the instance."""
# add multiple inventories
self.m.add_inventory(self.metadata_folders[0])
with HidePrints():
self.m.remove_inventory('metadata_not_existing')
self.assertIn(self.metadata_folders[0], self.m.inventories)
def isEmpty(self, metadata):
"""Asserts if the given metadata object is empty"""
self.assertDictEqual({}, metadata.inventory_files)
self.assertDictEqual({}, metadata.seed_ids)
self.assertEqual([], metadata.inventories)
class TestMetadata_read_single_file(unittest.TestCase):
def setUp(self):
self.station_id = 'BW.WETR..HHZ'
self.metadata_folders = (os.path.join('test_data', 'dless_multiple_files', 'metadata1'),
os.path.join('test_data', 'dless_multiple_files', 'metadata2'))
self.metadata_paths = []
self.m = Metadata()
def test_read_single_file(self):
"""Test if reading a single file works"""
fname = os.path.join(self.metadata_folders[0], 'DATALESS.'+self.station_id)
with HidePrints():
res = self.m.read_single_file(fname)
# method should return true if file is successfully read
self.assertTrue(res)
# list of inventories (folders) should be empty
self.assertEqual([], self.m.inventories)
# list of inventory files should contain the added file
self.assertIn(fname, self.m.inventory_files.keys())
self.assertEqual({}, self.m.seed_ids)
def test_read_single_file_invalid_path(self):
"""Test if reading from a non existing file fails. The filename should not be
added to the metadata object"""
fname = os.path.join("this", "path", "doesnt", "exist")
with HidePrints():
res = self.m.read_single_file(fname)
# method should return None if file reading fails
self.assertIsNone(res)
# list of inventories (folders) should be empty
self.assertEqual([], self.m.inventories)
# list of inventory files should not contain the added file
self.assertNotIn(fname, self.m.inventory_files.keys())
self.assertEqual({}, self.m.seed_ids)
def test_read_single_file_multiple_times(self):
"""Test if reading a file twice doesnt add it twice to the metadata object"""
fname = os.path.join(self.metadata_folders[0], 'DATALESS.'+self.station_id)
with HidePrints():
res1 = self.m.read_single_file(fname)
res2 = self.m.read_single_file(fname)
self.assertTrue(res1)
self.assertIsNone(res2)
self.assertItemsEqual([fname], self.m.inventory_files.keys())
class TestMetadataMultipleTime(unittest.TestCase):
"""Test if stations with multiple metadata entries in a single file are handled correctly.
The user must specify the time where he wants to get metadata.
The station ROTT changed has metadata available at multiple times
LE.ROTT..HNE | 200.00 Hz | Titan 4g-EDR-209, Very Low gain, 200 sps | 2015-01-08 - 2015-03-19 | Lat: 49.1, Lng: 8.1
LE.ROTT..HNE | 200.00 Hz | Titan 4g-EDR-209, Very Low gain, 200 sps | 2015-03-19 - | Lat: 49.1, Lng: 8.1
LE.ROTT..HNN | 200.00 Hz | Titan 4g-EDR-209, Very Low gain, 200 sps | 2015-01-08 - 2015-03-19 | Lat: 49.1, Lng: 8.1
LE.ROTT..HNN | 200.00 Hz | Titan 4g-EDR-209, Very Low gain, 200 sps | 2015-03-19 - | Lat: 49.1, Lng: 8.1
LE.ROTT..HNZ | 200.00 Hz | Titan 4g-EDR-209, Very Low gain, 200 sps | 2015-01-08 - 2015-03-19 | Lat: 49.1, Lng: 8.1
LE.ROTT..HNZ | 200.00 Hz | Titan 4g-EDR-209, Very Low gain, 200 sps | 2015-03-19 - | Lat: 49.1, Lng: 8.1
"""
def setUp(self):
self.seed_id = 'LE.ROTT..HN'
path = os.path.dirname(__file__) # gets path to currently running script
metadata = os.path.join('test_data', 'dless_multiple_times', 'MAGS2_LE_ROTT.dless') # specific subfolder of test data
metadata_path = os.path.join(path, metadata)
self.m = Metadata(metadata_path)
self.p = Parser(metadata_path)
def test_get_metadata_works_without_datetime(self):
"""Test if get_metadata works if multiple metadata entries are available but no time is
specified."""
for channel in ('Z', 'N', 'E'):
with HidePrints():
md = self.m.get_metadata(self.seed_id + channel)
self.assertDictEqual(md['data'].get_inventory(), self.p.get_inventory())
def test_get_metadata_works_with_first_datetime(self):
"""Test if get_metadata works if multiple metadata entries are available and the older time is specified."""
t = UTCDateTime('2015-02-08')
for channel in ('Z', 'N', 'E'):
with HidePrints():
md = self.m.get_metadata(self.seed_id + channel, t)
self.assertDictEqual(md['data'].get_inventory(), self.p.get_inventory())
def test_get_metadata_fails_when_time_before_starttime(self):
"""Tests if get_metadata returns None when given a data that is before the start date
of the metadata"""
with HidePrints():
md = self.m.get_metadata(self.seed_id, UTCDateTime('1960-07-20'))
self.assertIs(md, None)
def test_get_metadata_invalid_seed_id(self):
"""Tes if get metadata returns none when asked for a seed id that does not exist"""
with HidePrints():
res = self.m.get_metadata("this.doesnt..exist")
self.assertIsNone(res)
class TestMetadataMultipleEntries(unittest.TestCase):
"""
The station KB.TMO07 has changed instruments multiple times.
Networks:
KB (KB network)
Stations:
KB.TMO07 (Karlsruhe GPI)
Channels:
KB.TMO07.00.BHE | 50.00 Hz | Streckeisen KABBA-STS-2 | 2004-12-06 - 2005-04-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHE | 50.00 Hz | Streckeisen KABBA-STS-2 | 2005-04-18 - 2006-07-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHE | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2006-10-10 - 2006-11-14 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHE | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2006-11-24 - 2007-01-12 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHE | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-01-18 - 2007-03-15 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHE | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-10-25 - 2007-11-21 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHE | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-11-21 - 2008-01-17 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Streckeisen KABBA-STS-2 | 2004-12-06 - 2005-04-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Streckeisen KABBA-STS-2 | 2005-04-18 - 2006-07-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2006-10-10 - 2006-11-14 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2006-11-24 - 2007-01-12 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-01-18 - 2007-03-15 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-10-25 - 2007-11-21 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHN | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-11-21 - 2008-01-17 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Streckeisen KABBA-STS-2 | 2004-12-06 - 2005-04-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Streckeisen KABBA-STS-2 | 2005-04-18 - 2006-07-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2006-10-10 - 2006-11-14 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2006-11-24 - 2007-01-12 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-01-18 - 2007-03-15 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-10-25 - 2007-11-21 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.BHZ | 50.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-11-21 - 2008-01-17 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-01-12 - 2007-01-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-10-10 - 2007-10-25 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Streckeisen KABBA-STS-2 | 2008-07-11 - 2008-12-05 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Streckeisen KABBA-STS-2 | 2009-05-12 - 2010-02-15 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Streckeisen KABBA-STS-2 | 2010-02-15 - 2010-04-07 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Lennartz KABBA-LE-3D/1 | 2010-04-07 - 2010-08-03 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 200.00 Hz | Streckeisen KABBA-STS-2 | 2010-08-05 - 2010-12-20 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 100.00 Hz | Streckeisen KABBA-STS-2 | 2010-12-20 - 2010-12-22 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 200.00 Hz | Streckeisen KABBA-STS-2 | 2010-12-22 - 2011-04-02 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 200.00 Hz | Streckeisen KABBA-STS-2 | 2011-04-15 - 2012-05-07 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHE | 200.00 Hz | Streckeisen KABBA-STS-2 | 2012-05-07 - | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-01-12 - 2007-01-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-10-10 - 2007-10-25 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Streckeisen KABBA-STS-2 | 2008-07-11 - 2008-12-05 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Streckeisen KABBA-STS-2 | 2009-05-12 - 2010-02-15 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Streckeisen KABBA-STS-2 | 2010-02-15 - 2010-04-07 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Lennartz KABBA-LE-3D/1 | 2010-04-07 - 2010-08-03 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 200.00 Hz | Streckeisen KABBA-STS-2 | 2010-08-05 - 2010-12-20 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 100.00 Hz | Streckeisen KABBA-STS-2 | 2010-12-20 - 2010-12-22 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 200.00 Hz | Streckeisen KABBA-STS-2 | 2010-12-22 - 2011-04-02 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 200.00 Hz | Streckeisen KABBA-STS-2 | 2011-04-15 - 2012-05-07 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHN | 200.00 Hz | Streckeisen KABBA-STS-2 | 2012-05-07 - | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-01-12 - 2007-01-18 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Lennartz KABBA-LE-3D/5 | 2007-10-10 - 2007-10-25 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Streckeisen KABBA-STS-2 | 2008-07-11 - 2008-12-05 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Streckeisen KABBA-STS-2 | 2009-05-12 - 2010-02-15 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Streckeisen KABBA-STS-2 | 2010-02-15 - 2010-04-07 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Lennartz KABBA-LE-3D/1 | 2010-04-07 - 2010-08-03 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 200.00 Hz | Streckeisen KABBA-STS-2 | 2010-08-05 - 2010-12-20 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 100.00 Hz | Streckeisen KABBA-STS-2 | 2010-12-20 - 2010-12-22 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 200.00 Hz | Streckeisen KABBA-STS-2 | 2010-12-22 - 2011-04-02 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 200.00 Hz | Streckeisen KABBA-STS-2 | 2011-04-15 - 2012-05-07 | Lat: 49.0, Lng: 8.4
KB.TMO07.00.HHZ | 200.00 Hz | Streckeisen KABBA-STS-2 | 2012-05-07 - | Lat: 49.0, Lng: 8.4
"""
def setUp(self):
self.seed_id = 'KB.TMO07.00.HHZ'
path = os.path.dirname(__file__) # gets path to currently running script
metadata = os.path.join('test_data', 'dless_multiple_instruments', 'MAGS2_KB_TMO07.dless') # specific subfolder of test data
metadata_path = os.path.join(path, metadata)
self.m = Metadata(metadata_path)
self.p = Parser(metadata_path)
def test_get_paz_current_time(self):
"""Test if getting the paz from the metadata object with the current time works"""
t = UTCDateTime()
with HidePrints():
pazm = self.m.get_paz(self.seed_id, t)
pazp = self.p.get_paz(self.seed_id, t)
self.assertEqual(pazm, pazp)
def test_get_paz_past(self):
"""Test if getting paz from metadata object with a time in the past works"""
t = UTCDateTime('2007-01-13')
with HidePrints():
pazm = self.m.get_paz(self.seed_id, t)
pazp = self.p.get_paz(self.seed_id, t)
self.assertEqual(pazm, pazp)
def test_get_paz_time_not_exisiting(self):
"""Test if getting paz from metadata at a time where there is no metadata
available fails correctly"""
with self.assertRaises(SEEDParserException):
with HidePrints():
self.m.get_paz(self.seed_id, UTCDateTime('1990-1-1'))
def test_get_paz_seed_id_not_existing(self):
"""Test if getting paz from a non existing seed id returns None as expected."""
with HidePrints():
res = self.m.get_paz('This.doesnt..exist', UTCDateTime)
self.assertIsNone(res)

1
tests/test_Metadata/test_data/dless_multiple_files/metadata1/DATALESS.BW.WETR..HHE Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata1/DATALESS.BW.WETR..HHN Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata1/DATALESS.BW.WETR..HHZ Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata2/DATALESS.GR.GRA1..HHE Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata2/DATALESS.GR.GRA1..HHN Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata2/DATALESS.GR.GRA1..HHZ Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata2/DATALESS.GR.GRA2..HHE Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata2/DATALESS.GR.GRA2..HHN Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_files/metadata2/DATALESS.GR.GRA2..HHZ Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_instruments/MAGS2_KB_TMO07.dless Normal file
View File

File diff suppressed because one or more lines are too long

1
tests/test_Metadata/test_data/dless_multiple_times/MAGS2_LE_ROTT.dless Normal file
View File

File diff suppressed because one or more lines are too long

52
tests/utils.py Normal file
View File

@ -0,0 +1,52 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Utilities/helpers for testing"""
import sys
import os
class HidePrints:
"""
Context manager that hides all standard output within its body.
The optional hide_prints argument can be used to quickly enable printing during debugging of tests.
Use (This will result in all console output of noisy_function to be suppressed):
from tests.utils import HidePrints
with HidePrints():
noise_function()
"""
@staticmethod
def hide(func, *args, **kwargs):
"""Decorator that hides all prints of the decorated function.
Use:
from tests.utils import HidePrints
@HidePrints.hide
def noise()
print("NOISE")
"""
def silencer(*args, **kwargs):
with HidePrints():
func(*args, **kwargs)
return silencer
def __init__(self, hide_prints=True):
"""Create object with hide_prints=False to disable print hiding"""
self.hide = hide_prints
def __enter__(self):
"""Redirect stdout to /dev/null, save old stdout"""
if self.hide:
self._original_stdout = sys.stdout
devnull = open(os.devnull, "w")
sys.stdout = devnull
def __exit__(self, exc_type, exc_val, exc_tb):
"""Reinstate old stdout"""
if self.hide:
sys.stdout = self._original_stdout