More flexible selection of vertical and horizontal channels (#505)

* let user choose allowed vertical can horizontal chans in template_gen

* add all_vert option

* add options to pick on all vertical/horizontal channels

* add changelog entries

* consistently add channel 3 to default horizontal channels

* add test with OBS data for all vertical channels-options

* hanning window now only accesible as hann in scipy

* limit scipy version until obspy is updated

* make sure all horizontal chan defaults exclude channel '3'

* adjust edge case with horizontal channels named '3'

Co-authored-by: Calum Chamberlain <calum.chamberlain@vuw.ac.nz>

Author: flixha

CHANGES.md

@@ -30,9 +30,17 @@
 * core.lag_calc._xcorr_interp
  - CC-interpolation replaced with resampling (more robust), old method
    deprecated. Use new method with use_new_resamp_method=True as **kwarg.
-* core.lag_calc:
+* core.lag_calc
+ - Added new option all_vert to transfer P-picks to all channels defined as
+   vertical_chans.
+ - Made usage of all_vert, all_horiz consistent across the lag_calc.
  - Fixed bug where minimum CC defined via min_cc_from_mean_cc_factor was not
    set correctly for negative correlation sums.
+* core.template_gen
+ - Added new option all_vert to transfer P-picks to all channels defined as
+   vertical_chans.
+ - Made handling of horizontal_chans and vertical_chans consistent so that user
+   can freely choose relevant channels.
 * utils.correlate
  - Fast Matched Filter now supported natively for version >= 1.4.0
  - Only full correlation stacks are returned now (e.g. where fewer than than

eqcorrscan/core/lag_calc.py

@@ -223,8 +223,9 @@ def _concatenate_and_correlate(streams, template, cores):
 
 def xcorr_pick_family(family, stream, shift_len=0.2, min_cc=0.4,
                       min_cc_from_mean_cc_factor=None,
+                      all_vert=False, all_horiz=False, vertical_chans=['Z'],
                       horizontal_chans=['E', 'N', '1', '2'],
-                      vertical_chans=['Z'], cores=1, interpolate=False,
+                      cores=1, interpolate=False,
                       plot=False, plotdir=None, export_cc=False, cc_dir=None,
                       **kwargs):
     """
@@ -281,7 +282,9 @@ def xcorr_pick_family(family, stream, shift_len=0.2, min_cc=0.4,
     picked_dict = {}
     delta = family.template.st[0].stats.delta
     detect_streams_dict = _prepare_data(
-        family=family, detect_data=stream, shift_len=shift_len)
+        family=family, detect_data=stream, shift_len=shift_len,
+        all_vert=all_vert, all_horiz=all_horiz, vertical_chans=vertical_chans,
+        horizontal_chans=horizontal_chans)
     detection_ids = list(detect_streams_dict.keys())
     detect_streams = [detect_streams_dict[detection_id]
                       for detection_id in detection_ids]
@@ -398,7 +401,9 @@ def xcorr_pick_family(family, stream, shift_len=0.2, min_cc=0.4,
     return picked_dict
 
 
-def _prepare_data(family, detect_data, shift_len):
+def _prepare_data(family, detect_data, shift_len, all_vert=False,
+                  all_horiz=False, vertical_chans=['Z'],
+                  horizontal_chans=['E', 'N', '1', '2']):
     """
     Prepare data for lag_calc - reduce memory here.
 
@@ -425,7 +430,9 @@ def _prepare_data(family, detect_data, shift_len):
                     "samples".format(length))
     prepick = shift_len + family.template.prepick
     detect_streams_dict = family.extract_streams(
-        stream=detect_data, length=length, prepick=prepick)
+        stream=detect_data, length=length, prepick=prepick,
+        all_vert=all_vert, all_horiz=all_horiz, vertical_chans=vertical_chans,
+        horizontal_chans=horizontal_chans)
     for key, detect_stream in detect_streams_dict.items():
         # Split to remove trailing or leading masks
         for i in range(len(detect_stream) - 1, -1, -1):
@@ -453,6 +460,7 @@ def _prepare_data(family, detect_data, shift_len):
 
 def lag_calc(detections, detect_data, template_names, templates,
              shift_len=0.2, min_cc=0.4, min_cc_from_mean_cc_factor=None,
+             all_vert=False, all_horiz=False,
              horizontal_chans=['E', 'N', '1', '2'],
              vertical_chans=['Z'], cores=1, interpolate=False,
              plot=False, plotdir=None, export_cc=False, cc_dir=None, **kwargs):
@@ -599,6 +607,7 @@ def lag_calc(detections, detect_data, template_names, templates,
             template_dict = xcorr_pick_family(
                 family=family, stream=detect_data, min_cc=min_cc,
                 min_cc_from_mean_cc_factor=min_cc_from_mean_cc_factor,
+                all_vert=all_vert, all_horiz=all_horiz,
                 horizontal_chans=horizontal_chans,
                 vertical_chans=vertical_chans, interpolate=interpolate,
                 cores=cores, shift_len=shift_len, plot=plot, plotdir=plotdir,

eqcorrscan/core/match_filter/detection.py

@@ -350,7 +350,9 @@ def _calculate_event(self, template=None, template_st=None,
         self.event = ev
         return self
 
-    def extract_stream(self, stream, length, prepick):
+    def extract_stream(self, stream, length, prepick, all_vert=False,
+                       all_horiz=False, vertical_chans=['Z'],
+                       horizontal_chans=['E', 'N', '1', '2']):
         """
         Extract a cut stream of a given length around the detection.
 
@@ -376,7 +378,17 @@ def extract_stream(self, stream, length, prepick):
             pick = [
                 p for p in self.event.picks
                 if p.waveform_id.station_code == station and
-                p.waveform_id.channel_code == channel]
+                p.waveform_id.channel_code[0:-1] == channel[0:-1]]
+            # Allow picks to be transferred to other vertical/horizontal chans
+            if all_vert and channel[-1] in vertical_chans:
+                pick = [p for p in pick
+                        if p.waveform_id.channel_code[-1] in vertical_chans]
+            elif all_horiz and channel[-1] in horizontal_chans:
+                pick = [p for p in pick
+                        if p.waveform_id.channel_code[-1] in horizontal_chans]
+            else:
+                pick = [p for p in pick
+                        if p.waveform_id.channel_code == channel]
             if len(pick) == 0:
                 Logger.info("No pick for {0}.{1}".format(station, channel))
                 continue

eqcorrscan/core/match_filter/family.py

@@ -507,8 +507,8 @@ def write(self, filename, format='tar', overwrite=False):
         return
 
     def lag_calc(self, stream, pre_processed, shift_len=0.2, min_cc=0.4,
-                 min_cc_from_mean_cc_factor=None,
-                 horizontal_chans=['E', 'N', '1', '2'], vertical_chans=['Z'],
+                 min_cc_from_mean_cc_factor=None, vertical_chans=['Z'],
+                 horizontal_chans=['E', 'N', '1', '2'],
                  cores=1, interpolate=False, plot=False, plotdir=None,
                  parallel=True, process_cores=None, ignore_length=False,
                  ignore_bad_data=False, export_cc=False, cc_dir=None,
@@ -779,7 +779,9 @@ def _process_streams(self, stream, pre_processed, process_cores=1,
             processed_stream = stream.merge()
         return processed_stream.split()
 
-    def extract_streams(self, stream, length, prepick):
+    def extract_streams(self, stream, length, prepick, all_vert=False,
+                        all_horiz=False, vertical_chans=['Z'],
+                        horizontal_chans=['E', 'N', '1', '2']):
         """
         Generate a dictionary of cut streams around detections.
 
@@ -797,7 +799,9 @@ def extract_streams(self, stream, length, prepick):
         """
         # Splitting and merging to remove trailing and leading masks
         return {d.id: d.extract_stream(
-            stream=stream, length=length, prepick=prepick).split().merge()
+            stream=stream, length=length, prepick=prepick, all_vert=all_vert,
+            all_horiz=all_horiz, vertical_chans=vertical_chans,
+            horizontal_chans=horizontal_chans).split().merge()
             for d in self.detections}
 
 

eqcorrscan/core/match_filter/party.py

@@ -833,8 +833,8 @@ def read(self, filename=None, read_detection_catalog=True,
         return self
 
     def lag_calc(self, stream, pre_processed, shift_len=0.2, min_cc=0.4,
-                 min_cc_from_mean_cc_factor=None,
-                 horizontal_chans=['E', 'N', '1', '2'], vertical_chans=['Z'],
+                 min_cc_from_mean_cc_factor=None, vertical_chans=['Z'],
+                 horizontal_chans=['E', 'N', '1', '2'],
                  cores=1, interpolate=False, plot=False, plotdir=None,
                  parallel=True, process_cores=None, ignore_length=False,
                  ignore_bad_data=False, export_cc=False, cc_dir=None,

eqcorrscan/core/template_gen.py

@@ -52,10 +52,11 @@ def __str__(self):
 
 def template_gen(method, lowcut, highcut, samp_rate, filt_order,
                  length, prepick, swin="all", process_len=86400,
-                 all_horiz=False, delayed=True, plot=False, plotdir=None,
-                 return_event=False, min_snr=None, parallel=False,
-                 num_cores=False, save_progress=False, skip_short_chans=False,
-                 **kwargs):
+                 all_vert=False, all_horiz=False, delayed=True, plot=False,
+                 plotdir=None, return_event=False, min_snr=None,
+                 parallel=False, num_cores=False, save_progress=False,
+                 skip_short_chans=False, vertical_chans=['Z'],
+                 horizontal_chans=['E', 'N', '1', '2'], **kwargs):
     """
     Generate processed and cut waveforms for use as templates.
 
@@ -82,6 +83,10 @@ def template_gen(method, lowcut, highcut, samp_rate, filt_order,
         :func:`eqcorrscan.core.template_gen.template_gen`
     :type process_len: int
     :param process_len: Length of data in seconds to download and process.
+    :type all_vert: bool
+    :param all_vert:
+        To use all channels defined in vertical_chans for P-arrivals even if
+        there is only a pick on one of them.  Defaults to False.
     :type all_horiz: bool
     :param all_horiz:
         To use both horizontal channels even if there is only a pick on one of
@@ -119,18 +124,26 @@ def template_gen(method, lowcut, highcut, samp_rate, filt_order,
         Whether to ignore channels that have insufficient length data or not.
         Useful when the quality of data is not known, e.g. when downloading
         old, possibly triggered data from a datacentre
+    :type vertical_chans: list
+    :param vertical_chans:
+        List of channel endings on which P-picks are accepted.
+    :type horizontal_chans: list
+    :param horizontal_chans:
+        List of channel endings for horizontal channels, on which S-picks are
+        accepted.
 
     :returns: List of :class:`obspy.core.stream.Stream` Templates
     :rtype: list
 
     .. note::
         By convention templates are generated with P-phases on the
-        vertical channel and S-phases on the horizontal channels, normal
-        seismograph naming conventions are assumed, where Z denotes vertical
-        and N, E, R, T, 1 and 2 denote horizontal channels, either oriented
-        or not.  To this end we will **only** use Z channels if they have a
-        P-pick, and will use one or other horizontal channels **only** if
-        there is an S-pick on it.
+        vertical channel [can be multiple, e.g., Z (vertical) and H
+        (hydrophone) for an ocean bottom seismometer] and S-phases on the
+        horizontal channels. By default, normal seismograph naming conventions
+        are assumed, where Z denotes vertical and N, E, 1 and 2 denote
+        horizontal channels, either oriented or not.  To this end we will
+        **only** use vertical channels if they have a P-pick, and will use one
+        or other horizontal channels **only** if there is an S-pick on it.
 
     .. warning::
         If there is no phase_hint included in picks, and swin=all, all channels
@@ -388,8 +401,9 @@ def template_gen(method, lowcut, highcut, samp_rate, filt_order,
             # Cut and extract the templates
             template = _template_gen(
                 event.picks, st, length, swin, prepick=prepick, plot=plot,
-                all_horiz=all_horiz, delayed=delayed, min_snr=min_snr,
-                plotdir=plotdir)
+                all_vert=all_vert, all_horiz=all_horiz, delayed=delayed,
+                min_snr=min_snr, vertical_chans=vertical_chans,
+                horizontal_chans=horizontal_chans, plotdir=plotdir)
             process_lengths.append(len(st[0].data) / samp_rate)
             temp_list.append(template)
             catalog_out += event
@@ -582,9 +596,10 @@ def _rms(array):
     return np.sqrt(np.mean(np.square(array)))
 
 
-def _template_gen(picks, st, length, swin='all', prepick=0.05,
+def _template_gen(picks, st, length, swin='all', prepick=0.05, all_vert=False,
                   all_horiz=False, delayed=True, plot=False, min_snr=None,
-                  plotdir=None):
+                  plotdir=None, vertical_chans=['Z'],
+                  horizontal_chans=['E', 'N', '1', '2']):
     """
     Master function to generate a multiplexed template for a single event.
 
@@ -607,6 +622,10 @@ def _template_gen(picks, st, length, swin='all', prepick=0.05,
     :param prepick:
         Length in seconds to extract before the pick time default is 0.05
         seconds.
+    :type all_vert: bool
+    :param all_vert:
+        To use all channels defined in vertical_chans for P-arrivals even if
+        there is only a pick on one of them.  Defaults to False.
     :type all_horiz: bool
     :param all_horiz:
         To use both horizontal channels even if there is only a pick on one
@@ -630,18 +649,26 @@ def _template_gen(picks, st, length, swin='all', prepick=0.05,
     :param plotdir:
         The path to save plots to. If `plotdir=None` (default) then the figure
         will be shown on screen.
+    :type vertical_chans: list
+    :param vertical_chans:
+        List of channel endings on which P-picks are accepted.
+    :type horizontal_chans: list
+    :param horizontal_chans:
+        List of channel endings for horizontal channels, on which S-picks are
+        accepted.
 
     :returns: Newly cut template.
     :rtype: :class:`obspy.core.stream.Stream`
 
     .. note::
         By convention templates are generated with P-phases on the
-        vertical channel and S-phases on the horizontal channels, normal
-        seismograph naming conventions are assumed, where Z denotes vertical
-        and N, E, R, T, 1 and 2 denote horizontal channels, either oriented
-        or not.  To this end we will **only** use Z channels if they have a
-        P-pick, and will use one or other horizontal channels **only** if
-        there is an S-pick on it.
+        vertical channel [can be multiple, e.g., Z (vertical) and H
+        (hydrophone) for an ocean bottom seismometer] and S-phases on the
+        horizontal channels. By default, normal seismograph naming conventions
+        are assumed, where Z denotes vertical and N, E, 1 and 2 denote
+        horizontal channels, either oriented or not.  To this end we will
+        **only** use vertical channels if they have a P-pick, and will use one
+        or other horizontal channels **only** if there is an S-pick on it.
 
     .. note::
         swin argument: Setting to `P` will return only data for channels
@@ -735,13 +762,18 @@ def _template_gen(picks, st, length, swin='all', prepick=0.05,
                     continue
                 starttime.update({'picks': s_pick})
             elif _swin == 'all':
-                if all_horiz and tr.stats.channel[-1] in ['1', '2', '3',
-                                                          'N', 'E']:
+                if all_vert and tr.stats.channel[-1] in vertical_chans:
+                    # Get all picks on vertical channels
+                    channel_pick = [
+                        pick for pick in station_picks
+                        if pick.waveform_id.channel_code[-1] in
+                        vertical_chans]
+                elif all_horiz and tr.stats.channel[-1] in horizontal_chans:
                     # Get all picks on horizontal channels
                     channel_pick = [
                         pick for pick in station_picks
                         if pick.waveform_id.channel_code[-1] in
-                        ['1', '2', '3', 'N', 'E']]
+                        horizontal_chans]
                 else:
                     channel_pick = [
                         pick for pick in station_picks
@@ -751,16 +783,19 @@ def _template_gen(picks, st, length, swin='all', prepick=0.05,
                 starttime.update({'picks': channel_pick})
             elif _swin == 'P':
                 p_pick = [pick for pick in station_picks
-                          if pick.phase_hint.upper()[0] == 'P' and
-                          pick.waveform_id.channel_code == tr.stats.channel]
+                          if pick.phase_hint.upper()[0] == 'P']
+                if not all_vert:
+                    p_pick = [pick for pick in p_pick
+                              if pick.waveform_id.channel_code ==
+                              tr.stats.channel]
                 if len(p_pick) == 0:
                     Logger.debug(
                         f"No picks with phase_hint P "
                         f"found for {tr.stats.station}.{tr.stats.channel}")
                     continue
                 starttime.update({'picks': p_pick})
             elif _swin == 'S':
-                if tr.stats.channel[-1] in ['Z', 'U']:
+                if tr.stats.channel[-1] in vertical_chans:
                     continue
                 s_pick = [pick for pick in station_picks
                           if pick.phase_hint.upper()[0] == 'S']

eqcorrscan/tests/template_gen_test.py

@@ -458,7 +458,8 @@ def test_swin_P(self):
 
     def test_swin_all_and_all_horiz(self):
         template = _template_gen(self.picks, self.st.copy(), 10, swin='all',
-                                 all_horiz=True)
+                                 all_horiz=True,
+                                 horizontal_chans=['E', 'N', '1', '2', '3'])
         for pick in self.picks:
             if pick.phase_hint == 'S':
                 self.assertGreaterEqual(
@@ -487,6 +488,57 @@ def test_triggered_data(self):
         self.assertEqual(len(templates), 0)
 
 
+class TestEdgeGenObs(unittest.TestCase):
+    @classmethod
+    # Extra test case with OBS data with hydrophone channels (HDH) and T-phases
+    def setUpClass(cls):
+        import eqcorrscan
+        cls.testing_path = os.path.dirname(eqcorrscan.__file__) + '/tests'
+        log = logging.getLogger(template_gen_module.__name__)
+        cls._log_handler = MockLoggingHandler(level='DEBUG')
+        log.addHandler(cls._log_handler)
+        cls.log_messages = cls._log_handler.messages
+        cls.st = read(os.path.join(
+            cls.testing_path, 'test_data', 'WAV', 'TEST_',
+            '2019-08-09-1558-47M.NNSN__038'))
+        # for tr in cls.st:
+        #     tr.stats.channel = tr.stats.channel[0] + tr.stats.channel[-1]
+        # Sfile in New Nordic format
+        event = read_events(os.path.join(
+            cls.testing_path, 'test_data', 'REA', 'TEST_',
+            '09-1558-48R.S201908'))[0]
+        cat = filter_picks(
+            Catalog([event]), stations=['KBS', 'OBIN1', 'OBIN2', 'SPA0',
+                                        'NOR', 'DAG', 'HOPEN', 'HSPB'])
+        cls.picks = cat[0].picks
+
+    def setUp(self):
+        self._log_handler.reset()
+
+    def test_swin_all_and_all_vert_and_all_horiz(self):
+        # Test that the hydrophone channel on an OBS is included in the
+        # creation of the vertical channel (P-arrival) template.
+        template = _template_gen(self.picks, self.st.copy(), 20, swin='all',
+                                 all_horiz=True, all_vert=True,
+                                 vertical_chans=['Z', 'H'])
+        for pick in self.picks:
+            if pick.phase_hint and pick.phase_hint[0] == 'P':
+                self.assertGreaterEqual(
+                    len(template.select(
+                        station=pick.waveform_id.station_code,
+                        channel='??[ZH]')), 1)
+            if pick.phase_hint and pick.phase_hint[0] == 'S':
+                self.assertGreaterEqual(
+                    len(template.select(
+                        station=pick.waveform_id.station_code,
+                        channel='??[NE12]')), 2)
+            if pick.phase_hint and pick.phase_hint[0] == 'T':
+                self.assertGreaterEqual(
+                    len(template.select(
+                        station=pick.waveform_id.station_code,
+                        channel='??[ZH]')), 2)
+
+
 class TestDayLong(unittest.TestCase):
     @classmethod
     def setUpClass(cls):