diff --git a/autoPyLoT.py b/autoPyLoT.py index df94e03e..76e00649 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -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)) diff --git a/pylot/RELEASE-VERSION b/pylot/RELEASE-VERSION index 7733445e..3f8a4f95 100644 --- a/pylot/RELEASE-VERSION +++ b/pylot/RELEASE-VERSION @@ -1 +1 @@ -2628-dirty +6563-dirty diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 50231a51..3e68bed1 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -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 diff --git a/pylot/core/io/inputs.py b/pylot/core/io/inputs.py index a200e056..7c57f5ea 100644 --- a/pylot/core/io/inputs.py +++ b/pylot/core/io/inputs.py @@ -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) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index a062992f..12365c00 100644 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -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) - input_tuples.append((topick, param, 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,9 +250,9 @@ 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, - minsiglength / 2, - nfacsl, minpercent, iplot) + Pflag, fig_dict['slength'] = checksignallength(zne, aicpick.getpick(), tsnrz, + minsiglength / 2, + nfacsl, minpercent, iplot) else: # filter and taper horizontal traces trH1_filt = edat.copy() @@ -253,9 +267,9 @@ 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, - minsiglength, - nfacsl, minpercent, iplot) + Pflag, fig_dict['slength'] = checksignallength(zne, aicpick.getpick(), tsnrz, + minsiglength, + nfacsl, minpercent, iplot) if Pflag == 1: # check for spuriously picked S onset @@ -315,13 +329,18 @@ 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, - mpickP, iplot) + 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 [SNRP, SNRPdB, Pnoiselevel] = getSNR(z_copy, tsnrz, mpickP) @@ -342,7 +361,10 @@ def autopickstation(wfstream, pickparam, verbose=False): # get first motion of P onset # certain quality required if Pweight <= minfmweight and SNRP >= minFMSNR: - FM = fmpicker(zdat, z_copy, fmpickwin, mpickP, iplot) + 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,20 +520,31 @@ 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, - tsnrh, - mpickS, iplot) - + 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, - tsnrh, - mpickS, iplot) + 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: if algoS == 'ARH': # get earliest pick of both earliest possible picks @@ -603,62 +636,58 @@ 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 ' - 'Polarity: %s' % (tr_filt.stats.station, - tr_filt.stats.channel, - Pweight, - SNRP, - SNRPdB, - FM)) + 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, + SNRP, + SNRPdB, + FM)) else: - plt.legend([p1, p2], ['Data', 'CF1']) - plt.title('%s, P Weight=%d, SNR=None, ' - 'SNRdB=None' % (tr_filt.stats.channel, Pweight)) + 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 diff --git a/pylot/core/pick/picker.py b/pylot/core/pick/picker.py index cc6ee7d9..446a245a 100644 --- a/pylot/core/pick/picker.py +++ b/pylot/core/pick/picker.py @@ -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,41 +252,36 @@ 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, - 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 = 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)) + 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): @@ -296,7 +290,7 @@ class PragPicker(AutoPicker): ''' def calcPick(self): - + if self.getpick1() is not None: print('PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...') @@ -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!') diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 3a45fc55..6ead7072 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -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,41 +279,37 @@ 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([]) - return FM + if iplot: + return FM, fig + else: + return FM def crossings_nonzero_all(data): @@ -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]]], - [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, - ndat[0].stats.channel]) - plt.title('CheckZ4S, Station %s' % zdat[0].stats.station) - plt.show() - raw_input() + 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, 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]) + ax.set_title('CheckZ4S, Station %s' % zdat[0].stats.station) + ax.legend() return returnflag