Skip to content

Commit

Permalink
More Flexible Output Generation
Browse files Browse the repository at this point in the history 
* Outputs from the RF workflow can now be restricted to selected
  stations, discriminable by location codes.
  • Loading branch information
@geojunky
geojunky committed Feb 28, 2023
1 parent aff76f7 commit 68b8ab7
Show file tree
Hide file tree
Showing 5 changed files with 57 additions and 31 deletions.
3 changes: 2 additions & 1 deletion seismic/bulk_station_orientations.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,8 @@ def plot_summary(station_coords_dict, corrections_dict, output_fn):
@click.argument('network', type=str, required=True)
@click.option('--output-basename', type=click.Path(dir_okay=False),
help='Output file basename to store results in JSON format and plots in pdf format')
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) to process.', type=str,
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) or a text file '
'with station names in each row, w/wo location codes.', type=str,
show_default=True)
def main(src_h5_event_file, network, output_basename, station_list):
"""
Expand Down
3 changes: 2 additions & 1 deletion seismic/receiver_fn/bulk_rf_report.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -247,7 +247,8 @@ def _plot_hk_solution_point(axes, k, h, idx, phase_time_dict=None):
@click.argument('output-file', type=click.Path(dir_okay=False), required=True)
@click.option('--network-list', default='*', help='A space-separated list of networks (within quotes) to process.', type=str,
show_default=True)
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) to process.', type=str,
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) or a text file '
'with station names in each row, w/wo location codes.', type=str,
show_default=True)
@click.option('--event-mask-folder', type=click.Path(dir_okay=True, exists=True, file_okay=False),
help='Folder containing event masks to use to filter traces. Such masks are generated '
Expand Down
3 changes: 2 additions & 1 deletion seismic/receiver_fn/generate_rf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,7 +267,8 @@ def event_waveforms_to_rf(input_file, output_file, config, network_list='*', sta
@click.argument('output-file', type=click.Path(dir_okay=False), required=True)
@click.option('--network-list', default='*', help='A space-separated list of networks (within quotes) to process.', type=str,
show_default=True)
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) to process.', type=str,
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) or a text file '
'with station names in each row, w/wo location codes.', type=str,
show_default=True)
@click.option('--config-file', type=click.Path(dir_okay=False), default=None,
show_default=True, help="Run configuration file in JSON format")
Expand Down
31 changes: 16 additions & 15 deletions seismic/receiver_fn/rf_inversion_export.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
ch_traces = rf.RFStream([tr for tr in ch_traces if len(tr) == most_common_len])
after = len(ch_traces)
if after < before:
print('{}.{}: {}/{} traces dropped to make them stackable!'.format(network_code, sta,
print('{}.{}.{}: {}/{} traces dropped to make them stackable!'.format(network_code, sta, loc,
before-after, after))
# end if

Expand All @@ -180,21 +180,21 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
if tr.stats.slope_ratio >= min_slope_ratio])
after = len(ch_traces)

print('{}.{}: {}/{} traces dropped with min-slope-ratio filter..'.format(network_code, sta,
before - after,
print('{}.{}.{}: {}/{} traces dropped with min-slope-ratio filter..'.format(network_code, sta,
loc, before - after,
before))
# end if

if (len(ch_traces) == 0):
print('{}.{}: no traces left to process..'.format(network_code, sta))
print('{}.{}.{}: no traces left to process..'.format(network_code, sta, loc))
continue
# end if

# Apply de-reverberation filter if specified
if(dereverberate):
has_reverberations = rf_corrections.has_reverberations(ch_traces)
if (has_reverberations):
print('{}.{}: removing reverberations..'.format(network_code, sta))
print('{}.{}.{}: removing reverberations..'.format(network_code, sta, loc))
ch_traces = rf_corrections.apply_reverberation_filter(ch_traces)
# end if
# end if
Expand All @@ -206,7 +206,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
if((tr.stats.back_azimuth >= baz_range[0]) and
(tr.stats.back_azimuth <= baz_range[1]))])
after = len(ch_traces)
print('{}.{}: {}/{} traces dropped with baz-range filter..'.format(network_code, sta,
print('{}.{}.{}: {}/{} traces dropped with baz-range filter..'.format(network_code, sta, loc,
before - after,
before))
# end if
Expand All @@ -217,8 +217,8 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
ch_traces = rf_util.filter_crosscorr_coeff(rf.RFStream(ch_traces), time_window=trim_window,
apply_moveout=True)
after = len(ch_traces)
print('{}.{}: {}/{} traces dropped with trace-similarity filter..'.format(network_code, sta,
before - after,
print('{}.{}.{}: {}/{} traces dropped with trace-similarity filter..'.format(network_code, sta,
loc, before - after,
before))
# end if

Expand All @@ -228,14 +228,14 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
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,
print('{}.{}.{}: {}/{} RF traces with amplitudes > 1.0 or troughs around onset time dropped ..'.format(
network_code, sta, loc,
before - after,
before))
# end if

if (len(ch_traces) == 0):
print('{}.{}: No traces left to stack. Moving along..'.format(network_code, sta))
print('{}.{}.{}: No traces left to stack. Moving along..'.format(network_code, sta, loc))
continue
# end if

Expand All @@ -244,7 +244,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
# end if

# report stats for traces included in stack
print('{}.{}: Traces included in stack ({}): '.format(network_code, sta, len(ch_traces)))
print('{}.{}.{}: Traces included in stack ({}): '.format(network_code, sta, loc, len(ch_traces)))
for strc in ch_traces:
print('\t Event id, time, lon, lat, baz: {}, {}, {:6.2f}, {:6.2f}, {:6.2f}'.format(
strc.stats.event_id,
Expand All @@ -256,7 +256,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_

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

stack = phase_weighted_stack(ch_traces, phase_weight=pw_exponent)

Expand All @@ -278,7 +278,7 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
plt.savefig(fn2)
# end if
else:
print('{}.{}: Computing a linear stack..'.format(network_code, sta))
print('{}.{}.{}: Computing a linear stack..'.format(network_code, sta, loc))
stack = ch_traces.stack()
# end if
trace = stack[0]
Expand Down Expand Up @@ -309,7 +309,8 @@ def rf_inversion_export(input_h5_file, output_folder, network_list="*", station_
@click.argument('output-folder', type=click.Path(dir_okay=True), required=True)
@click.option('--network-list', default='*', help='A space-separated list of networks (within quotes) to process.', type=str,
show_default=True)
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) to process.', type=str,
@click.option('--station-list', default='*', help='A space-separated list of stations (within quotes) or a text file '
'with station names in each row, w/wo location codes.', type=str,
show_default=True)
@click.option('--station-weights', type=str, default=None, show_default=True,
help='A comma-separated text file containing network and station in the first two columns, '
Expand Down
48 changes: 35 additions & 13 deletions seismic/receiver_fn/rf_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,17 @@ def split_list(lst, npartitions):
return [lst[i * k + min(i, m):(i + 1) * k + min(i + 1, m)] for i in range(npartitions)]
# end func

def trim_hdf_keys(hdf_key_list, networks_string, stations_string):
def trim_hdf_keys(hdf_key_list:[str], networks_string:str, stations_string:str) -> [str]:
"""
Trims a list of hdf_keys, filtering out unwanted networks and stations.
@param hdf_key_list:
@param networks_string: a space-separated list of networks. '*' includes all.
@param stations_string: a space-separated list of stations or a text file
with station names in each row, w/wo location codes.
'*' includes all.
@return: trimmed list
"""

network_list = []
station_list = []

Expand All @@ -40,10 +50,31 @@ def trim_hdf_keys(hdf_key_list, networks_string, stations_string):
if(stations_string=='*'):
station_list = []
else:
stations = []
if(os.path.isfile(stations_string)):
for iline, line in enumerate(open(stations_string, 'r').readlines()):
if(not line.strip()): continue
else: stations.append(line)
# end for

stations_string = ' '.join(stations)
# end if

station_list = re.findall('\S+', stations_string)
assert len(station_list), 'Invalid station list. Aborting..'
# end if

# sanity checks
for net in network_list:
if net not in [hdf_key.split('.')[0] for hdf_key in hdf_key_list]:
assert 0, 'Network {} not found in input dataset. Aborting..'.format(net)
# end for

for sta in station_list:
if sta.split('.')[0] not in [hdf_key.split('.')[1] for hdf_key in hdf_key_list]:
assert 0, 'Station {} not found in input dataset. Aborting..'.format(sta)
# end for

net_subset = [] # filter networks
if(len(network_list)):
for hdf_key in hdf_key_list:
Expand All @@ -57,25 +88,16 @@ def trim_hdf_keys(hdf_key_list, networks_string, stations_string):

sta_subset = [] # filter stations
if(len(station_list)):
for hdf_key in hdf_key_list:
for hdf_key in net_subset:
net, sta, loc = hdf_key.split('.')

if(sta in station_list): sta_subset.append(hdf_key)
if (sta in station_list): sta_subset.append(hdf_key)
if ('.'.join([sta, loc]) in station_list): sta_subset.append(hdf_key)
# end for
else:
sta_subset = net_subset.copy()
# end if

# sanity checks
for net in network_list:
if net not in [hdf_key.split('.')[0] for hdf_key in hdf_key_list]:
assert 0, 'Network {} not found in input dataset. Aborting..'.format(net)
# end for
for sta in station_list:
if sta not in [hdf_key.split('.')[1] for hdf_key in hdf_key_list]:
assert 0, 'Station {} not found in input dataset. Aborting..'.format(sta)
# end for

return sta_subset
# end func

Expand Down

0 comments on commit 68b8ab7

Please sign in to comment.