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

Cosmetics, changed print commands to keep compatibility to Python 3.

Browse Source
This commit is contained in:
Ludger Küperkoch 2015-09-03 14:55:25 +02:00
parent bf1194ec3b
commit 60b9f176f0
1 changed files with 43 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -7,7 +7,6 @@
:author: Ludger Kueperkoch / MAGS2 EP3 working group :author: Ludger Kueperkoch / MAGS2 EP3 working group
""" """
import numpy as np import numpy as np
import scipy as sc import scipy as sc
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -44,7 +43,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
LPick = None LPick = None
EPick = None EPick = None
PickError = None PickError = None
print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...' print ("earllatepicker: Get earliest and latest possible pick relative to most likely pick ...")
x = X[0].data x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
@ -60,8 +59,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
ilup, = np.where(x[isignal] > nlevel) ilup, = np.where(x[isignal] > nlevel)
ildown, = np.where(x[isignal] < -nlevel) ildown, = np.where(x[isignal] < -nlevel)
if not ilup.size and not ildown.size: if not ilup.size and not ildown.size:
print 'earllatepicker: Signal lower than noise level!' print ("earllatepicker: Signal lower than noise level!")
print 'Skip this trace!' print ("Skip this trace!")
return LPick, EPick, PickError return LPick, EPick, PickError
il = min(np.min(ilup) if ilup.size else float('inf'), il = min(np.min(ilup) if ilup.size else float('inf'),
np.min(ildown) if ildown.size else float('inf')) np.min(ildown) if ildown.size else float('inf'))
@ -143,7 +142,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
FM = None FM = None
if Pick is not None: if Pick is not None:
print 'fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...' print ("fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...")
xraw = Xraw[0].data xraw = Xraw[0].data
xfilt = Xfilt[0].data xfilt = Xfilt[0].data
@ -182,15 +181,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
else: else:
li1 = index1[0] li1 = index1[0]
if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 0: if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 0:
print 'fmpicker: Onset on unfiltered trace too emergent for first motion determination!' print ("fmpicker: Onset on unfiltered trace too emergent for first motion determination!")
P1 = None P1 = None
else: else:
imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]])) imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]]))
if imax1 == 0: if imax1 == 0:
imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][index1[1]]])) imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][index1[1]]]))
if imax1 == 0: if imax1 == 0:
print 'fmpicker: Zero crossings too close!' print ("fmpicker: Zero crossings too close!")
print 'Skip first motion determination!' print ("Skip first motion determination!")
return FM return FM
islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1])) islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1]))
@ -224,15 +223,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
else: else:
li2 = index2[0] li2 = index2[0]
if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 0: if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 0:
print 'fmpicker: Onset on filtered trace too emergent for first motion determination!' print ("fmpicker: Onset on filtered trace too emergent for first motion determination!")
P2 = None P2 = None
else: else:
imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]])) imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]]))
if imax2 == 0: if imax2 == 0:
imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][index2[1]]])) imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][index2[1]]]))
if imax2 == 0: if imax2 == 0:
print 'fmpicker: Zero crossings too close!' print ("fmpicker: Zero crossings too close!")
print 'Skip first motion determination!' print ("Skip first motion determination!")
return FM return FM
islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2])) islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2]))
@ -256,7 +255,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
elif P1[0] > 0 and P2[0] <= 0: elif P1[0] > 0 and P2[0] <= 0:
FM = '+' FM = '+'
print 'fmpicker: Found polarity %s' % FM print ("fmpicker: Found polarity %s" % FM)
if iplot > 1: if iplot > 1:
plt.figure(iplot) plt.figure(iplot)
@ -331,10 +330,10 @@ def getSNR(X, TSNR, t1):
# get signal window # get signal window
isignal = getsignalwin(t, t1, TSNR[2]) isignal = getsignalwin(t, t1, TSNR[2])
if np.size(inoise) < 1: if np.size(inoise) < 1:
print 'getSNR: Empty array inoise, check noise window!' print ("getSNR: Empty array inoise, check noise window!")
return return
elif np.size(isignal) < 1: elif np.size(isignal) < 1:
print 'getSNR: Empty array isignal, check signal window!' print ("getSNR: Empty array isignal, check signal window!")
return return
# demean over entire waveform # demean over entire waveform
@ -372,7 +371,7 @@ def getnoisewin(t, t1, tnoise, tgap):
inoise, = np.where((t <= max([t1 - tgap, 0])) \ inoise, = np.where((t <= max([t1 - tgap, 0])) \
& (t >= max([t1 - tnoise - tgap, 0]))) & (t >= max([t1 - tnoise - tgap, 0])))
if np.size(inoise) < 1: if np.size(inoise) < 1:
print 'getnoisewin: Empty array inoise, check noise window!' print ("getnoisewin: Empty array inoise, check noise window!")
return inoise return inoise
@ -396,7 +395,7 @@ def getsignalwin(t, t1, tsignal):
isignal, = np.where((t <= min([t1 + tsignal, len(t)])) \ isignal, = np.where((t <= min([t1 + tsignal, len(t)])) \
& (t >= t1)) & (t >= t1))
if np.size(isignal) < 1: if np.size(isignal) < 1:
print 'getsignalwin: Empty array isignal, check signal window!' print ("getsignalwin: Empty array isignal, check signal window!")
return isignal return isignal
@ -483,8 +482,8 @@ def wadaticheck(pickdic, dttolerance, iplot):
# calculate vp/vs ratio before check # calculate vp/vs ratio before check
vpvsr = p1[0] + 1 vpvsr = p1[0] + 1
print '###############################################' print ("###############################################")
print 'wadaticheck: Average Vp/Vs ratio before check:', vpvsr print ("wadaticheck: Average Vp/Vs ratio before check:", vpvsr)
checkedPpicks = [] checkedPpicks = []
checkedSpicks = [] checkedSpicks = []
@ -521,18 +520,18 @@ def wadaticheck(pickdic, dttolerance, iplot):
# calculate vp/vs ratio after check # calculate vp/vs ratio after check
cvpvsr = p2[0] + 1 cvpvsr = p2[0] + 1
print 'wadaticheck: Average Vp/Vs ratio after check:', cvpvsr print ("wadaticheck: Average Vp/Vs ratio after check:", cvpvsr)
print 'wadatacheck: Skipped %d S pick(s).' % ibad print ("wadatacheck: Skipped %d S pick(s)." % ibad)
else: else:
print '###############################################' print ("###############################################")
print 'wadatacheck: Not enough checked S-P times available!' print ("wadatacheck: Not enough checked S-P times available!")
print 'Skip Wadati check!' print ("Skip Wadati check!")
checkedonsets = pickdic checkedonsets = pickdic
else: else:
print 'wadaticheck: Not enough S-P times available for reliable regression!' print ("wadaticheck: Not enough S-P times available for reliable regression!")
print 'Skip wadati check!' print ("Skip wadati check!")
wfitflag = 1 wfitflag = 1
# plot results # plot results
@ -592,7 +591,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
assert isinstance(X, Stream), "%s is not a stream object" % str(X) assert isinstance(X, Stream), "%s is not a stream object" % str(X)
print 'Checking signal length ...' print ("Checking signal length ...")
x = X[0].data x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
@ -601,8 +600,8 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
# generate envelope function from Hilbert transform # generate envelope function from Hilbert transform
y = np.imag(sc.signal.hilbert(x)) y = np.imag(sc.signal.hilbert(x))
e = np.sqrt(np.power(x, 2) + np.power(y, 2)) e = np.sqrt(np.power(x, 2) + np.power(y, 2))
# get noise window # get noise window in front of pick plus saftey gap
inoise = getnoisewin(t, pick, TSNR[0], TSNR[1]) inoise = getnoisewin(t, pick - 0.5, TSNR[0], TSNR[1])
# get signal window # get signal window
isignal = getsignalwin(t, pick, TSNR[2]) isignal = getsignalwin(t, pick, TSNR[2])
# calculate minimum adjusted signal level # calculate minimum adjusted signal level
@ -613,12 +612,12 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
numoverthr = len(np.where(e[isignal] >= minsiglevel)[0]) numoverthr = len(np.where(e[isignal] >= minsiglevel)[0])
if numoverthr >= minnum: if numoverthr >= minnum:
print 'checksignallength: Signal reached required length.' print ("checksignallength: Signal reached required length.")
returnflag = 1 returnflag = 1
else: else:
print 'checksignallength: Signal shorter than required minimum signal length!' print ("checksignallength: Signal shorter than required minimum signal length!")
print 'Presumably picked noise peak, pick is rejected!' print ("Presumably picked noise peak, pick is rejected!")
print '(min. signal length required:', minsiglength, 's)' print ("(min. signal length required:', minsiglength, 's)'")
returnflag = 0 returnflag = 0
if iplot == 2: if iplot == 2:
@ -629,7 +628,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
p2, = plt.plot(t[inoise], e[inoise]) p2, = plt.plot(t[inoise], e[inoise])
p3, = plt.plot(t[isignal],e[isignal], 'r') p3, = plt.plot(t[isignal],e[isignal], 'r')
p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \ p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \
[minsiglevel, minsiglevel], 'g') [minsiglevel, minsiglevel], 'g', linewidth=2)
p5, = plt.plot([pick, pick], [min(x), max(x)], 'b', linewidth=2) p5, = plt.plot([pick, pick], [min(x), max(x)], 'b', linewidth=2)
plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \ plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \
'Envelope Signal Window', 'Minimum Signal Level', \ 'Envelope Signal Window', 'Minimum Signal Level', \
@ -675,8 +674,8 @@ def checkPonsets(pickdic, dttolerance, iplot):
stations.append(key) stations.append(key)
# apply jackknife bootstrapping on variance of P onsets # apply jackknife bootstrapping on variance of P onsets
print '###############################################' print ("###############################################")
print 'checkPonsets: Apply jackknife bootstrapping on P-onset times ...' print ("checkPonsets: Apply jackknife bootstrapping on P-onset times ...")
[xjack,PHI_pseudo,PHI_sub] = jackknife(Ppicks, 'VAR', 1) [xjack,PHI_pseudo,PHI_sub] = jackknife(Ppicks, 'VAR', 1)
# get pseudo variances smaller than average variances # get pseudo variances smaller than average variances
# (times safety factor), these picks passed jackknife test # (times safety factor), these picks passed jackknife test
@ -684,7 +683,7 @@ def checkPonsets(pickdic, dttolerance, iplot):
# these picks did not pass jackknife test # these picks did not pass jackknife test
badjk = np.where(PHI_pseudo > 2 * xjack) badjk = np.where(PHI_pseudo > 2 * xjack)
badjkstations = np.array(stations)[badjk] badjkstations = np.array(stations)[badjk]
print 'checkPonsets: %d pick(s) did not pass jackknife test!' % len(badjkstations) print ("checkPonsets: %d pick(s) did not pass jackknife test!" % len(badjkstations))
# calculate median from these picks # calculate median from these picks
pmedian = np.median(np.array(Ppicks)[ij]) pmedian = np.median(np.array(Ppicks)[ij])
@ -696,9 +695,9 @@ def checkPonsets(pickdic, dttolerance, iplot):
goodstations = np.array(stations)[igood] goodstations = np.array(stations)[igood]
badstations = np.array(stations)[ibad] badstations = np.array(stations)[ibad]
print 'checkPonsets: %d pick(s) deviate too much from median!' % len(ibad) print ("checkPonsets: %d pick(s) deviate too much from median!" % len(ibad))
print 'checkPonsets: Skipped %d P pick(s) out of %d' % (len(badstations) \ print ("checkPonsets: Skipped %d P pick(s) out of %d" % (len(badstations) \
+ len(badjkstations), len(stations)) + len(badjkstations), len(stations)))
goodmarker = 'goodPonsetcheck' goodmarker = 'goodPonsetcheck'
badmarker = 'badPonsetcheck' badmarker = 'badPonsetcheck'
@ -765,8 +764,8 @@ def jackknife(X, phi, h):
g = len(X) / h g = len(X) / h
if type(g) is not int: if type(g) is not int:
print 'jackknife: Cannot divide quantity X in equal sized subgroups!' print ("jackknife: Cannot divide quantity X in equal sized subgroups!")
print 'Choose another size for subgroups!' print ("Choose another size for subgroups!")
return PHI_jack, PHI_pseudo, PHI_sub return PHI_jack, PHI_pseudo, PHI_sub
else: else:
# estimator of undisturbed spot check # estimator of undisturbed spot check
@ -834,7 +833,7 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
assert isinstance(X, Stream), "%s is not a stream object" % str(X) assert isinstance(X, Stream), "%s is not a stream object" % str(X)
print 'Check for spuriously picked S onset instead of P onset ...' print ("Check for spuriously picked S onset instead of P onset ...")
returnflag = 0 returnflag = 0
@ -875,9 +874,9 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
# vertical P-coda level must exceed horizontal P-coda level # vertical P-coda level must exceed horizontal P-coda level
# zfac times encodalevel # zfac times encodalevel
if zcodalevel < minsiglevel: if zcodalevel < minsiglevel:
print 'checkZ4S: Maybe S onset? Skip this P pick!' print ("checkZ4S: Maybe S onset? Skip this P pick!")
else: else:
print 'checkZ4S: P onset passes checkZ4S test!' print ("checkZ4S: P onset passes checkZ4S test!")
returnflag = 1 returnflag = 1
if iplot > 1: if iplot > 1: