Making RF Plotting Routines More Flexible

* Exposed more parameters from the migration routines to the CLI
* Minor fix in bulk_rf_report.py

seismic/receiver_fn/bulk_rf_report.py

@@ -30,6 +30,7 @@
 import uuid
 from shutil import rmtree
 from collections import defaultdict
+from seismic.misc import split_list
 
 logging.basicConfig()
 
@@ -352,7 +353,7 @@ def match_stations(sta, sta_list):
         proc_hdfkeys = rf_util.trim_hdf_keys(proc_hdfkeys, network_list, station_list)
 
         # split work-load over all procs
-        proc_hdfkeys = rf_util.split_list(proc_hdfkeys, nproc)
+        proc_hdfkeys = split_list(proc_hdfkeys, nproc)
 
         tempdir = os.path.join(os.path.dirname(output_file), str(uuid.uuid4()))
         os.makedirs(tempdir, exist_ok=True)

seismic/receiver_fn/plot_ccp.py

@@ -109,7 +109,7 @@ def validate_start_end(line):
               help='Horizontal distance-step (km) between start and end locations of a profile')
 @click.option('--dz', type=float, default=0.5, show_default=True,
               help='Depth-step (km)')
-@click.option('--max-depth', type=click.FloatRange(0, 150), default=100, show_default=True,
+@click.option('--max-depth', type=click.FloatRange(0, 750), default=100, show_default=True,
               help='Maximum depth (km) of profile')
 @click.option('--swath-width', type=float, default=40, show_default=True,
               help='CCP amplitudes are averaged over grid elements across a swath defined by '

seismic/receiver_fn/rf_3dmigrate.py

@@ -18,7 +18,7 @@
 import sys
 import logging
 import click
-from seismic.receiver_fn.rf_ccp_util import Migrate
+from seismic.receiver_fn.rf_ccp_util import Migrator
 
 logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)
 log = logging.getLogger('rf_3dmigrate')
@@ -27,18 +27,22 @@
 @click.command(context_settings=CONTEXT_SETTINGS)
 @click.argument('rf-h5-file', type=click.Path(exists=True, dir_okay=False), required=True)
 @click.argument('output_h5_file', type=click.Path(exists=False, dir_okay=False), required=True)
-@click.option('--min-slope-ratio', type=float, default=-1, show_default=True,
-              help='Apply filtering to the RFs based on the "slope_ratio" metric '
-                   'that indicates robustness of P-arrival. Typically, a minimum '
-                   'slope-ratio of 5 is able to pick out strong arrivals. The '
-                   'default value of -1 does not apply this filter')
+@click.option('--dz', type=float, default=0.1, show_default=True,
+              help='Depth-step (km)')
+@click.option('--max-depth', type=click.FloatRange(0, 750), default=150, show_default=True,
+              help='Maximum depth (km) of profile')
 @click.option('--fmin', type=float, default=None, show_default=True,
               help="Lowest frequency for bandpass filter; default is None."
                    "If only --fmin in provided, a highpass filter is aplied.")
 @click.option('--fmax', type=float, default=None, show_default=True,
               help="Highest frequency for bandpass filter; default is None."
                    "If only --fmax is provided, a lowpass filter is applied.")
-def main(rf_h5_file, output_h5_file, fmin, fmax, min_slope_ratio):
+@click.option('--min-slope-ratio', type=float, default=-1, show_default=True,
+              help='Apply filtering to the RFs based on the "slope_ratio" metric '
+                   'that indicates robustness of P-arrival. Typically, a minimum '
+                   'slope-ratio of 5 is able to pick out strong arrivals. The '
+                   'default value of -1 does not apply this filter')
+def main(rf_h5_file, output_h5_file, dz, max_depth, fmin, fmax, min_slope_ratio):
     """Perform 3D migration of RFs
     RF_H5_FILE : Path to RFs in H5 format
     OUTPUT_H5_FILE: H5 output file name
@@ -48,7 +52,8 @@ def main(rf_h5_file, output_h5_file, fmin, fmax, min_slope_ratio):
     """
     log.setLevel(logging.DEBUG)
 
-    m = Migrate(rf_filename=rf_h5_file, min_slope_ratio=min_slope_ratio, logger=log)
+    m = Migrator(rf_filename=rf_h5_file, dz=dz, max_depth=max_depth,
+                 min_slope_ratio=min_slope_ratio, logger=log)
     m.process_streams(output_h5_file, fmin=fmin, fmax=fmax)
 # end
 

seismic/receiver_fn/rf_ccp_util.py

@@ -34,7 +34,7 @@
 from affine import Affine
 import struct
 from seismic.stream_io import get_obspyh5_index
-from seismic.misc import rtp2xyz
+from seismic.misc import rtp2xyz, split_list
 
 class Gravity:
     def __init__(self, gravity_grid_fn):
@@ -87,9 +87,11 @@ def query(self, lon, lat):
     # end func
 # end class
 
-class Migrate:
-    def __init__(self, rf_filename, min_slope_ratio=-1, earth_radius=6371, logger=None):
+class Migrator:
+    def __init__(self, rf_filename, dz, max_depth, min_slope_ratio=-1, earth_radius=6371, logger=None):
         self._rf_filename = rf_filename
+        self._dz = dz
+        self._max_depth = max_depth
         self._min_slope_ratio = min_slope_ratio
         self._earth_radius = earth_radius #km
         self._logger = logger
@@ -110,7 +112,7 @@ def process_streams(self, output_file, fmin=None, fmax=None, primary_component='
             assert len(hkeys), 'No hdf5-groups found in file {}. Aborting..'.format(self._rf_filename)
 
             # split workload
-            proc_hkeys = rf_util.split_list(hkeys, self._nproc)
+            proc_hkeys = split_list(hkeys, self._nproc)
         # end if
         # broadcast workload to all procs
         proc_hkeys = self._comm.bcast(proc_hkeys, root=0)
@@ -231,8 +233,8 @@ def process_streams(self, output_file, fmin=None, fmax=None, primary_component='
             #===============================================================
             tck = splrep(times, depths, k=3, s=0)
 
-            NZ = 1500
-            ZMAX = 150 # km
+            ZMAX = self._max_depth # km
+            NZ = np.int_(np.ceil(ZMAX/self._dz)) + 1
             dnew = np.linspace(1e-5, ZMAX, NZ)
             tnew = np.zeros(dnew.shape)
             for idx in np.arange(dnew.shape[0]):
@@ -645,10 +647,10 @@ def plot(self, ax, amp_min=-0.2, amp_max=0.2, gax=None, gravity=None):
         # end if
 
         tickstep_x = 50
-        tickstep_y = 5
+        tickstep_y = 10
 
-        ax.set_xlabel('Distance [km]', fontsize=12)
-        ax.set_ylabel('Depth [km]', fontsize=12)
+        ax.set_xlabel('Distance [km]', fontsize=10)
+        ax.set_ylabel('Depth [km]', fontsize=10)
         ax.tick_params(direction="in", labelleft=True, labelright=True)
 
         ax.set_xticks(np.arange(0, np.max(self._gx) * 1.0001, tickstep_x))