Skip to content

Commit

Permalink
Streamlining RF plotting and inversion-export routines
Browse files Browse the repository at this point in the history
  • Loading branch information
@geojunky
geojunky committed Dec 13, 2022
1 parent 8fc4904 commit 1aebe58
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 66 deletions.
114 changes: 70 additions & 44 deletions seismic/receiver_fn/bulk_rf_report.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,7 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
t_channel = list(channel)
t_channel[-1] = 'T'
t_channel = ''.join(t_channel)
transverse_data = station_db[t_channel]

rf_stream = rf.RFStream(channel_data).sort(['back_azimuth'])
if event_mask_dict and full_code in event_mask_dict:
Expand All @@ -395,6 +396,10 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
rf_stream = rf.RFStream(
[tr for tr in rf_stream if tr.stats.event_id in event_mask]).sort(['back_azimuth'])
# end if

###############################################################################
# Apply amplitude filter
###############################################################################
if apply_amplitude_filter:
# Label and filter quality
rf_util.label_rf_quality_simple_amplitude(rf_type, rf_stream)
Expand All @@ -406,25 +411,77 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
# end if
# end if

###############################################################################
# Apply slope-ratio filter
###############################################################################
if(min_slope_ratio>0):
rf_stream = rf.RFStream([tr for tr in rf_stream \
if tr.stats.slope_ratio > min_slope_ratio]).sort(['back_azimuth'])
# end if

###############################################################################
# Apply similarity filter
###############################################################################
if apply_similarity_filter and len(rf_stream) >= 3:
rf_stream = rf_util.filter_crosscorr_coeff(rf_stream)
# end if

if not rf_stream:
continue

# Find matching T-component data
###############################################################################
# Filter rf_stream if needed
###############################################################################
if(hk_hpf_freq and hk_hpf_freq>0):
rf_stream.filter(type='highpass', freq=hk_hpf_freq,
corners=1, zerophase=True)
# end if

rf_stream_raw = rf_stream.copy()
###############################################################################
# Apply reverberation filter if needed
###############################################################################
reverberations_removed = False
# Apply reverberation filter if needed
if(rf_corrections.has_reverberations(rf_stream)):
rf_stream = rf_corrections.apply_reverberation_filter(rf_stream)
reverberations_removed = True
# end if

###############################################################################
# RF amplitudes should not exceed 1.0 and should peak around onset time --
# otherwise, such traces are deemed problematic and excluded while computing
# H-k stacks.
###############################################################################
before = len(rf_stream)
rf_stream = rf_util.filter_invalid_radial_component(rf_stream)
after = len(rf_stream)

if(before > after):
print("{}: {}/{} RF traces with amplitudes > 1.0 or troughs around onset time dropped "
"before computing H-k stack ..".format(full_code,
before - after,
before))
# end if

if not rf_stream:
continue

############################################
# Collate streams to plot
############################################
events = [tr.stats.event_id for tr in rf_stream]
transverse_data = station_db[t_channel]
num_traces = len(rf_stream)

rf_stream_raw = rf.RFStream(
[tr for tr in rf_stream_raw if tr.stats.event_id in events]).sort(['back_azimuth'])
t_stream = rf.RFStream(
[tr for tr in transverse_data if tr.stats.event_id in events]).sort(['back_azimuth'])
assert len(t_stream) == num_traces or not t_stream

############################################
# Plot psd
############################################
fig, ax = plt.subplots()
fig.set_size_inches(paper_size_A4[1], paper_size_A4[0])
fig.suptitle('.'.join([nsl, channel]))
Expand All @@ -434,8 +491,10 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
pdf.savefig(dpi=300, orientation='landscape')
plt.close()

############################################
# Plot pinwheel of primary and transverse components
fig = rf_plot_utils.plot_rf_wheel([rf_stream, t_stream], fontscaling=0.8)
############################################
fig = rf_plot_utils.plot_rf_wheel([rf_stream_raw, t_stream], fontscaling=0.8)
fig.set_size_inches(*paper_size_A4)
plt.tight_layout()
plt.subplots_adjust(hspace=0.15, top=0.95, bottom=0.15)
Expand All @@ -445,18 +504,13 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
pdf.savefig(dpi=300, orientation='portrait')
plt.close()

num_traces = len(rf_stream)
assert len(t_stream) == num_traces or not t_stream

# Filter rf_stream if needed
if(hk_hpf_freq and hk_hpf_freq>0):
rf_stream.filter(type='highpass', freq=hk_hpf_freq,
corners=1, zerophase=True)
# end if

# Plot RF stack of primary component
############################################
# Plot RF waveform stacks
############################################
trace_ht = min(total_trace_height_inches/num_traces, max_trace_height)
fig = rf_plot_utils.plot_rf_stack(rf_stream, trace_height=trace_ht, stack_height=fixed_stack_height_inches,
# Plot raw RF stack of primary component
fig = rf_plot_utils.plot_rf_stack(rf_stream_raw, trace_height=trace_ht,
stack_height=fixed_stack_height_inches,
fig_width=paper_size_A4[0])
fig.suptitle("Channel {}".format(rf_stream[0].stats.channel))
# Customize layout to pack to top of page while preserving RF plots aspect ratios
Expand All @@ -465,17 +519,10 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
pdf.savefig(dpi=300, orientation='portrait')
plt.close()

reverberations_removed = False
# Apply reverberation filter if needed
if(rf_corrections.has_reverberations(rf_stream)):
rf_stream = rf_corrections.apply_reverberation_filter(rf_stream)
reverberations_removed = True
# end if

if(reverberations_removed):
# Plot reverberation-filtered RF stack
trace_ht = min(total_trace_height_inches/num_traces, max_trace_height)
fig = rf_plot_utils.plot_rf_stack(rf_stream, trace_height=trace_ht, stack_height=fixed_stack_height_inches,
fig = rf_plot_utils.plot_rf_stack(rf_stream, trace_height=trace_ht,
stack_height=fixed_stack_height_inches,
fig_width=paper_size_A4[0])
fig.suptitle("Channel {} (Reverberations removed)".format(rf_stream[0].stats.channel))
# Customize layout to pack to top of page while preserving RF plots aspect ratios
Expand All @@ -498,27 +545,6 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
plt.close()
# end if

###############################################################################
# RF amplitudes should not exceed 1.0 and should peak around onset time --
# otherwise, such traces are deemed problematic and discarded before computing
# H-k stacks
###############################################################################
before = len(rf_stream)
rf_stream = rf.RFStream([tr for tr in rf_stream if \
((np.max(tr.data) <= 1.0) and \
(np.max(tr.data) > -np.min(tr.data)))
])
after = len(rf_stream)
if(before > after):
print("""{}: {}/{} RF traces with amplitudes > 1.0 or troughs around onset time dropped
before computing H-k stack ..""".format(full_code,
before - after,
before))
# end if

if not rf_stream:
continue

############################################
# Plot H-k stack using primary RF component
############################################
Expand Down
32 changes: 14 additions & 18 deletions seismic/receiver_fn/rf_inversion_export.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,22 +128,6 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
# Select component
data = data.select(component=component)

# RF amplitudes should not exceed 1.0 and should peak around onset time --
# otherwise, such traces are deemed problematic and discarded
before = len(data)
data = rf.RFStream([tr for tr in data if \
((np.max(tr.data) <= 1.0) and \
(np.max(tr.data) > -np.min(tr.data)))
])
after = len(data)
if(before > after):
print('{}.{}: {}/{} RF traces with amplitudes > 1.0 or troughs around onset time dropped ..'.format(
net,
sta,
before - after,
before))
# end if

if apply_amplitude_filter:
# Label and filter quality
rf_util.label_rf_quality_simple_amplitude('ZRT', data)
Expand Down Expand Up @@ -238,6 +222,18 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
before))
# end if

# RF amplitudes should not exceed 1.0 and should peak around onset time --
# otherwise, such traces are deemed problematic and discarded
before = len(ch_traces)
ch_traces = rf_util.filter_invalid_radial_component(ch_traces)
after = len(ch_traces)
if (before > after):
print('{}.{}: {}/{} RF traces with amplitudes > 1.0 or troughs around onset time dropped ..'.format(
network_code, sta,
before - after,
before))
# end if

if (len(ch_traces) == 0):
print('{}.{}: No traces left to stack. Moving along..'.format(network_code, sta))
continue
Expand All @@ -260,7 +256,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_

stack = None
if(apply_phase_weighting):
print('{}.{}: Applying a phase-weighted stack..'.format(network_code, sta))
print('{}.{}: Computing a phase-weighted stack..'.format(network_code, sta))

stack = phase_weighted_stack(ch_traces, phase_weight=pw_exponent)

Expand All @@ -282,7 +278,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
plt.savefig(fn2)
# end if
else:
print('{}.{}: Applying a linear stack..'.format(network_code, sta))
print('{}.{}: Computing a linear stack..'.format(network_code, sta))
stack = ch_traces.stack()
# end if
trace = stack[0]
Expand Down
8 changes: 4 additions & 4 deletions seismic/receiver_fn/rf_plot_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ def plot_rf_stack(rf_stream, time_window=(-10.0, 30.0), trace_height=0.2, stack_
:param rf_stream: RFStream to plot
:type rf_stream: rf.RFStream
:param time_window: Time window to plot, defaults to (-10.0, 25.0)
:param time_window: Time window to plot, defaults to (-10.0, 30.0)
:type time_window: tuple, optional
:param trace_height: Height of a single trace (reduce to cram RFs closer together), defaults to 0.2
:type trace_height: float, optional
Expand All @@ -102,8 +102,9 @@ def plot_rf_stack(rf_stream, time_window=(-10.0, 30.0), trace_height=0.2, stack_
logging.warning('Removed {} traces from RF plot to make it stackable!'.format(num_stackable))
# end if

fig = stackable_stream.plot_rf(fillcolors=('#000000', '#a0a0a0'), trim=time_window, trace_height=trace_height,
stack_height=stack_height, fname=save_file, show_vlines=True, **kwargs)
fig = stackable_stream.plot_rf(fillcolors=('#000000', '#a0a0a0'), trim=time_window,
trace_height=trace_height, stack_height=stack_height,
fname=save_file, show_vlines=True, **kwargs)
return fig
# end func

Expand Down Expand Up @@ -382,7 +383,6 @@ def plot_rf_wheel(rf_stream, max_time=15.0, deg_per_unit_amplitude=45.0, plt_col
return fig
# end func


def plot_iir_filter_response(filter_band_hz, sampling_rate_hz, corners):
"""Plot one-way bandpass filter response in the frequency domain. If filter is used as zero-phase,
the attenuation will be twice what is computed here.
Expand Down
24 changes: 24 additions & 0 deletions seismic/receiver_fn/rf_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -568,3 +568,27 @@ def filter_crosscorr_coeff(rf_stream, time_window=(-2, 25), threshold_cc=0.70, m
kept_data = rf.RFStream([tr for i, tr in enumerate(rf_stream) if keep_trace_mask[i]])
return kept_data
# end func

def filter_invalid_radial_component(rf_stream):
"""
Filter out invalid radial RFs with amplitudes > 1 or troughs around onset time
:param rf_stream: Stream of RF traces to filter, should be **for a single component of a single station**
:type rf_stream: rf.RFStream
:return: Filtered RF stream
:rtype: rf.RFStream
"""
if(len(rf_stream) == 0): return rf_stream

assert rf_stream[0].stats.channel[-1] in ('R', 'Q'), 'Invalid cahnnel; must be either R or Q. Aborting..'
assert_homogenous_stream(rf_stream, filter_invalid_radial_component.__name__)

rf_stream_out = []
for trc in rf_stream:
if ((np.max(trc.data) <= 1.0) and \
(np.max(trc.data) > -np.min(trc.data))):
rf_stream_out.append(trc)
# end if
# end for

return rf.RFStream(rf_stream_out)
# end func

0 comments on commit 1aebe58

Please sign in to comment.