Skip to content

Commit 26ca96f

Browse files
geojunkygeojunky
committed
Incorporating ANT Velocities in RF-CCP Workflow
* Velocities from ANT models can now be incorporated while migrating RFs. * Added options to relax RF quality checks * Other miscellaneous changes
1 parent 0dd9319 commit 26ca96f

File tree

5 files changed

+174
-82
lines changed

5 files changed

+174
-82
lines changed

‎seismic/receiver_fn/bulk_rf_report.py‎

Lines changed: 106 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,6 @@
66
from mpi4py import MPI
77
import os
88
import re
9-
import logging
109
import itertools
1110

1211
import numpy as np
@@ -31,11 +30,9 @@
3130
from shutil import rmtree
3231
from collections import defaultdict
3332
from seismic.misc import split_list
34-
35-
logging.basicConfig()
33+
from seismic.misc import setup_logger
3634

3735
paper_size_A4 = (8.27, 11.69) # inches
38-
3936
DEFAULT_HK_SOLN_LABEL = 'global'
4037

4138
def _get_aspect(ax):
@@ -96,7 +93,7 @@ def _rf_layout_A4(fig):
9693
def _produce_hk_stacking(channel_data, Vp=rf_stacking.DEFAULT_Vp,
9794
weighting=rf_stacking.DEFAULT_WEIGHTS,
9895
semblance_weighted=True, labelling=DEFAULT_HK_SOLN_LABEL,
99-
depth_colour_range=(20, 70)):
96+
depth_colour_range=(20, 70), log=None):
10097
"""Helper function to produce H-k stacking figure."""
10198

10299
k_grid, h_grid, hk_stack = rf_stacking.compute_hk_stack(channel_data,
@@ -130,10 +127,10 @@ def _produce_hk_stacking(channel_data, Vp=rf_stacking.DEFAULT_Vp,
130127
if labelling == 'global':
131128
h_max, k_max = rf_stacking.find_global_hk_maximum(k_grid, h_grid, hk_stack)
132129
soln = [(h_max, k_max)]
133-
log.info("Numerical solution (H, k) = ({:.3f}, {:.3f})".format(*soln[0]))
130+
if(log): log.info("Numerical solution (H, k) = ({:.3f}, {:.3f})".format(*soln[0]))
134131
elif labelling == 'local':
135132
soln = rf_stacking.find_local_hk_maxima(k_grid, h_grid, hk_stack)
136-
log.info("Numerical solutions (H, k) = {}".format(soln))
133+
if(log): log.info("Numerical solutions (H, k) = {}".format(soln))
137134
# end if
138135

139136
# ==========================================================================
@@ -169,7 +166,7 @@ def compute_arrival_times(p, h, k):
169166
# end func
170167

171168
def _produce_sediment_hk_stacking(channel_data, H_c, k_c, Vp=rf_stacking.DEFAULT_Vp,
172-
labelling=DEFAULT_HK_SOLN_LABEL):
169+
labelling=DEFAULT_HK_SOLN_LABEL, log=None):
173170
"""Helper function to produce H-k stacking figure."""
174171

175172
k_grid, h_grid, hk_stack, weighting = rf_stacking.compute_sediment_hk_stack(channel_data,
@@ -200,10 +197,10 @@ def _produce_sediment_hk_stacking(channel_data, H_c, k_c, Vp=rf_stacking.DEFAULT
200197
if labelling == 'global':
201198
h_max, k_max = rf_stacking.find_global_hk_maximum(k_grid, h_grid, hk_stack)
202199
soln = [(h_max, k_max)]
203-
log.info("Numerical solution (H, k) = ({:.3f}, {:.3f})".format(*soln[0]))
200+
if(log): log.info("Numerical solution (H, k) = ({:.3f}, {:.3f})".format(*soln[0]))
204201
elif labelling == 'local':
205202
soln = rf_stacking.find_local_hk_maxima(k_grid, h_grid, hk_stack)
206-
log.info("Numerical solutions (H, k) = {}".format(soln))
203+
if(log): log.info("Numerical solutions (H, k) = {}".format(soln))
207204
# end if
208205

209206
# Plot the local maxima
@@ -254,8 +251,11 @@ def _plot_hk_solution_point(axes, k, h, idx, phase_time_dict=None):
254251
@click.option('--event-mask-folder', type=click.Path(dir_okay=True, exists=True, file_okay=False),
255252
help='Folder containing event masks to use to filter traces. Such masks are generated '
256253
'using rf_handpick_tool')
257-
@click.option('--apply-amplitude-filter', is_flag=True, default=False, show_default=True,
258-
help='Apply RF amplitude filtering to the RFs. The default filtering logic includes: '
254+
@click.option('--relax-sanity-checks', is_flag=True, default=False, show_default=True,
255+
help='RF traces with amplitudes > 1.0 or troughs around onset time are dropped by default. '
256+
'This option allows RF traces with amplitudes > 1.0 to pass through')
257+
@click.option('--apply-simple-quality-filter', is_flag=True, default=False, show_default=True,
258+
help='Apply RF snr and amplitude filtering to the RFs. The default filtering logic includes: '
259259
'Signal SNR >= 2.0 '
260260
'RMS amplitude of signal < 0.2 '
261261
'Maximum amplitude of signal < 1.0')
@@ -291,8 +291,8 @@ def _plot_hk_solution_point(axes, k, h, idx, phase_time_dict=None):
291291
@click.option('--disable-semblance-weighting', is_flag=True, default=False, show_default=True,
292292
help='Disables default semblance-weighting applied to H-k stacks')
293293
def main(input_file, output_file, network_list='*', station_list='*', event_mask_folder='',
294-
apply_amplitude_filter=False, apply_similarity_filter=False, min_slope_ratio=-1,
295-
force_dereverberate=None, filter_by_distance=None,
294+
relax_sanity_checks=False, apply_simple_quality_filter=False, apply_similarity_filter=False,
295+
min_slope_ratio=-1, force_dereverberate=None, filter_by_distance=None,
296296
hk_vp=rf_stacking.DEFAULT_Vp, hk_weights=rf_stacking.DEFAULT_WEIGHTS,
297297
hk_solution_labels=DEFAULT_HK_SOLN_LABEL, depth_colour_range=(20, 70),
298298
hk_hpf_freq=None, hk_hpf_after_dereverberation=False, disable_semblance_weighting=False):
@@ -303,8 +303,13 @@ def main(input_file, output_file, network_list='*', station_list='*', event_mask
303303
OUTPUT_FILE : Output pdf file name
304304
"""
305305

306+
log = setup_logger('__func__', log_file=os.path.splitext(output_file)[0] + '.log')
306307

307-
log.setLevel(logging.INFO)
308+
# check parameter compatibility
309+
if(relax_sanity_checks and apply_simple_quality_filter):
310+
raise RuntimeError('Options --relax-sanity-checks and --apply-simple-quality-filter '
311+
'are mutually exclusive. Aborting..')
312+
# end if
308313

309314
comm = MPI.COMM_WORLD
310315
nproc = comm.Get_size()
@@ -416,8 +421,6 @@ def match_stations(sta, sta_list):
416421
total_trace_height_inches = paper_size_A4[1] - fixed_stack_height_inches - y_pad_inches
417422
max_trace_height = 0.2
418423

419-
log.setLevel(logging.WARNING)
420-
421424
curr_output_file = os.path.join(tempdir, '{}.pdf'.format(nsl))
422425

423426
with PdfPages(curr_output_file) as pdf:
@@ -443,35 +446,91 @@ def match_stations(sta, sta_list):
443446
transverse_data = station_db[t_channel]
444447

445448
rf_stream = rf.RFStream(channel_data).sort(['back_azimuth'])
449+
450+
451+
###############################################################################
452+
# Remove streams from masked events
453+
###############################################################################
446454
if event_mask_dict and full_code in event_mask_dict:
447455
# Select events from external source
448456
event_mask = event_mask_dict[full_code]
457+
before = len(rf_stream)
449458
rf_stream = rf.RFStream(
450459
[tr for tr in rf_stream if tr.stats.event_id in event_mask]).sort(['back_azimuth'])
460+
after = len(rf_stream)
461+
462+
if (before > after):
463+
log.info("{}: {}/{} RF traces from masked events dropped..".
464+
format(nsl,
465+
before - after,
466+
before))
467+
# end if
468+
if (after == 0):
469+
log.warning("All RF traces for {} were from masked events..".format(nsl))
470+
# end if
471+
# end if
472+
473+
if not rf_stream: continue
474+
475+
###############################################################################
476+
# RF amplitudes should not exceed 1.0 and should peak around onset time --
477+
# otherwise, such traces are deemed problematic and excluded while computing
478+
# H-k stacks. Only applies for P RFs.
479+
###############################################################################
480+
if(rf_type == 'prf'):
481+
before = len(rf_stream)
482+
rf_stream = rf_util.\
483+
filter_invalid_radial_component(rf_stream,
484+
allow_rfs_over_unity=relax_sanity_checks)
485+
after = len(rf_stream)
486+
487+
if(before > after):
488+
log.info("{}: {}/{} RF traces with amplitudes > 1.0 or troughs around "
489+
"onset time dropped before computing H-k stack ..".
490+
format(nsl,
491+
before - after,
492+
before))
493+
# end if
494+
if (after == 0):
495+
log.warning("Sanity check filter has removed all traces for {}..".format(nsl))
496+
# end if
451497
# end if
452498

499+
if not rf_stream: continue
500+
453501
###############################################################################
454502
# Restrict RFs by distance if specified
455503
###############################################################################
456504
if st in filter_by_distance_dict.keys():
457505
min_dist, max_dist = filter_by_distance_dict[st]
458506

507+
before = len(rf_stream)
459508
rf_stream = rf_util.filter_by_distance(rf_stream, min_dist, max_dist)
460-
if (len(rf_stream) == 0):
509+
after = len(rf_stream)
510+
log.info('{}: ({}/{}) traces dropped using distance filter..'.format(nsl,
511+
before-after,
512+
before))
513+
if(after == 0):
461514
log.warning("Filter-by-distance has removed all traces for {}.".format(nsl))
462515
# end if
463516
# end if
464517

465518
###############################################################################
466519
# Apply amplitude filter
467520
###############################################################################
468-
if apply_amplitude_filter:
521+
if apply_simple_quality_filter:
469522
# Label and filter quality
470523
rf_util.label_rf_quality_simple_amplitude(rf_rot, rf_stream)
524+
before = len(rf_stream)
471525
rf_stream = rf.RFStream(
472526
[tr for tr in rf_stream if tr.stats.predicted_quality == 'a']).sort(['back_azimuth'])
473-
if(len(rf_stream) == 0):
474-
log.warning("Amplitude filter has removed all traces for {}. "
527+
after = len(rf_stream)
528+
log.info('{}: ({}/{}) traces dropped using predicted quality filter..'.
529+
format(nsl,
530+
before-after,
531+
before))
532+
if(after == 0):
533+
log.warning("Predicted quality filter has removed all traces for {}. "
475534
"Ensure rf_quality_filter was run beforehand..".format(nsl))
476535
# end if
477536
# end if
@@ -480,15 +539,35 @@ def match_stations(sta, sta_list):
480539
# Apply slope-ratio filter
481540
###############################################################################
482541
if(min_slope_ratio>0):
542+
before = len(rf_stream)
483543
rf_stream = rf.RFStream([tr for tr in rf_stream \
484544
if tr.stats.slope_ratio > min_slope_ratio]).sort(['back_azimuth'])
545+
after = len(rf_stream)
546+
log.info('{}: ({}/{}) traces dropped using min-slope-ratio filter..'.
547+
format(nsl,
548+
before - after,
549+
before))
550+
if (after == 0):
551+
log.warning("Min-slope-ratio filter has removed all traces for {}..".format(nsl))
552+
# end if
485553
# end if
486554

555+
if not rf_stream: continue
556+
487557
###############################################################################
488558
# Apply similarity filter
489559
###############################################################################
490560
if apply_similarity_filter and len(rf_stream) >= 3:
561+
before = len(rf_stream)
491562
rf_stream = rf_util.filter_crosscorr_coeff(rf_stream)
563+
after = len(rf_stream)
564+
log.info('{}: ({}/{}) traces dropped using similarity filter..'.
565+
format(nsl,
566+
before - after,
567+
before))
568+
if (after == 0):
569+
log.warning("Similarity filter has removed all traces for {}..".format(nsl))
570+
# end if
492571
# end if
493572

494573
if not rf_stream: continue
@@ -509,31 +588,11 @@ def match_stations(sta, sta_list):
509588
reverberations_removed = False
510589
# Apply reverberation filter if needed (or forced). Only applies for P RFs
511590
if((rf_corrections.has_reverberations(rf_stream) or st in fd_sta_list) and
512-
(rf_type == 'prf')):
591+
(rf_type == 'prf')):
513592
rf_stream = rf_corrections.apply_reverberation_filter(rf_stream)
514593
reverberations_removed = True
515594
# end if
516595

517-
###############################################################################
518-
# RF amplitudes should not exceed 1.0 and should peak around onset time --
519-
# otherwise, such traces are deemed problematic and excluded while computing
520-
# H-k stacks. Only applies for P RFs.
521-
###############################################################################
522-
if(rf_type == 'prf'):
523-
before = len(rf_stream)
524-
rf_stream = rf_util.filter_invalid_radial_component(rf_stream)
525-
after = len(rf_stream)
526-
527-
if(before > after):
528-
print("{}: {}/{} RF traces with amplitudes > 1.0 or troughs around "
529-
"onset time dropped before computing H-k stack ..".format(full_code,
530-
before - after,
531-
before))
532-
# end if
533-
# end if
534-
535-
if not rf_stream: continue
536-
537596
###############################################################################
538597
# Filter rf_stream if needed, after dereverberation, as requested
539598
###############################################################################
@@ -632,7 +691,8 @@ def match_stations(sta, sta_list):
632691
fig, maxima = _produce_hk_stacking(rf_stream, Vp=hk_vp, weighting=hk_weights,
633692
semblance_weighted=(not disable_semblance_weighting),
634693
labelling=hk_solution_labels,
635-
depth_colour_range=depth_colour_range)
694+
depth_colour_range=depth_colour_range,
695+
log=log)
636696
hk_soln[nsl] = maxima
637697
station_coords[nsl] = (channel_data[0].stats.station_longitude, channel_data[0].stats.station_latitude)
638698

@@ -647,7 +707,7 @@ def match_stations(sta, sta_list):
647707
if(len(hk_soln[nsl])):
648708
H_c = hk_soln[nsl][0][0]
649709
k_c = hk_soln[nsl][0][1]
650-
fig, maxima = _produce_sediment_hk_stacking(rf_stream, H_c=H_c, k_c=k_c, Vp=hk_vp)
710+
fig, maxima = _produce_sediment_hk_stacking(rf_stream, H_c=H_c, k_c=k_c, Vp=hk_vp, log=log)
651711
sediment_hk_soln[nsl] = maxima
652712

653713
sediment_station_coords[nsl] = (channel_data[0].stats.station_longitude,
@@ -732,12 +792,11 @@ def flatten_dict_list(dict_list):
732792
if(rank == 0):
733793
rmtree(tempdir)
734794

735-
print("Finishing...")
736-
print("bulk_rf_report SUCCESS!")
795+
log.info("Finishing...")
796+
log.info("bulk_rf_report SUCCESS!")
737797
# end if
738798
# end main
739799

740800
if __name__ == "__main__":
741-
log = logging.getLogger(__name__)
742801
main() # pylint: disable=no-value-for-parameter
743802
# end if

‎seismic/receiver_fn/plot_ccp.py‎

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
import click
2020
import json
2121
import matplotlib.pyplot as plt
22+
import pandas as pd
2223
from PIL.PngImagePlugin import PngImageFile, PngInfo
2324

2425
from seismic.receiver_fn.rf_ccp_util import CCPVolume, CCP_VerticalProfile, CCP_DepthProfile, Gravity
@@ -95,6 +96,10 @@ def validate_start_end(line):
9596
@click.command(name='vertical', context_settings=CONTEXT_SETTINGS)
9697
@click.argument('ccp-h5-volume', type=click.Path(exists=True, dir_okay=False), required=True)
9798
@click.argument('profile-def', type=click.Path(exists=True, dir_okay=False), required=True)
99+
@click.option('--hk-estimates', type=click.Path(exists=True, dir_okay=False), default=None, show_default=True,
100+
help='H-K stacking results in csv format, with the following columns: '
101+
'Net.Sta.loc, Lon, Lat, H0, k0, H1, k1, H2, k2. Depth estimates are plotted '
102+
'on CCP profiles, if available for a given station.')
98103
@click.option('--gravity-grid', type=click.Path(exists=True, dir_okay=False), default=None, show_default=True,
99104
help='Gravity grid in tif/ers format. If provided, a gravity line-plot is added to each'
100105
' vertical profile. Note that this parameter is ignored if the output-format is set '
@@ -145,9 +150,9 @@ def validate_start_end(line):
145150
help='Output format')
146151
@click.option('--plot-aspect', type=click.Choice(['auto', 'equal']), default='auto', show_default=True,
147152
help='Aspect ratio of CCP plot. Default is "auto", which applies a vertical exaggeration.')
148-
def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx, dz, max_depth, swath_width, ds,
149-
extend, cell_radius, idw_exponent, pw_exponent, amplitude_min, amplitude_max, max_station_dist,
150-
output_folder, output_format, plot_aspect):
153+
def vertical(ccp_h5_volume, profile_def, hk_estimates, gravity_grid, mt_sgrid, mt_utm_zone, dx, dz, max_depth,
154+
swath_width, ds, extend, cell_radius, idw_exponent, pw_exponent, amplitude_min, amplitude_max,
155+
max_station_dist, output_folder, output_format, plot_aspect):
151156
""" Plot CCP vertical profile \n
152157
CCP_H5_VOLUME: Path to CCP volume in H5 format (output of rf_3dmigrate.py) or a text file containing paths to CCP volumes in H5 format. The latter is useful for generating profiles spanning multiple neighboring networks \n
153158
PROFILE_DEF: text file containing start and end locations of each vertical profile as\n
@@ -171,6 +176,12 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx
171176
log.info('Loading profile definition file {}..'.format(profile_def))
172177
profiles = read_profile_def(profile_def, vol, 'vertical')
173178

179+
hk = None
180+
if(hk_estimates):
181+
log.info('Loading hk estimates {}..'.format(hk_estimates))
182+
hk = pd.read_csv(hk_estimates, delimiter=',', skipinitialspace=True)
183+
# end if
184+
174185
gravity = None
175186
if(gravity_grid):
176187
log.info('Loading gravity grid {}..'.format(gravity_grid))
@@ -237,7 +248,7 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx
237248
ax.set_aspect(plot_aspect)
238249

239250
# plot profile
240-
vprof.plot(ax, amp_min=amplitude_min, amp_max=amplitude_max, gax=gax, gravity=gravity)
251+
vprof.plot(ax, amp_min=amplitude_min, amp_max=amplitude_max, gax=gax, hk=hk, gravity=gravity)
241252

242253
fig.savefig(fname)
243254

0 commit comments

Comments
 (0)