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

changed whole autoPyLoT figure structure to fig/ax structure to keep connections to figures in memory, iPlot now part of autoPyLoT function call parameters (not yet in argparser)

Browse Source
This commit is contained in:
Marcel Paffrath 2017-05-08 15:38:41 +02:00
parent 6563b01293
commit c784d46521
7 changed files with 290 additions and 296 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
autoPyLoT.py
View File

@ -29,7 +29,7 @@ from pylot.core.util.version import get_git_version as _getVersionString
__version__ = _getVersionString()
def autoPyLoT(inputfile, fnames=None, savepath=None):
def autoPyLoT(inputfile, fnames=None, savepath=None, iplot=0):
"""
Determine phase onsets automatically utilizing the automatic picking
algorithms by Kueperkoch et al. 2010/2012.
@ -152,7 +152,7 @@ def autoPyLoT(inputfile, fnames=None, savepath=None):
print(data)
##########################################################
# !automated picking starts here!
picks = autopickevent(wfdat, parameter)
picks, mainFig = autopickevent(wfdat, parameter, iplot=iplot)
##########################################################
# locating
if locflag == 1:
@ -192,14 +192,14 @@ def autoPyLoT(inputfile, fnames=None, savepath=None):
moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'),
parameter.get('Qp'),
parameter.get('rho'), True, \
parameter.get('iplot'))
iplot)
# update pick with moment property values (w0, fc, Mo)
for station, props in moment_mag.moment_props.items():
picks[station]['P'].update(props)
evt = moment_mag.updated_event()
local_mag = RichterMagnitude(corr_dat, evt,
parameter.get('sstop'), True,\
parameter.get('iplot'))
iplot)
for station, amplitude in local_mag.amplitudes.items():
picks[station]['S']['Ao'] = amplitude.generic_amplitude
evt = local_mag.updated_event()
@ -219,7 +219,7 @@ def autoPyLoT(inputfile, fnames=None, savepath=None):
print("autoPyLoT: Number of maximum iterations reached, stop iterative picking!")
break
print("autoPyLoT: Starting with iteration No. %d ..." % nlloccounter)
picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
picks, _ = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
# write phases to NLLoc-phase file
nll.export(picks, phasefile, parameter)
# remove actual NLLoc-location file to keep only the last
@ -244,14 +244,14 @@ def autoPyLoT(inputfile, fnames=None, savepath=None):
moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'),
parameter.get('Qp'),
parameter.get('rho'), True, \
parameter.get('iplot'))
iplot)
# update pick with moment property values (w0, fc, Mo)
for station, props in moment_mag.moment_props.items():
picks[station]['P'].update(props)
evt = moment_mag.updated_event()
local_mag = RichterMagnitude(corr_dat, evt,
parameter.get('sstop'), True, \
parameter.get('iplot'))
iplot)
for station, amplitude in local_mag.amplitudes.items():
picks[station]['S']['Ao'] = amplitude.generic_amplitude
evt = local_mag.updated_event()
@ -303,6 +303,7 @@ def autoPyLoT(inputfile, fnames=None, savepath=None):
The Python picking and Location Tool\n
************************************'''.format(version=_getVersionString())
print(endsp)
return picks, mainFig
if __name__ == "__main__":
@ -318,6 +319,8 @@ if __name__ == "__main__":
parser.add_argument('-f', '-F', '--fnames', type=str,
action='store',
help='''optional, list of data file names''')
# parser.add_argument('-p', '-P', '--plot', action='store',
# help='show interactive plots')
parser.add_argument('-s', '-S', '--spath', type=str,
action='store',
help='''optional, save path for autoPyLoT output''')
@ -327,4 +330,4 @@ if __name__ == "__main__":
cla = parser.parse_args()
autoPyLoT(str(cla.inputfile), str(cla.fnames), str(cla.spath))
picks, mainFig = autoPyLoT(str(cla.inputfile), str(cla.fnames), str(cla.spath))

2
pylot/RELEASE-VERSION
View File

@ -1 +1 @@
2628-dirty
6563-dirty

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

@ -203,20 +203,18 @@ class RichterMagnitude(Magnitude):
# check for plot flag (for debugging only)
if self.plot_flag > 1:
st.plot()
f = plt.figure(2)
plt.plot(th, sqH)
plt.plot(th[iwin], sqH[iwin], 'g')
plt.plot([t0, t0], [0, max(sqH)], 'r', linewidth=2)
plt.title(
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(th, sqH)
ax.plot(th[iwin], sqH[iwin], 'g')
ax.plot([t0, t0], [0, max(sqH)], 'r', linewidth=2)
ax.title(
'Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \
% (st[0].stats.station, wapp))
plt.xlabel('Time [s]')
plt.ylabel('Displacement [mm]')
plt.show()
raw_input()
plt.close(f)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Displacement [mm]')
return wapp
return wapp, fig
def calc(self):
for a in self.arrivals:
@ -234,7 +232,7 @@ class RichterMagnitude(Magnitude):
continue
delta = degrees2kilometers(a.distance)
onset = pick.time
a0 = self.peak_to_peak(wf, onset)
a0, self.p2p_fig = self.peak_to_peak(wf, onset)
amplitude = ope.Amplitude(generic_amplitude=a0 * 1e-3)
amplitude.unit = 'm'
amplitude.category = 'point'
@ -581,9 +579,6 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [m/Hz]')
plt.grid()
plt.show()
raw_input()
plt.close(f1)
return w0, fc
@ -685,7 +680,7 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
"fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0))
if iplot > 1:
plt.figure(iplot)
plt.figure()#iplot)
plt.loglog(f, S, 'k')
plt.loglog([f[0], fc], [w0, w0], 'g')
plt.loglog([fc, fc], [w0 / 100, w0], 'g')
@ -694,7 +689,7 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [m/Hz]')
plt.grid()
plt.figure(iplot + 1)
plt.figure()#iplot + 1)
plt.subplot(311)
plt.plot(f[il:ir], STD, '*')
plt.title('Common Standard Deviations')
@ -707,8 +702,5 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
plt.plot(f[il:ir], stdfc, '*')
plt.title('Standard Deviations of Corner Frequencies')
plt.xlabel('Corner Frequencies [Hz]')
plt.show()
raw_input()
plt.close()
return w0, fc

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

@ -160,7 +160,7 @@ class AutoPickParameter(object):
fid_out = open(fnout, 'w')
lines = []
for key, value in self.iteritems():
lines.append('{key}\t{value}'.format(key=key, value=value))
lines.append('{key}\t{value}\n'.format(key=key, value=value))
fid_out.writelines(lines)

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

@ -21,16 +21,16 @@ from pylot.core.util.utils import getPatternLine, gen_Pool
from pylot.core.io.data import Data
def autopickevent(data, param):
def autopickevent(data, param, iplot=0):
stations = []
all_onsets = {}
fig_dict = {}
input_tuples = []
# get some parameters for quality control from
# parameter input file (usually autoPyLoT.in).
wdttolerance = param.get('wdttolerance')
mdttolerance = param.get('mdttolerance')
iplot = param.get('iplot')
apverbose = param.get('apverbose')
for n in range(len(data)):
station = data[n].stats.station
@ -41,18 +41,27 @@ def autopickevent(data, param):
for station in stations:
topick = data.select(station=station)
#all_onsets[station] = autopickstation(topick, param, verbose=apverbose)
if not iplot:
input_tuples.append((topick, param, apverbose))
if iplot>0:
all_onsets[station], fig_dict[station] = autopickstation(topick, param, verbose=apverbose, iplot=iplot)
if iplot>0:
print('iPlot Flag active: NO MULTIPROCESSING possible.')
return all_onsets, fig_dict # changing structure of autopicking and figure generation MP MP
pool = gen_Pool()
result = pool.map(call_autopickstation, input_tuples)
pool.close()
for pick in result:
for pick, fig_dict in result:
station = pick['station']
pick.pop('station')
all_onsets[station] = pick
return all_onsets, fig_dict # changing structure of autopicking and figure generation MP MP
# quality control
# median check and jackknife on P-onset times
jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, iplot)
@ -62,10 +71,11 @@ def autopickevent(data, param):
def call_autopickstation(input_tuple):
wfstream, pickparam, verbose = input_tuple
return autopickstation(wfstream, pickparam, verbose)
#multiprocessing not possible with interactive plotting
return autopickstation(wfstream, pickparam, verbose, iplot=0)
def autopickstation(wfstream, pickparam, verbose=False):
def autopickstation(wfstream, pickparam, verbose=False, iplot=0):
"""
:param wfstream: `~obspy.core.stream.Stream` containing waveform
:type wfstream: obspy.core.stream.Stream
@ -82,8 +92,9 @@ def autopickstation(wfstream, pickparam, verbose=False):
# read your autoPyLoT.in for details!
# special parameters for P picking
iplot = iplot
algoP = pickparam.get('algoP')
iplot = pickparam.get('iplot')
pstart = pickparam.get('pstart')
pstop = pickparam.get('pstop')
thosmw = pickparam.get('tlta')
@ -161,6 +172,8 @@ def autopickstation(wfstream, pickparam, verbose=False):
Ao = None # Wood-Anderson peak-to-peak amplitude
picker = 'auto' # type of picks
fig_dict = {}
# split components
zdat = wfstream.select(component="Z")
if len(zdat) == 0: # check for other components
@ -223,6 +236,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
# get prelimenary onset time from AIC-HOS-CF using subclass AICPicker
# of class AutoPicking
aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP)
fig_dict['aicFig'] = aicpick.fig
##############################################################
if aicpick.getpick() is not None:
# check signal length to detect spuriously picked noise peaks
@ -236,7 +250,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
'Decreasing minsiglengh from {0} to ' \
'{1}'.format(minsiglength, minsiglength / 2)
if verbose: print(msg)
Pflag = checksignallength(zne, aicpick.getpick(), tsnrz,
Pflag, fig_dict['slength'] = checksignallength(zne, aicpick.getpick(), tsnrz,
minsiglength / 2,
nfacsl, minpercent, iplot)
else:
@ -253,7 +267,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
trH2_filt.taper(max_percentage=0.05, type='hann')
zne += trH1_filt
zne += trH2_filt
Pflag = checksignallength(zne, aicpick.getpick(), tsnrz,
Pflag, fig_dict['slength'] = checksignallength(zne, aicpick.getpick(), tsnrz,
minsiglength,
nfacsl, minpercent, iplot)
@ -315,12 +329,17 @@ def autopickstation(wfstream, pickparam, verbose=False):
algoP=algoP)
refPpick = PragPicker(cf2, tsnrz, pickwinP, iplot, ausP, tsmoothP,
aicpick.getpick())
fig_dict['refPpick'] = refPpick.fig
mpickP = refPpick.getpick()
#############################################################
if mpickP is not None:
# quality assessment
# get earliest/latest possible pick and symmetrized uncertainty
[epickP, lpickP, Perror] = earllatepicker(z_copy, nfacP, tsnrz,
if iplot:
epickP, lpickP, Perror, fig_dict['el_Ppick'] = earllatepicker(z_copy, nfacP, tsnrz,
mpickP, iplot)
else:
epickP, lpickP, Perror = earllatepicker(z_copy, nfacP, tsnrz,
mpickP, iplot)
# get SNR
@ -342,6 +361,9 @@ def autopickstation(wfstream, pickparam, verbose=False):
# get first motion of P onset
# certain quality required
if Pweight <= minfmweight and SNRP >= minFMSNR:
if iplot:
FM, fig_dict['fm_picker'] = fmpicker(zdat, z_copy, fmpickwin, mpickP, iplot)
else:
FM = fmpicker(zdat, z_copy, fmpickwin, mpickP, iplot)
else:
FM = 'N'
@ -498,18 +520,29 @@ def autopickstation(wfstream, pickparam, verbose=False):
# get refined onset time from CF2 using class Picker
refSpick = PragPicker(arhcf2, tsnrh, pickwinS, iplot, ausS,
tsmoothS, aicarhpick.getpick())
fig_dict['refSpick'] = refSpick.fig
mpickS = refSpick.getpick()
#############################################################
if mpickS is not None:
# quality assessment
# get earliest/latest possible pick and symmetrized uncertainty
h_copy[0].data = trH1_filt.data
[epickS1, lpickS1, Serror1] = earllatepicker(h_copy, nfacS,
if iplot:
epickS1, lpickS1, Serror1, fig_dict['el_S1pick'] = earllatepicker(h_copy, nfacS,
tsnrh,
mpickS, iplot)
else:
epickS1, lpickS1, Serror1 = earllatepicker(h_copy, nfacS,
tsnrh,
mpickS, iplot)
h_copy[0].data = trH2_filt.data
[epickS2, lpickS2, Serror2] = earllatepicker(h_copy, nfacS,
if iplot:
epickS2, lpickS2, Serror2, fig_dict['el_S2pick'] = earllatepicker(h_copy, nfacS,
tsnrh,
mpickS, iplot)
else:
epickS2, lpickS2, Serror2 = earllatepicker(h_copy, nfacS,
tsnrh,
mpickS, iplot)
if epickS1 is not None and epickS2 is not None:
@ -603,39 +636,36 @@ def autopickstation(wfstream, pickparam, verbose=False):
##############################################################
if iplot > 0:
# plot vertical trace
plt.figure()
plt.subplot(3, 1, 1)
fig = plt.figure()
ax1 = fig.add_subplot(311)
tdata = np.arange(0, zdat[0].stats.npts / tr_filt.stats.sampling_rate,
tr_filt.stats.delta)
# check equal length of arrays, sometimes they are different!?
wfldiff = len(tr_filt.data) - len(tdata)
if wfldiff < 0:
tdata = tdata[0:len(tdata) - abs(wfldiff)]
p1, = plt.plot(tdata, tr_filt.data / max(tr_filt.data), 'k')
ax1.plot(tdata, tr_filt.data / max(tr_filt.data), 'k', label='Data')
if Pweight < 4:
p2, = plt.plot(cf1.getTimeArray(), cf1.getCF() / max(cf1.getCF()),
'b')
ax1.plot(cf1.getTimeArray(), cf1.getCF() / max(cf1.getCF()),
'b', label='CF1')
if aicPflag == 1:
p3, = plt.plot(cf2.getTimeArray(),
cf2.getCF() / max(cf2.getCF()), 'm')
p4, = plt.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1],
'r')
plt.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5],
ax1.plot(cf2.getTimeArray(),
cf2.getCF() / max(cf2.getCF()), 'm', label='CF2')
ax1.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1],
'r', label='Initial P Onset')
ax1.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5],
[1, 1], 'r')
plt.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5],
ax1.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5],
[-1, -1], 'r')
p5, = plt.plot([refPpick.getpick(), refPpick.getpick()],
[-1.3, 1.3], 'r', linewidth=2)
plt.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5],
ax1.plot([refPpick.getpick(), refPpick.getpick()],
[-1.3, 1.3], 'r', linewidth=2, label='Final P Pick')
ax1.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5],
[1.3, 1.3], 'r', linewidth=2)
plt.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5],
ax1.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5],
[-1.3, -1.3], 'r', linewidth=2)
plt.plot([lpickP, lpickP], [-1.1, 1.1], 'r--')
plt.plot([epickP, epickP], [-1.1, 1.1], 'r--')
plt.legend([p1, p2, p3, p4, p5],
['Data', 'CF1', 'CF2', 'Initial P Onset',
'Final P Pick'])
plt.title('%s, %s, P Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f '
ax1.plot([lpickP, lpickP], [-1.1, 1.1], 'r--', label='lpp')
ax1.plot([epickP, epickP], [-1.1, 1.1], 'r--', label='epp')
ax1.set_title('%s, %s, P Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f '
'Polarity: %s' % (tr_filt.stats.station,
tr_filt.stats.channel,
Pweight,
@ -643,22 +673,21 @@ def autopickstation(wfstream, pickparam, verbose=False):
SNRPdB,
FM))
else:
plt.legend([p1, p2], ['Data', 'CF1'])
plt.title('%s, P Weight=%d, SNR=None, '
ax1.set_title('%s, P Weight=%d, SNR=None, '
'SNRdB=None' % (tr_filt.stats.channel, Pweight))
else:
plt.title('%s, %s, P Weight=%d' % (tr_filt.stats.station,
ax1.set_title('%s, %s, P Weight=%d' % (tr_filt.stats.station,
tr_filt.stats.channel,
Pweight))
plt.yticks([])
plt.ylim([-1.5, 1.5])
plt.ylabel('Normalized Counts')
plt.suptitle(tr_filt.stats.starttime)
ax1.legend()
ax1.set_yticks([])
ax1.set_ylim([-1.5, 1.5])
ax1.set_ylabel('Normalized Counts')
fig.suptitle(tr_filt.stats.starttime)
if len(edat[0]) > 1 and len(ndat[0]) > 1 and Sflag == 1:
# plot horizontal traces
plt.subplot(3, 1, 2)
ax2 = fig.add_subplot(312)
th1data = np.arange(0,
trH1_filt.stats.npts /
trH1_filt.stats.sampling_rate,
@ -667,50 +696,47 @@ def autopickstation(wfstream, pickparam, verbose=False):
wfldiff = len(trH1_filt.data) - len(th1data)
if wfldiff < 0:
th1data = th1data[0:len(th1data) - abs(wfldiff)]
p21, = plt.plot(th1data, trH1_filt.data / max(trH1_filt.data), 'k')
ax2.plot(th1data, trH1_filt.data / max(trH1_filt.data), 'k', label='Data')
if Pweight < 4:
p22, = plt.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b')
ax2.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b', label='CF1')
if aicSflag == 1:
p23, = plt.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(arhcf2.getCF()), 'm')
p24, = plt.plot(
ax2.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(arhcf2.getCF()), 'm', label='CF2')
ax2.plot(
[aicarhpick.getpick(), aicarhpick.getpick()],
[-1, 1], 'g')
plt.plot(
[-1, 1], 'g', label='Initial S Onset')
ax2.plot(
[aicarhpick.getpick() - 0.5,
aicarhpick.getpick() + 0.5],
[1, 1], 'g')
plt.plot(
ax2.plot(
[aicarhpick.getpick() - 0.5,
aicarhpick.getpick() + 0.5],
[-1, -1], 'g')
p25, = plt.plot([refSpick.getpick(), refSpick.getpick()],
[-1.3, 1.3], 'g', linewidth=2)
plt.plot(
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)
plt.plot(
ax2.plot(
[refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[-1.3, -1.3], 'g', linewidth=2)
plt.plot([lpickS, lpickS], [-1.1, 1.1], 'g--')
plt.plot([epickS, epickS], [-1.1, 1.1], 'g--')
plt.legend([p21, p22, p23, p24, p25],
['Data', 'CF1', 'CF2', 'Initial S Onset',
'Final S Pick'])
plt.title('%s, S Weight=%d, SNR=%7.2f, SNR[dB]=%7.2f' % (
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:
plt.legend([p21, p22], ['Data', 'CF1'])
plt.title('%s, S Weight=%d, SNR=None, SNRdB=None' % (
ax2.set_title('%s, S Weight=%d, SNR=None, SNRdB=None' % (
trH1_filt.stats.channel, Sweight))
plt.yticks([])
plt.ylim([-1.5, 1.5])
plt.ylabel('Normalized Counts')
plt.suptitle(trH1_filt.stats.starttime)
ax2.legend()
ax2.set_yticks([])
ax2.set_ylim([-1.5, 1.5])
ax2.set_ylabel('Normalized Counts')
fig.suptitle(trH1_filt.stats.starttime)
plt.subplot(3, 1, 3)
ax3 = fig.add_subplot(313)
th2data = np.arange(0,
trH2_filt.stats.npts /
trH2_filt.stats.sampling_rate,
@ -719,47 +745,41 @@ def autopickstation(wfstream, pickparam, verbose=False):
wfldiff = len(trH2_filt.data) - len(th2data)
if wfldiff < 0:
th2data = th2data[0:len(th2data) - abs(wfldiff)]
plt.plot(th2data, trH2_filt.data / max(trH2_filt.data), 'k')
ax3.plot(th2data, trH2_filt.data / max(trH2_filt.data), 'k', label='Data')
if Pweight < 4:
p22, = plt.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b')
p22, = ax3.plot(arhcf1.getTimeArray(),
arhcf1.getCF() / max(arhcf1.getCF()), 'b', label='CF1')
if aicSflag == 1:
p23, = plt.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(arhcf2.getCF()), 'm')
p24, = plt.plot(
ax3.plot(arhcf2.getTimeArray(),
arhcf2.getCF() / max(arhcf2.getCF()), 'm', label='CF2')
ax3.plot(
[aicarhpick.getpick(), aicarhpick.getpick()],
[-1, 1], 'g')
plt.plot(
[-1, 1], 'g', label='Initial S Onset')
ax3.plot(
[aicarhpick.getpick() - 0.5,
aicarhpick.getpick() + 0.5],
[1, 1], 'g')
plt.plot(
ax3.plot(
[aicarhpick.getpick() - 0.5,
aicarhpick.getpick() + 0.5],
[-1, -1], 'g')
p25, = plt.plot([refSpick.getpick(), refSpick.getpick()],
[-1.3, 1.3], 'g', linewidth=2)
plt.plot(
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)
plt.plot(
ax3.plot(
[refSpick.getpick() - 0.5, refSpick.getpick() + 0.5],
[-1.3, -1.3], 'g', linewidth=2)
plt.plot([lpickS, lpickS], [-1.1, 1.1], 'g--')
plt.plot([epickS, epickS], [-1.1, 1.1], 'g--')
plt.legend([p21, p22, p23, p24, p25],
['Data', 'CF1', 'CF2', 'Initial S Onset',
'Final S Pick'])
else:
plt.legend([p21, p22], ['Data', 'CF1'])
plt.yticks([])
plt.ylim([-1.5, 1.5])
plt.xlabel('Time [s] after %s' % tr_filt.stats.starttime)
plt.ylabel('Normalized Counts')
plt.title(trH2_filt.stats.channel)
plt.show()
raw_input()
plt.close()
ax3.plot([lpickS, lpickS], [-1.1, 1.1], 'g--', label='lpp')
ax3.plot([epickS, epickS], [-1.1, 1.1], 'g--', label='epp')
ax3.legend()
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)
fig_dict['mainFig'] = fig
##########################################################################
# calculate "real" onset times
if lpickP is not None and lpickP == mpickP:
@ -806,7 +826,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
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
return picks, fig_dict
def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter):
@ -884,7 +904,7 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter):
print("zfac: %f => %f" % (zfac_old, pickparameter.get('zfac')))
# repick station
newpicks = autopickstation(wf2pick, pickparameter)
newpicks, fig = autopickstation(wf2pick, pickparameter)
# replace old dictionary with new one
picks[badpicks[i][0]] = newpicks
@ -899,4 +919,4 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter):
pickparameter.setParam(noisefactor=noisefactor_old)
pickparameter.setParam(zfac=zfac_old)
return picks
return picks, fig

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

@ -72,7 +72,7 @@ class AutoPicker(object):
self.setaus(aus)
self.setTsmooth(Tsmooth)
self.setpick1(Pick1)
self.calcPick()
self.fig = self.calcPick()
def __str__(self):
return '''\n\t{name} object:\n
@ -152,6 +152,7 @@ class AICPicker(AutoPicker):
self.Pick = None
self.slope = None
self.SNR = None
fig = None
# find NaN's
nn = np.isnan(self.cf)
if len(nn) > 1:
@ -225,18 +226,16 @@ 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:
p = plt.figure(self.iplot)
fig = plt.figure() #self.iplot) ### WHY? MP MP
ax = fig.add_subplot(111)
x = self.Data[0].data
p1, = plt.plot(self.Tcf, x / max(x), 'k')
p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r')
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.yticks([])
plt.title(self.Data[0].stats.station)
plt.show()
raw_input()
plt.close(p)
return
ax.plot(self.Tcf, x / max(x), 'k', legend='(HOS-/AR-) Data')
ax.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', legend='Smoothed AIC-CF')
ax.legend()
ax.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax.set_yticks([])
ax.set_title(self.Data[0].stats.station)
return fig
islope = islope[0][0:imax]
dataslope = self.Data[0].data[islope]
# calculate slope as polynomal fit of order 1
@ -253,42 +252,37 @@ class AICPicker(AutoPicker):
self.slope = None
if self.iplot > 1:
p = plt.figure(self.iplot)
fig = plt.figure()#self.iplot)
ax1 = fig.add_subplot(211)
x = self.Data[0].data
p1, = plt.plot(self.Tcf, x / max(x), 'k')
p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r')
ax1.plot(self.Tcf, x / max(x), 'k', label='(HOS-/AR-) Data')
ax1.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
if self.Pick is not None:
p3, = plt.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b', linewidth=2)
plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick'])
else:
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.yticks([])
plt.title(self.Data[0].stats.station)
ax1.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b', linewidth=2, label='AIC-Pick')
ax1.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax1.set_yticks([])
ax1.set_title(self.Data[0].stats.station)
ax1.legend()
if self.Pick is not None:
plt.figure(self.iplot + 1)
p11, = plt.plot(self.Tcf, x, 'k')
p12, = plt.plot(self.Tcf[inoise], self.Data[0].data[inoise])
p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r')
p14, = plt.plot(self.Tcf[islope], dataslope, 'g--')
p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2)
plt.legend([p11, p12, p13, p14, p15],
['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'],
loc='best')
plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
ax2 = fig.add_subplot(212)
ax2.plot(self.Tcf, x, 'k', label='Data')
ax2.plot(self.Tcf[inoise], self.Data[0].data[inoise], label='Noise Window')
ax2.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r', label='Signal Window')
ax2.plot(self.Tcf[islope], dataslope, 'g--', label='Slope Window')
ax2.plot(self.Tcf[islope], datafit, 'g', linewidth=2, label='Slope')
ax2.set_title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station,
self.SNR, self.slope))
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.ylabel('Counts')
plt.yticks([])
plt.show()
raw_input()
plt.close(p)
ax2.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax2.set_ylabel('Counts')
ax2.set_yticks([])
ax2.legend()
if self.Pick == None:
print('AICPicker: Could not find minimum, picking window too short?')
return fig
class PragPicker(AutoPicker):
'''
@ -380,18 +374,17 @@ class PragPicker(AutoPicker):
pickflag = 0
if self.getiplot() > 1:
p = plt.figure(self.getiplot())
p1, = plt.plot(Tcfpick, cfipick, 'k')
p2, = plt.plot(Tcfpick, cfsmoothipick, 'r')
fig = plt.figure()#self.getiplot())
ax = fig.add_subplot(111)
ax.plot(Tcfpick, cfipick, 'k', label='CF')
ax.plot(Tcfpick, cfsmoothipick, 'r', label='Smoothed CF')
if pickflag > 0:
p3, = plt.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', linewidth=2)
plt.legend([p1, p2, p3], ['CF', 'Smoothed CF', 'Pick'])
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.yticks([])
plt.title(self.Data[0].stats.station)
plt.show()
raw_input()
plt.close(p)
ax.plot([self.Pick, self.Pick], [min(cfipick), max(cfipick)], 'b', 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()
return fig
else:
print('PragPicker: No initial onset time given! Check input!')

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

@ -105,35 +105,33 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False):
PickError = symmetrize_error(diffti_te, diffti_tl)
if iplot > 1:
p = plt.figure(iplot)
p1, = plt.plot(t, x, 'k')
p2, = plt.plot(t[inoise], x[inoise])
p3, = plt.plot(t[isignal], x[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
p5, = plt.plot(t[isignal[zc]], np.zeros(len(zc)), '*g',
markersize=14)
plt.legend([p1, p2, p3, p4, p5],
['Data', 'Noise Window', 'Signal Window', 'Noise Level',
'Zero Crossings'],
loc='best')
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2)
plt.plot([LPick, LPick], [max(x) / 2, -max(x) / 2], '--k')
plt.plot([EPick, EPick], [max(x) / 2, -max(x) / 2], '--k')
plt.plot([Pick1 + PickError, Pick1 + PickError],
fig = plt.figure()#iplot)
ax = fig.add_subplot(111)
ax.plot(t, x, 'k', label='Data')
ax.plot(t[inoise], x[inoise], label='Noise Window')
ax.plot(t[isignal], x[isignal], 'r', label='Signal Window')
ax.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k', label='Noise Level')
ax.plot(t[isignal[zc]], np.zeros(len(zc)), '*g',
markersize=14, label='Zero Crossings')
ax.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
ax.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2, label='mpp')
ax.plot([LPick, LPick], [max(x) / 2, -max(x) / 2], '--k', label='lpp')
ax.plot([EPick, EPick], [max(x) / 2, -max(x) / 2], '--k', label='epp')
ax.plot([Pick1 + PickError, Pick1 + PickError],
[max(x) / 2, -max(x) / 2], 'r--', label='spe')
ax.plot([Pick1 - PickError, Pick1 - PickError],
[max(x) / 2, -max(x) / 2], 'r--')
plt.plot([Pick1 - PickError, Pick1 - PickError],
[max(x) / 2, -max(x) / 2], 'r--')
plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
plt.yticks([])
plt.title(
ax.set_xlabel('Time [s] since %s' % X[0].stats.starttime)
ax.set_yticks([])
ax.set_title(
'Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' %
X[0].stats.station)
plt.show()
raw_input()
plt.close(p)
ax.legend()
return EPick, LPick, PickError
if iplot:
return EPick, LPick, PickError, fig
else:
return Epick, LPick, PickError
def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
@ -281,40 +279,36 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
print ("fmpicker: Found polarity %s" % FM)
if iplot > 1:
plt.figure(iplot)
plt.subplot(2, 1, 1)
plt.plot(t, xraw, 'k')
p1, = plt.plot([Pick, Pick], [max(xraw), -max(xraw)], 'b', linewidth=2)
fig = plt.figure()#iplot)
ax1 = fig.add_subplot(211)
ax1.plot(t, xraw, 'k')
ax1.plot([Pick, Pick], [max(xraw), -max(xraw)], 'b', linewidth=2, label='Pick')
if P1 is not None:
p2, = plt.plot(t[islope1], xraw[islope1])
p3, = plt.plot(zc1, np.zeros(len(zc1)), '*g', markersize=14)
p4, = plt.plot(t[islope1], datafit1, '--g', linewidth=2)
plt.legend([p1, p2, p3, p4],
['Pick', 'Slope Window', 'Zero Crossings', 'Slope'],
loc='best')
plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax = plt.gca()
plt.yticks([])
plt.title('First-Motion Determination, %s, Unfiltered Data' % Xraw[
ax1.plot(t[islope1], xraw[islope1], label='Slope Window')
ax1.plot(zc1, np.zeros(len(zc1)), '*g', markersize=14, label='Zero Crossings')
ax1.plot(t[islope1], datafit1, '--g', linewidth=2)
ax1.legend()
ax1.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax1.set_yticks([])
ax1.set_title('First-Motion Determination, %s, Unfiltered Data' % Xraw[
0].stats.station)
plt.subplot(2, 1, 2)
plt.title('First-Motion Determination, Filtered Data')
plt.plot(t, xfilt, 'k')
p1, = plt.plot([Pick, Pick], [max(xfilt), -max(xfilt)], 'b',
ax2=fig.add_subplot(212)
ax2.set_title('First-Motion Determination, Filtered Data')
ax2.plot(t, xfilt, 'k')
ax2.plot([Pick, Pick], [max(xfilt), -max(xfilt)], 'b',
linewidth=2)
if P2 is not None:
p2, = plt.plot(t[islope2], xfilt[islope2])
p3, = plt.plot(zc2, np.zeros(len(zc2)), '*g', markersize=14)
p4, = plt.plot(t[islope2], datafit2, '--g', linewidth=2)
plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax = plt.gca()
plt.xlabel('Time [s] since %s' % Xraw[0].stats.starttime)
plt.yticks([])
plt.show()
raw_input()
plt.close(iplot)
ax2.plot(t[islope2], xfilt[islope2])
ax2.plot(zc2, np.zeros(len(zc2)), '*g', markersize=14)
ax2.plot(t[islope2], datafit2, '--g', linewidth=2)
ax2.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14)
ax2.set_xlabel('Time [s] since %s' % Xraw[0].stats.starttime)
ax2.set_yticks([])
if iplot:
return FM, fig
else:
return FM
@ -606,7 +600,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
# plot results
if iplot > 1:
plt.figure(iplot)
plt.figure()#iplot)
f1, = plt.plot(Ppicks, SPtimes, 'ro')
if wfitflag == 0:
f2, = plt.plot(Ppicks, wdfit, 'k')
@ -621,9 +615,6 @@ def wadaticheck(pickdic, dttolerance, iplot):
plt.ylabel('S-P Times [s]')
plt.xlabel('P Times [s]')
plt.show()
raw_input()
plt.close(iplot)
return checkedonsets
@ -700,25 +691,21 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
returnflag = 0
if iplot == 2:
plt.figure(iplot)
p1, = plt.plot(t, rms, 'k')
p2, = plt.plot(t[inoise], rms[inoise], 'c')
p3, = plt.plot(t[isignal], rms[isignal], 'r')
p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal) - 1]]],
[minsiglevel, minsiglevel], 'g', linewidth=2)
p5, = plt.plot([pick, pick], [min(rms), max(rms)], 'b', linewidth=2)
plt.legend([p1, p2, p3, p4, p5], ['RMS Data', 'RMS Noise Window',
'RMS Signal Window', 'Minimum Signal Level',
'Onset'], loc='best')
plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
plt.ylabel('Counts')
plt.title('Check for Signal Length, Station %s' % X[0].stats.station)
plt.yticks([])
plt.show()
raw_input()
plt.close(iplot)
fig = plt.figure()#iplot)
ax = fig.add_subplot(111)
ax.plot(t, rms, 'k', label='RMS Data')
ax.plot(t[inoise], rms[inoise], 'c', label='RMS Noise Window')
ax.plot(t[isignal], rms[isignal], 'r', label='RMS Signal Window')
ax.plot([t[isignal[0]], t[isignal[len(isignal) - 1]]],
[minsiglevel, minsiglevel], 'g', linewidth=2, label='Minimum Signal Level')
ax.plot([pick, pick], [min(rms), max(rms)], 'b', linewidth=2, label='Onset')
ax.legend()
ax.set_xlabel('Time [s] since %s' % X[0].stats.starttime)
ax.set_ylabel('Counts')
ax.set_title('Check for Signal Length, Station %s' % X[0].stats.station)
ax.set_yticks([])
return returnflag
return returnflag, fig
def checkPonsets(pickdic, dttolerance, iplot):
@ -808,8 +795,6 @@ def checkPonsets(pickdic, dttolerance, iplot):
plt.legend([p1, p2, p3], ['Skipped P Picks', 'Good P Picks', 'Median'],
loc='best')
plt.title('Check P Onsets')
plt.show()
raw_input()
return checkedonsets
@ -962,22 +947,23 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
edat[0].stats.delta)
tn = np.arange(0, ndat[0].stats.npts / ndat[0].stats.sampling_rate,
ndat[0].stats.delta)
plt.plot(tz, z / max(z), 'k')
plt.plot(tz[isignal], z[isignal] / max(z), 'r')
plt.plot(te, edat[0].data / max(edat[0].data) + 1, 'k')
plt.plot(te[isignal], edat[0].data[isignal] / max(edat[0].data) + 1, 'r')
plt.plot(tn, ndat[0].data / max(ndat[0].data) + 2, 'k')
plt.plot(tn[isignal], ndat[0].data[isignal] / max(ndat[0].data) + 2, 'r')
plt.plot([tz[isignal[0]], tz[isignal[len(isignal) - 1]]],
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(tz, z / max(z), 'k')
ax.plot(tz[isignal], z[isignal] / max(z), 'r')
ax.plot(te, edat[0].data / max(edat[0].data) + 1, 'k')
ax.plot(te[isignal], edat[0].data[isignal] / max(edat[0].data) + 1, 'r')
ax.plot(tn, ndat[0].data / max(ndat[0].data) + 2, 'k')
ax.plot(tn[isignal], ndat[0].data[isignal] / max(ndat[0].data) + 2, 'r')
ax.plot([tz[isignal[0]], tz[isignal[len(isignal) - 1]]],
[minsiglevel / max(z), minsiglevel / max(z)], 'g',
linewidth=2)
plt.xlabel('Time [s] since %s' % zdat[0].stats.starttime)
plt.ylabel('Normalized Counts')
plt.yticks([0, 1, 2], [zdat[0].stats.channel, edat[0].stats.channel,
linewidth=2, label='Minimum Signal Level')
ax.set_xlabel('Time [s] since %s' % zdat[0].stats.starttime)
ax.set_ylabel('Normalized Counts')
ax.set_yticks([0, 1, 2], [zdat[0].stats.channel, edat[0].stats.channel,
ndat[0].stats.channel])
plt.title('CheckZ4S, Station %s' % zdat[0].stats.station)
plt.show()
raw_input()
ax.set_title('CheckZ4S, Station %s' % zdat[0].stats.station)
ax.legend()
return returnflag