[remove] old autopickstation code

This commit is contained in:
Darius Arnold 2018-08-13 22:35:38 +02:00
parent b7d3568498
commit dcd0bc40d7

View File

@ -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.