Extending CCP-profile Plotting Functionality

* CCP profiles can now be plotted based on multiple migrated volumes
* Minor changes to help and error messages in plot_ccp.py

seismic/receiver_fn/plot_ccp.py

@@ -46,7 +46,10 @@ def validate_start_end(line):
             start, end = [item.strip() for item in line.split('-')]
             for item in [start, end]:
                 if (item not in ccpVolume._meta.keys()):
-                    raise ValueError('{} not found in {}. Aborting..'.format(item, ccpVolume._fn))
+                    raise ValueError('{} not found in {}. '
+                                     'Available stations are: {}. '
+                                     'Aborting..'.format(item, ccpVolume._fn_list,
+                                                         list(ccpVolume._meta.keys())))
                 # end if
             # end for
         except Exception as e:
@@ -146,7 +149,7 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx
              extend, cell_radius, idw_exponent, pw_exponent, amplitude_min, amplitude_max, max_station_dist,
              output_folder, output_format, plot_aspect):
     """ Plot CCP vertical profile \n
-    CCP_H5_VOLUME: Path to CCP volume in H5 format (output of rf_3dmigrate.py)\n
+    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
     PROFILE_DEF: text file containing start and end locations of each vertical profile as\n
                  NET.STA.LOC-NET.STA.LOC in each line\n\n
 
@@ -155,6 +158,9 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx
     """
     log.setLevel(logging.DEBUG)
 
+    # create output folder if necessary
+    if not os.path.exists(output_folder): os.makedirs(output_folder)
+
     # sanity checks
     if (mt_sgrid and not mt_utm_zone): assert 0, 'UTM zone for {} not specified. ' \
                                                 'See help for --mt-utm-zone'.format(mt_sgrid)
@@ -295,7 +301,7 @@ def vertical(ccp_h5_volume, profile_def, gravity_grid, mt_sgrid, mt_utm_zone, dx
 def depth(ccp_h5_volume, profile_def, mt_sgrid, mt_utm_zone, dx, dy, dz, extend, cell_radius,
           idw_exponent, pw_exponent, amplitude_min, amplitude_max, output_folder, output_format):
     """ Plot CCP depth profile \n
-    CCP_H5_VOLUME: Path to CCP volume in H5 format (output of rf_3dmigrate.py)\n
+    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
     PROFILE_DEF: text file containing start and end locations, defining the rectangular region spanned by each
                  depth profile, in the first line as NET.STA.LOC-NET.STA.LOC, followed by depth(s) (km) in
                  subsequent lines \n\n
@@ -305,6 +311,9 @@ def depth(ccp_h5_volume, profile_def, mt_sgrid, mt_utm_zone, dx, dy, dz, extend,
     """
     log.setLevel(logging.DEBUG)
 
+    # create output folder if necessary
+    if not os.path.exists(output_folder): os.makedirs(output_folder)
+
     # sanity checks
     if (mt_sgrid and not mt_utm_zone): assert 0, 'UTM zone for {} not specified. ' \
                                                 'See help for --mt-utm-zone'.format(mt_sgrid)

seismic/receiver_fn/rf_ccp_util.py

@@ -349,44 +349,63 @@ def process_streams(self, output_file, fmin=None, fmax=None, primary_component='
 
 class CCPVolume():
     def __init__(self, fn):
-        self._fn = fn
+        def get_volume_filenames(input_string):
+            vfns = []
+            if(not h5py.is_hdf5(input_string)):
+                for iline, line in enumerate(open(input_string, 'r').readlines()):
+                    line = line.strip()
+                    if(len(line)): vfns.append(line)
+                # end for
+            else:
+                vfns.append(input_string)
+            # end if
+
+            for fn in vfns:
+                assert h5py.is_hdf5(fn), 'Invalid H5 file: {}. Aborting..'.format(fn)
+            # end for
+
+            return vfns
+        # end func
+
+        self._fn_list = get_volume_filenames(fn)
         self._meta = defaultdict(list)
         self._earth_radius = None
-        self._data = None
+        self._data = []
         self._tree = None
 
-        hf = None
-        try:
-            hf = h5py.File(self._fn, 'r')
-        except Exception as e:
-            print(str(e))
+        for fn in self._fn_list:
+            hf = None
+            try:
+                hf = h5py.File(fn, 'r')
+            except Exception as e:
+                print(str(e))
 
-            assert 0, 'Failed to load {}. Aborting..'.format(self._fn)
-        # end try
+                assert 0, 'Failed to load {}. Aborting..'.format(fn)
+            # end try
 
-        # read metadata
-        for k in hf.attrs.keys():
-            if (k == 'earth_radius'):
-                self._earth_radius = hf.attrs[k]
-        # end for
+            # read metadata
+            for k in hf.attrs.keys():
+                if (k == 'earth_radius'):
+                    self._earth_radius = hf.attrs[k]
+            # end for
 
-        # read station metadata
-        for dk in hf.keys():
-            for sk in hf[dk].attrs.keys():
-                self._meta[sk] = hf[dk].attrs[sk] # station -> lon, lat, elevation, hasReverberations
+            # read station metadata
+            for dk in hf.keys():
+                for sk in hf[dk].attrs.keys():
+                    self._meta[sk] = hf[dk].attrs[sk] # station -> lon, lat, elevation, hasReverberations
+                # end for
             # end for
-        # end for
 
-        # read ccp volume
-        self._data = []
-        for k in hf.keys():
-            self._data.append(np.array(hf[k]))
+            # read ccp volume
+            for k in hf.keys():
+                self._data.append(np.array(hf[k]))
+            # end for
+
+            hf.close()
         # end for
-        self._data = np.vstack(self._data)
 
+        self._data = np.vstack(self._data)
         self._tree = cKDTree(self._data[:, :3], balanced_tree=False)
-
-        hf.close()
     # end func
 # end class