metrics

Functions in metrics.py:

process_trace(tr, inv):

• calls qcTrace(), which adds tr.stats.metrics
• calls trace_quality_factor() to compute tr.stats.quality_factor, and modifies this based on num_gaps, num_overlaps and percent_availability in tr.stats.metrics
• adds tr.stats.metrics.twin
• calls clean_trace()
• calls qcTrace() again

qcTrace(tr):

• writes tr to MiniSEED file
• calls obspy.signal.quality_control.MSEEDMetadata() on that file, saved as mseedqc
• sets tr.stats.metrics to mseedqc.meta

trace_quality_factor(tr):

• initialize quality factor to 1.0
• detects bad trace and exits if:

• < 100 samples (too short)
• sampling_rate < 20 (wrong sampling rate)
• all samples are 0 (dead trace)
• less than 10 unique amplitudes (bit level noise) (otherwise, improves quality_factor)
• calls _check0andMinus1(tr.data) (sequences of 0 or -1 are bad)
• a flat trace (any amplitude) lasts more than 1 second (by calling _get_islands and _FindMaxLength)

• calls _detectClipping() to check if at least 10 samples hit upper or lower bounds of Trace amplitude, and if so, halves quality_factor
• calls _mad_based_outlier() and adds 1 to quality_factor if no outliers found

ampengfft(tr, outdir):

• detrend tr if not already detrended in tr.stats.history
• add peakamp, peaktime and energy to tr.stats.metrics (initialize this if necessary)
• add snr, signal_level and noise_level to tr.stats.metrics by calling signaltonoise(tr)
• add skewness and kurtosis to tr.stats.metrics by calling scipy.stats.describe
• if tr.stats.spectrum exists then IceWeb.icewebSpectrogram() has been called, and we can:

• call _ssam() which averages spectrogram in 1 Hz bands from 0-1, ..., 15-16 Hz, and returns this as tr.stats.ssam (a dict of 'f' and 'A'). This is meant to be used for continuous waveform data.
• call _band_ratio for 1-6 vs 6-11 Hz (appends a dict to tr.stats.bandratio with 'freqlims', 'RSAM_high', 'RSAM_low', 'RSAM_ratio', with the latter being a log2 ratio)
• call _band_ratio for 0.8-4.0 vs 4.0-16.0 Hz
• call _save_esam(), which rather like _ssam() produces an average spectrum but here interpolated to 0.1 Hz resolution from 0-20 Hz. For each event, a line is generated in ESAMYYYYMMDD.csv for each tr with columns id, time, and 20 amplitude values for matching frequencies from 0.0-0.1, ..., 19.9-20.0. This is meant to be used for event waveform data. NOTE: Could add code from eventStatistics() to compute peak time/index

signaltonoise(tr):

• estimates the signal-to-noise ratio
• tr should be cleaned/detrended/filtered before passing to this routine
• checks tr is at least 1-s long, and computes absolute values
• computes the maximum of each 1-s of data, call this time series M
• computes 95th and 5th percentile of M, call these M95 and M5
• estimate signal-to-noise ratio as M95/M5
• adds snr, signal_level and noise_level to tr.stats.metrics
• NOTE: Why not instead just take the ratio of the amplitudes of the "loudest" second and the "quietest" second?

choose_best_traces(st):

• for each tr in st

• sets priority = tr.stats.quality_factor
• optional flags to eliminate seismic and/or infrasound tr
• optional flag to eliminate uncorrected tr
• priority*=2 for seismic Z channel or infrasound channel (effectively downweighting seismic N and E channels)

• chooses the 8 Trace objects with highest priority (change MAX_TRACES for different than 8)
• returns the corresponding indices with st

select_by_index_list(st, chosen): subsets st using indices in chosen. returns new Stream object

peak_amplitudes(st), returns 3 DataFrames:

• seismic3d: vector PGD, PGV and PGA values with corresponding calib and units (of the velocity trace)
• seismic1d: scalar PGD, PGV and PGA values with corresponding calib and units (of the velocity trace)
• infrasound: PP (peak pressure) and PPF (1-20 Hz filtered peak pressure) with corresponding calib and units

max_3c(st):

• detrend st
• compute peak vector amplitude of st

eventStatistics(st) detrends st and returns a DataFrame (one row for each tr in st) which contains:

• id - tr.id
• peakamp - peak amplitude value
• sample - index of peak amplitude
• time - time of peak amplitude
• energy - energy of the tr NOTE: This was written for MiamiLakes, but has some overlap with ampengfft(). That could be a step towards creating PyMSEC for Miami Lakes project events from the Seisan database.