From dcd0bc40d7b24ba57adc248ab15cc329705b0829 Mon Sep 17 00:00:00 2001 From: Darius Arnold Date: Mon, 13 Aug 2018 22:35:38 +0200 Subject: [PATCH] [remove] old autopickstation code --- pylot/core/pick/autopick.py | 977 ------------------------------------ 1 file changed, 977 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 2e466cf9..c64e4a38 100644 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -1184,983 +1184,6 @@ def autopickstation(wfstream, pickparam, verbose=False, iplot=0, fig_dict=None, return None, station_name -def nautopickstation(wfstream, pickparam, verbose=False, - iplot=0, fig_dict=None, metadata=None, origin=None): - """ - picks a single station - :param wfstream: stream object containing waveform of all traces - :type wfstream: ~obspy.core.stream.Stream - :param pickparam: container of picking parameters from input file, usually pylot.in - :type pickparam: pylot.core.io.inputs.PylotParameter - :param verbose: used to control output to log during picking. True = more information printed - :type verbose: bool - :param iplot: logical variable for plotting: 0=none, 1=partial, 2=all - :type iplot: int, (Boolean or String) - :param fig_dict: dictionary containing Matplotlib figures used for plotting picking results during tuning - :type fig_dict: dict - :param metadata: tuple containing metadata type string and Parser object read from inventory file - :type metadata: tuple (str, ~obspy.io.xseed.parser.Parser) - :param origin: list containing origin objects representing origins for all events - :type origin: list(~obspy.core.event.origin) - :return: dictionary containing P pick, S pick and station name - :rtype: dict - """ - - # declaring pickparam variables (only for convenience) - # read your pylot.in for details! - plt_flag = 0 - - # special parameters for P picking - algoP = pickparam.get('algoP') - pstart = pickparam.get('pstart') - pstop = pickparam.get('pstop') - thosmw = pickparam.get('tlta') - tsnrz = pickparam.get('tsnrz') - hosorder = pickparam.get('hosorder') - bpz1 = pickparam.get('bpz1') - bpz2 = pickparam.get('bpz2') - pickwinP = pickparam.get('pickwinP') - aictsmoothP = pickparam.get('aictsmooth') - tsmoothP = pickparam.get('tsmoothP') - ausP = pickparam.get('ausP') - 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') - minAICPslope = pickparam.get('minAICPslope') - minAICPSNR = pickparam.get('minAICPSNR') - timeerrorsP = pickparam.get('timeerrorsP') - # special parameters for S picking - algoS = pickparam.get('algoS') - sstart = pickparam.get('sstart') - sstop = pickparam.get('sstop') - use_taup = real_Bool(pickparam.get('use_taup')) - taup_model = pickparam.get('taup_model') - bph1 = pickparam.get('bph1') - bph2 = pickparam.get('bph2') - tsnrh = pickparam.get('tsnrh') - pickwinS = pickparam.get('pickwinS') - tpred1h = pickparam.get('tpred1h') - tdet1h = pickparam.get('tdet1h') - tpred2h = pickparam.get('tpred2h') - tdet2h = pickparam.get('tdet2h') - Sarorder = pickparam.get('Sarorder') - aictsmoothS = pickparam.get('aictsmoothS') - tsmoothS = pickparam.get('tsmoothS') - ausS = pickparam.get('ausS') - minAICSslope = pickparam.get('minAICSslope') - minAICSSNR = pickparam.get('minAICSSNR') - Srecalcwin = pickparam.get('Srecalcwin') - nfacS = pickparam.get('nfacS') - timeerrorsS = pickparam.get('timeerrorsS') - # parameters for first-motion determination - minFMSNR = pickparam.get('minFMSNR') - fmpickwin = pickparam.get('fmpickwin') - minfmweight = pickparam.get('minfmweight') - # parameters for checking signal length - minsiglength = pickparam.get('minsiglength') - minpercent = pickparam.get('minpercent') - nfacsl = pickparam.get('noisefactor') - # parameter to check for spuriously picked S onset - zfac = pickparam.get('zfac') - # path to inventory-, dataless- or resp-files - - # initialize output - Pweight = 4 # weight for P onset - Sweight = 4 # weight for S onset - FM = 'N' # first motion (polarity) - SNRP = None # signal-to-noise ratio of P onset - SNRPdB = None # signal-to-noise ratio of P onset [dB] - SNRS = None # signal-to-noise ratio of S onset - SNRSdB = None # signal-to-noise ratio of S onset [dB] - mpickP = None # most likely P onset - lpickP = None # latest possible P onset - epickP = None # earliest possible P onset - mpickS = None # most likely S onset - lpickS = None # latest possible S onset - epickS = None # earliest possible S onset - Perror = None # symmetrized picking error P onset - Serror = None # symmetrized picking error S onset - - aicSflag = 0 - aicPflag = 0 - Pflag = 0 - Sflag = 0 - Pmarker = [] - Ao = None # Wood-Anderson peak-to-peak amplitude - picker = 'auto' # type of picks - - def get_components_from_waveformstream(waveformstream): - """ - Splits waveformstream into multiple components zdat, ndat, edat. For traditional orientation (ZNE) these contain - the vertical, north-south or east-west component. Otherwise they contain components numbered 123 with - orientation diverging from the traditional orientation. - :param waveformstream: Stream containing all three components for one station either by ZNE or 123 channel code - (mixture of both options is handled as well) - :type waveformstream: obspy.core.stream.Stream - :return: Tuple containing (z waveform, n waveform, e waveform) selected by the given channels - :rtype: (obspy.core.stream.Stream, obspy.core.stream.Stream, obspy.core.stream.Stream) - """ - - #TODO: get this order from the pylot preferences - channelorder_default = {'Z': 3, 'N': 1, 'E': 2} - waveform_data = {} - for key in channelorder_default: - waveform_data[key] = waveformstream.select(component=key) # try ZNE first - if len(waveform_data[key]) == 0: - waveform_data[key] = waveformstream.select(component=str(channelorder_default[key])) # use 123 as second option - return waveform_data['Z'], waveform_data['N'], waveform_data['E'] - - - def prepare_wfstream_component(wfstream, detrend_type='demean', filter_type='bandpass', freqmin=None, freqmax=None, zerophase=False, taper_max_percentage=0.05, taper_type='hann'): - """ - Prepare a waveformstream for picking by applying detrending, filtering and tapering. Creates a copy of the - waveform the leave the original unchanged. - :param wfstream: - :type wfstream: - :param detrend_type: - :type detrend_type: - :param filter_type: - :type filter_type: - :param freqmin: - :type freqmin: - :param freqmax: - :type freqmax: - :param zerophase: - :type zerophase: - :param taper_max_percentage: - :type taper_max_percentage: - :param taper_type: - :type taper_type: - :return: Tuple containing the changed waveform stream and the first trace of the stream - :rtype: (obspy.core.stream.Stream, obspy.core.trace.Trace) - """ - wfstream_copy = wfstream.copy() - trace_copy = wfstream[0].copy() - trace_copy.detrend(type=detrend_type) - trace_copy.filter(filter_type, freqmin=freqmin, freqmax=freqmax, zerophase=zerophase) - trace_copy.taper(max_percentage=taper_max_percentage, type=taper_type) - wfstream_copy[0].data = trace_copy.data - return wfstream_copy, trace_copy - - # split components - zdat, ndat, edat = get_components_from_waveformstream(wfstream) - - picks = {} - station = wfstream[0].stats.station - - if not zdat: - print('No z-component found for station {}. STOP'.format(station)) - return picks, station - - if p_params['algoP'] == 'HOS' or p_params['algoP'] == 'ARZ' and zdat is not None: - msg = '##################################################\nautopickstation:' \ - ' Working on P onset of station {station}\nFiltering vertical ' \ - 'trace ...\n{data}'.format(station=wfstream[0].stats.station, data=str(zdat)) - if verbose: print(msg) - z_copy, tr_filt = prepare_wfstream_component(zdat, freqmin=p_params['bpz1'][0], freqmax=p_params['bpz1'][1]) - ############################################################## - # check length of waveform and compare with cut times - - # for global seismology: use tau-p method for estimating travel times (needs source and station coords.) - # if not given: sets Lc to infinity to use full stream - if p_params['use_taup'] is True: - Lc = np.inf - print('autopickstation: use_taup flag active.') - 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() - station_coords = metadata.get_coordinates(station_id, time=wfstream[0].stats.starttime) - if station_coords and origin: - source_origin = origin[0] - model = TauPyModel(p_params['taup_model']) - arrivals = model.get_travel_times_geo( - source_origin.depth, - source_origin.latitude, - source_origin.longitude, - station_coords['latitude'], - station_coords['longitude'] - ) - phases = {'P': [], - 'S': []} - for arr in arrivals: - phases[identifyPhaseID(arr.phase.name)].append(arr) - - # get first P and S onsets from arrivals list - arrP, estFirstP = min([(arr, arr.time) for arr in phases['P']], key=lambda t: t[1]) - arrS, estFirstS = min([(arr, arr.time) for arr in phases['S']], key=lambda t: t[1]) - print('autopick: estimated first arrivals for P: {} s, S:{} s after event' - ' origin time using TauPy'.format(estFirstP, estFirstS)) - - # modifiy pstart and pstop relative to estimated first P arrival (relative to station time axis) - p_params['pstart'] += (source_origin.time + estFirstP) - zdat[0].stats.starttime - p_params['pstop']+= (source_origin.time + estFirstP) - zdat[0].stats.starttime - print('autopick: CF calculation times respectively:' - ' pstart: {} s, pstop: {} s'.format(p_params['pstart'], p_params['pstop'])) - elif not origin: - print('No source origins given!') - - # make sure pstart and pstop are inside zdat[0] - pstart = max(p_params['pstart'], 0) - pstop = min(p_params['pstop'], len(zdat[0])*zdat[0].stats.delta) - - if p_params['use_taup'] is False or origin: - Lc = p_params['pstop'] - p_params['pstart'] - - Lwf = zdat[0].stats.endtime - zdat[0].stats.starttime - if not Lwf > 0: - print('autopickstation: empty trace! Return!') - return picks, station - - Ldiff = Lwf - abs(Lc) - 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) - pstart = 0 - pstop = len(zdat[0].data) * zdat[0].stats.delta - cuttimes = [pstart, pstop] - cf1 = None - if p_params['algoP'] == 'HOS': - # calculate HOS-CF using subclass HOScf of class - # CharacteristicFunction - cf1 = HOScf(z_copy, cuttimes, p_params['tlta'], p_params['hosorder']) # instance of HOScf - elif p_params['algoP'] == 'ARZ': - # calculate ARZ-CF using subclass ARZcf of class - # CharcteristicFunction - cf1 = ARZcf(z_copy, cuttimes, p_params['tpred1z'], p_params['Parorder'], p_params['tdet1z'], - p_params['addnoise']) # instance of ARZcf - ############################################################## - # calculate AIC-HOS-CF using subclass AICcf of class - # CharacteristicFunction - # class needs stream object => build it - assert isinstance(cf1, CharacteristicFunction), 'cf2 is not set ' \ - 'correctly: maybe the algorithm name ({algoP}) is ' \ - 'corrupted'.format(algoP=p_params['algoP']) - tr_aic = tr_filt.copy() - tr_aic.data = cf1.getCF() - z_copy[0].data = tr_aic.data - aiccf = AICcf(z_copy, cuttimes) # instance of AICcf - ############################################################## - # get preliminary onset time from AIC-HOS-CF using subclass AICPicker - # of class AutoPicking - key = 'aicFig' - if fig_dict: - fig = fig_dict[key] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - aicpick = AICPicker(aiccf, p_params['tsnrz'], p_params['pickwinP'], iplot, Tsmooth=p_params['aictsmooth'], - fig=fig, linecolor=linecolor) - # add pstart and pstop to aic plot - if fig: - for ax in fig.axes: - ax.vlines(pstart, ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed', label='P start') - ax.vlines(pstop, ax.get_ylim()[0], ax.get_ylim()[1], color='c', linestyles='dashed', label='P stop') - ax.legend(loc=1) - ############################################################## - if aicpick.getpick() is not None: - # check signal length to detect spuriously picked noise peaks - # use all available components to avoid skipping correct picks - # on vertical traces with weak P coda - z_copy[0].data = tr_filt.data - zne = z_copy - if len(ndat) == 0 or len(edat) == 0: - msg = 'One or more horizontal component(s) missing!\n' \ - 'Signal length only checked on vertical component!\n' \ - 'Decreasing minsiglengh from {0} to {1}' \ - .format(signal_length_params['minsiglength'], signal_length_params['minsiglength'] / 2) - if verbose: print(msg) - key = 'slength' - if fig_dict: - fig = fig_dict[key] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - Pflag = checksignallength(zne, aicpick.getpick(), p_params['tsnrz'], - signal_length_params['minsiglength'] / 2, - signal_length_params['noisefactor'], signal_length_params['minpercent'], iplot, - fig, linecolor) - else: - trH1_filt, _ = prepare_wfstream_component(edat, freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1]) - trH2_filt, _ = prepare_wfstream_component(ndat, freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1]) - zne += trH1_filt - zne += trH2_filt - if fig_dict: - fig = fig_dict['slength'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - Pflag = checksignallength(zne, aicpick.getpick(), p_params['tsnrz'], - signal_length_params['minsiglength'], - signal_length_params['noisefactor'], signal_length_params['minpercent'], iplot, - fig, linecolor) - - if Pflag == 1: - # check for spuriously picked S onset - # both horizontal traces needed - if len(ndat) == 0 or len(edat) == 0: - msg = 'One or more horizontal components missing!\n' \ - 'Skipping control function checkZ4S.' - if verbose: print(msg) - else: - if iplot > 1: - if fig_dict: - fig = fig_dict['checkZ4s'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - Pflag = checkZ4S(zne, aicpick.getpick(), s_params['zfac'], - p_params['tsnrz'][2], iplot, fig, linecolor) - if Pflag == 0: - Pmarker = 'SinsteadP' - Pweight = 9 - else: - Pmarker = 'shortsignallength' - Pweight = 9 - ############################################################## - # go on with processing if AIC onset passes quality control - slope = aicpick.getSlope() - if not slope: - slope = 0 - if slope >= p_params['minAICPslope'] and aicpick.getSNR() >= p_params['minAICPSNR'] and Pflag == 1: - aicPflag = 1 - msg = 'AIC P-pick passes quality control: Slope: {0} counts/s, ' \ - 'SNR: {1}\nGo on with refined picking ...\n' \ - 'autopickstation: re-filtering vertical trace ' \ - '...'.format(aicpick.getSlope(), aicpick.getSNR()) - if verbose: print(msg) - # re-filter waveform with larger bandpass - z_copy, tr_filt = prepare_wfstream_component(zdat, freqmin=p_params['bpz2'][0], freqmax=p_params['bpz2'][1]) - ############################################################# - # re-calculate CF from re-filtered trace in vicinity of initial - # onset - cuttimes2 = [round(max([aicpick.getpick() - p_params['Precalcwin'], 0])), - round(min([len(zdat[0].data) * zdat[0].stats.delta, - aicpick.getpick() + p_params['Precalcwin']]))] - cf2 = None - if p_params['algoP'] == 'HOS': - # calculate HOS-CF using subclass HOScf of class - # CharacteristicFunction - cf2 = HOScf(z_copy, cuttimes2, p_params['tlta'], - p_params['hosorder']) # instance of HOScf - elif p_params['algoP'] == 'ARZ': - # calculate ARZ-CF using subclass ARZcf of class - # CharcteristicFunction - cf2 = ARZcf(z_copy, cuttimes2, tpred2z, Parorder, tdet2z, - addnoise) # instance of ARZcf - ############################################################## - # get refined onset time from CF2 using class Picker - assert isinstance(cf2, CharacteristicFunction), 'cf2 is not set ' \ - 'correctly: maybe the algorithm name ({algoP}) is ' \ - 'corrupted'.format(algoP=p_params['algoP']) - if fig_dict: - fig = fig_dict['refPpick'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - refPpick = PragPicker(cf2, p_params['tsnrz'], p_params['pickwinP'], iplot, p_params['ausP'], - p_params['tsmoothP'], aicpick.getpick(), fig, linecolor) - mpickP = refPpick.getpick() - ############################################################# - if mpickP is not None: - # quality assessment - # get earliest/latest possible pick and symmetrized uncertainty - if iplot: - if fig_dict: - fig = fig_dict['el_Ppick'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - epickP, lpickP, Perror = earllatepicker(z_copy, p_params['nfacP'], p_params['tsnrz'], - mpickP, iplot, fig=fig, - linecolor=linecolor) - - # get SNR - SNRP, SNRPdB, Pnoiselevel = getSNR(z_copy, p_params['tsnrz'], mpickP) - - # weight P-onset using symmetric error - Pweight = get_quality_class(Perror, p_params['timeerrorsP']) - - ############################################################## - # get first motion of P onset - # certain quality required - if Pweight <= first_motion_params['minfmweight'] and SNRP >= first_motion_params['minFMSNR']: - if iplot: - if fig_dict: - fig = fig_dict['fm_picker'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - FM = fmpicker(zdat, z_copy, first_motion_params['fmpickwin'], mpickP, iplot, fig, linecolor) - else: - FM = fmpicker(zdat, z_copy, first_motion_params['fmpickwin'], mpickP, iplot) - else: - FM = 'N' - - msg = "autopickstation: P-weight: {0}, " \ - "SNR: {1}, SNR[dB]: {2}, Polarity: {3}".format(Pweight, SNRP, SNRPdB, FM) - print(msg) - msg = 'autopickstation: Refined P-Pick: {} s | P-Error: {} s'.format(zdat[0].stats.starttime \ - + mpickP, Perror) - print(msg) - Sflag = 1 - - else: - msg = 'Bad initial (AIC) P-pick, skipping this onset!\n' \ - 'AIC-SNR={0}, AIC-Slope={1}counts/s\n' \ - '(min. AIC-SNR={2}, ' \ - 'min. AIC-Slope={3}counts/s)'.format(aicpick.getSNR(), - aicpick.getSlope(), - p_params['minAICPSNR'], - p_params['minAICPslope']) - if verbose: print(msg) - Sflag = 0 - - else: - print('autopickstation: No vertical component data available!, ' - 'Skipping station!') - - if ((len(edat) > 0 and len(ndat) == 0) or (len(ndat) > 0 and len(edat) == 0)) and Pweight < 4: - msg = 'Go on picking S onset ...\n' \ - '##################################################\n' \ - 'Only one horizontal component available!\n' \ - 'ARH prediction requires at least 2 components!\n' \ - 'Copying existing horizontal component ...' - if verbose: print(msg) - - # check which component is missing - if len(edat) == 0: - edat = ndat - else: - ndat = edat - - pickSonset = (edat is not None and ndat is not None and len(edat) > 0 and len( - ndat) > 0 and Pweight < 4) - - if pickSonset: - # determine time window for calculating CF after P onset - cuttimesh = [ - round(max([mpickP + s_params['sstart'], 0])), # MP MP relative time axis - round(min([ - mpickP + s_params['sstop'], - edat[0].stats.endtime - edat[0].stats.starttime, - ndat[0].stats.endtime - ndat[0].stats.starttime - ])) - ] - - if not cuttimesh[1] >= cuttimesh[0]: - print('Cut window for horizontal phases too small! Will not pick S onsets.') - pickSonset = False - - if pickSonset: - msg = 'Go on picking S onset ...\n' \ - '##################################################\n' \ - 'Working on S onset of station {0}\nFiltering horizontal ' \ - 'traces ...'.format(edat[0].stats.station) - if verbose: print(msg) - - if s_params['algoS'] == 'ARH': - # re-create stream object including both horizontal components - hdat = edat.copy() - hdat += ndat - h_copy = hdat.copy() - # filter and taper data - trH1_filt = hdat[0].copy() - trH2_filt = hdat[1].copy() - trH1_filt.detrend(type='demean') - trH2_filt.detrend(type='demean') - trH1_filt.filter('bandpass', freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1], - zerophase=False) - trH2_filt.filter('bandpass', freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1], - zerophase=False) - trH1_filt.taper(max_percentage=0.05, type='hann') - trH2_filt.taper(max_percentage=0.05, type='hann') - h_copy[0].data = trH1_filt.data - h_copy[1].data = trH2_filt.data - elif s_params['algoS'] == 'AR3': - # re-create stream object including all components - hdat = zdat.copy() - hdat += edat - hdat += ndat - h_copy = hdat.copy() - # filter and taper data - trH1_filt = hdat[0].copy() - trH2_filt = hdat[1].copy() - trH3_filt = hdat[2].copy() - trH1_filt.detrend(type='demean') - trH2_filt.detrend(type='demean') - trH3_filt.detrend(type='demean') - trH1_filt.filter('bandpass', freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1], - zerophase=False) - trH2_filt.filter('bandpass', freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1], - zerophase=False) - trH3_filt.filter('bandpass', freqmin=s_params['bph1'][0], freqmax=s_params['bph1'][1], - zerophase=False) - trH1_filt.taper(max_percentage=0.05, type='hann') - trH2_filt.taper(max_percentage=0.05, type='hann') - trH3_filt.taper(max_percentage=0.05, type='hann') - h_copy[0].data = trH1_filt.data - h_copy[1].data = trH2_filt.data - h_copy[2].data = trH3_filt.data - ############################################################## - if s_params['algoS'] == 'ARH': - # calculate ARH-CF using subclass ARHcf of class - # CharcteristicFunction - arhcf1 = ARHcf(h_copy, cuttimesh, s_params['tpred1h'], s_params['Sarorder'], s_params['tdet1h'], - p_params['addnoise']) # instance of ARHcf - elif s_params['algoS'] == 'AR3': - # calculate ARH-CF using subclass AR3cf of class - # CharcteristicFunction - arhcf1 = AR3Ccf(h_copy, cuttimesh, s_params['tpred1h'], s_params['Sarorder'], s_params['tdet1h'], - p_params['addnoise']) # instance of ARHcf - ############################################################## - # calculate AIC-ARH-CF using subclass AICcf of class - # CharacteristicFunction - # class needs stream object => build it - tr_arhaic = trH1_filt.copy() - tr_arhaic.data = arhcf1.getCF() - h_copy[0].data = tr_arhaic.data - # calculate ARH-AIC-CF - haiccf = AICcf(h_copy, cuttimesh) # instance of AICcf - ############################################################## - # get prelimenary onset time from AIC-HOS-CF using subclass AICPicker - # of class AutoPicking - if fig_dict: - fig = fig_dict['aicARHfig'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - aicarhpick = AICPicker(haiccf, s_params['tsnrh'], s_params['pickwinS'], iplot, None, - s_params['aictsmoothS'], fig=fig, linecolor=linecolor) - ############################################################### - # go on with processing if AIC onset passes quality control - slope = aicarhpick.getSlope() - if not slope: - slope = 0 - if (slope >= s_params['minAICSslope'] and - aicarhpick.getSNR() >= s_params['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' \ - 'autopickstation: re-filtering horizontal traces ' \ - '...'.format(aicarhpick.getSlope(), aicarhpick.getSNR()) - if verbose: print(msg) - # re-calculate CF from re-filtered trace in vicinity of initial - # onset - cuttimesh2 = [round(aicarhpick.getpick() - s_params['Srecalcwin']), - round(aicarhpick.getpick() + s_params['Srecalcwin'])] - # re-filter waveform with larger bandpass - h_copy = hdat.copy() - # filter and taper data - if s_params['algoS']== 'ARH': - trH1_filt = hdat[0].copy() - trH2_filt = hdat[1].copy() - trH1_filt.detrend(type='demean') - trH2_filt.detrend(type='demean') - trH1_filt.filter('bandpass', freqmin=s_params['bph2'][0], freqmax=s_params['bph2'][1], - zerophase=False) - trH2_filt.filter('bandpass', freqmin=s_params['bph2'][0], freqmax=s_params['bph2'][1], - zerophase=False) - trH1_filt.taper(max_percentage=0.05, type='hann') - trH2_filt.taper(max_percentage=0.05, type='hann') - h_copy[0].data = trH1_filt.data - h_copy[1].data = trH2_filt.data - ############################################################# - arhcf2 = ARHcf(h_copy, cuttimesh2, s_params['tpred2h'], s_params['Sarorder'], s_params['tdet2h'], - p_params['addnoise']) # instance of ARHcf - elif s_params['algoS'] == 'AR3': - trH1_filt = hdat[0].copy() - trH2_filt = hdat[1].copy() - trH3_filt = hdat[2].copy() - trH1_filt.detrend(type='demean') - trH2_filt.detrend(type='demean') - trH3_filt.detrend(type='demean') - trH1_filt.filter('bandpass', freqmin=s_params['bph2'][0], freqmax=s_params['bph2'][1], - zerophase=False) - trH2_filt.filter('bandpass', freqmin=s_params['bph2'][0], freqmax=s_params['bph2'][1], - zerophase=False) - trH3_filt.filter('bandpass', freqmin=s_params['bph2'][0], freqmax=s_params['bph2'][1], - zerophase=False) - trH1_filt.taper(max_percentage=0.05, type='hann') - trH2_filt.taper(max_percentage=0.05, type='hann') - trH3_filt.taper(max_percentage=0.05, type='hann') - h_copy[0].data = trH1_filt.data - h_copy[1].data = trH2_filt.data - h_copy[2].data = trH3_filt.data - ############################################################# - arhcf2 = AR3Ccf(h_copy, cuttimesh2, s_params['tpred2h'], s_params['Sarorder'], s_params['tdet2h'], - p_params['addnoise']) # instance of ARHcf - - # get refined onset time from CF2 using class Picker - if fig_dict: - fig = fig_dict['refSpick'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - refSpick = PragPicker(arhcf2, s_params['tsnrh'], s_params['pickwinS'], iplot, s_params['ausS'], - s_params['tsmoothS'], aicarhpick.getpick(), fig, linecolor) - 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 - if iplot: - if fig_dict: - fig = fig_dict['el_S1pick'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = 'k' - epickS1, lpickS1, Serror1 = earllatepicker(h_copy, s_params['nfacS'], s_params['tsnrh'], mpickS, - iplot, fig=fig, linecolor=linecolor) - else: - epickS1, lpickS1, Serror1 = earllatepicker(h_copy, s_params['nfacS'], s_params['tsnrh'], mpickS, iplot) - - h_copy[0].data = trH2_filt.data - if iplot: - if fig_dict: - fig = fig_dict['el_S2pick'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - else: - fig = None - linecolor = '' - epickS2, lpickS2, Serror2 = earllatepicker(h_copy, s_params['nfacS'], s_params['tsnrh'], mpickS, - iplot, fig=fig, linecolor=linecolor) - else: - epickS2, lpickS2, Serror2 = earllatepicker(h_copy, s_params['nfacS'], s_params['tsnrh'], mpickS, iplot) - if epickS1 is not None and epickS2 is not None: - if s_params['algoS'] == 'ARH': - # get earliest pick of both earliest possible picks - epick = [epickS1, epickS2] - lpick = [lpickS1, lpickS2] - pickerr = [Serror1, Serror2] - if epickS1 is None and epickS2 is not None: - ipick = 1 - elif epickS1 is not None and epickS2 is None: - ipick = 0 - elif epickS1 is not None and epickS2 is not None: - ipick = np.argmin([epickS1, epickS2]) - elif s_params['algoS'] == 'AR3': - [epickS3, lpickS3, Serror3] = earllatepicker(h_copy, s_params['nfacS'], s_params['tsnrh'], - mpickS, iplot) - # get earliest pick of all three picks - epick = [epickS1, epickS2, epickS3] - lpick = [lpickS1, lpickS2, lpickS3] - pickerr = [Serror1, Serror2, Serror3] - if epickS1 is None and epickS2 is not None \ - and epickS3 is not None: - ipick = np.argmin([epickS2, epickS3]) - elif epickS1 is not None and epickS2 is None \ - and epickS3 is not None: - ipick = np.argmin([epickS2, epickS3]) - elif epickS1 is not None and epickS2 is not None \ - and epickS3 is None: - ipick = np.argmin([epickS1, epickS2]) - elif epickS1 is not None and epickS2 is not None \ - and epickS3 is not None: - ipick = np.argmin([epickS1, epickS2, epickS3]) - - epickS = epick[ipick] - lpickS = lpick[ipick] - Serror = pickerr[ipick] - - msg = 'autopickstation: Refined S-Pick: {} s | S-Error: {} s'.format(hdat[0].stats.starttime \ - + mpickS, Serror) - print(msg) - - # get SNR - [SNRS, SNRSdB, Snoiselevel] = getSNR(h_copy, s_params['tsnrh'], mpickS) - - # weight S-onset using symmetric error - if Serror <= s_params['timeerrorsS'][0]: - Sweight = 0 - elif s_params['timeerrorsS'][0] < Serror <= s_params['timeerrorsS'][1]: - Sweight = 1 - elif Perror > s_params['timeerrorsS'][1] and Serror <= s_params['timeerrorsS'][2]: - Sweight = 2 - elif s_params['timeerrorsS'][2] < Serror <= s_params['timeerrorsS'][3]: - Sweight = 3 - elif Serror > s_params['timeerrorsS'][3]: - Sweight = 4 - - print('autopickstation: S-weight: {0}, SNR: {1}, ' - 'SNR[dB]: {2}\n' - '##################################################' - ''.format(Sweight, SNRS, SNRSdB)) - ################################################################ - # get Wood-Anderson peak-to-peak amplitude - # initialize Data object - # re-create stream object including both horizontal components - hdat = edat.copy() - hdat += ndat - else: - msg = 'Bad initial (AIC) S-pick, skipping this onset!\n' \ - 'AIC-SNR={0}, AIC-Slope={1}counts/s\n' \ - '(min. AIC-SNR={2}, ' \ - 'min. AIC-Slope={3}counts/s)\n' \ - '##################################################' \ - ''.format(aicarhpick.getSNR(), aicarhpick.getSlope(), s_params['minAICSSNR'], s_params['minAICSslope']) - if verbose: print(msg) - - ############################################################ - # get Wood-Anderson peak-to-peak amplitude - # initialize Data object - # 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!') - - ############################################################## - try: - iplot = int(iplot) - except ValueError: - if iplot is True or iplot == 'True': - iplot = 2 - else: - iplot = 0 - - if iplot > 0: - # plot vertical trace - if fig_dict is None or fig_dict == 'None': - fig = plt.figure() - plt_flag = 1 - linecolor = 'k' - else: - fig = fig_dict['mainFig'] - linecolor = fig_dict['plot_style']['linecolor']['rgba_mpl'] - fig._tight = True - 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)] - ax1.plot(tdata, tr_filt.data / max(tr_filt.data), color=linecolor, linewidth=0.7, label='Data') - if Pweight < 4: - ax1.plot(cf1.getTimeArray(), cf1.getCF() / max(cf1.getCF()), - 'b', label='CF1') - if aicPflag == 1: - 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') - ax1.plot([aicpick.getpick() - 0.5, aicpick.getpick() + 0.5], - [-1, -1], 'r') - 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) - ax1.plot([refPpick.getpick() - 0.5, refPpick.getpick() + 0.5], - [-1.3, -1.3], 'r', linewidth=2) - 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: - ax1.set_title('%s, P Weight=%d, SNR=None, ' - 'SNRdB=None' % (tr_filt.stats.channel, Pweight)) - else: - ax1.set_title('%s, %s, P Weight=%d' % (tr_filt.stats.station, - tr_filt.stats.channel, - Pweight)) - ax1.legend(loc=1) - ax1.set_yticks([]) - ax1.set_ylim([-1.5, 1.5]) - ax1.set_ylabel('Normalized Counts') - # fig.suptitle(tr_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) - ########################################################################## - # calculate "real" onset times - if lpickP is not None and lpickP == mpickP: - lpickP += zdat[0].stats.delta - if epickP is not None and epickP == mpickP: - epickP -= zdat[0].stats.delta - if mpickP is not None and epickP is not None and lpickP is not None: - lpickP = zdat[0].stats.starttime + lpickP - epickP = zdat[0].stats.starttime + epickP - mpickP = zdat[0].stats.starttime + mpickP - else: - # dummy values (start of seismic trace) in order to derive - # theoretical onset times for iteratve picking - lpickP = zdat[0].stats.starttime + p_params['timeerrorsP'][3] - epickP = zdat[0].stats.starttime - p_params['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: - # no horizontal components given - picks = dict(P=ppick) - return picks, station - - if lpickS is not None and lpickS == mpickS: - lpickS += hdat.stats.delta - if epickS is not None and epickS == mpickS: - epickS -= hdat.stats.delta - if mpickS is not None and epickS is not None and lpickS is not None: - lpickS = hdat.stats.starttime + lpickS - epickS = hdat.stats.starttime + epickS - mpickS = hdat.stats.starttime + mpickS - else: - # dummy values (start of seismic trace) in order to derive - # theoretical onset times for iteratve picking - lpickS = hdat.stats.starttime + s_params['timeerrorsS'][3] - epickS = hdat.stats.starttime - s_params['timeerrorsS'][3] - mpickS = hdat.stats.starttime - - # 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) - return picks, station - - def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter, fig_dict=None): """ Repicking of bad onsets. Uses theoretical onset times from NLLoc-location file.