Better handling of GPS-clock corrections

extract_event_traces.py was dropping half the traces due to masked
numpy arrays being returned when clock-corrections are enabled in
FederatedASDFDataSet. This fix ensures clock-corrected waveform
data are handled more robustly.

seismic/ASDFdatabase/FederatedASDFDataSet.py

@@ -76,6 +76,15 @@ def unique_coordinates(self):
 
     # end func
 
+    def corrections_enabled(self):
+        """
+
+        @return: whether GPS clock-corrections have been enabled by setting
+                 the environment variable GPS_CLOCK_CORRECTION=1
+        """
+        return self.fds.corrections_enabled
+    # end func
+
     def get_closest_stations(self, lon, lat, nn=1):
         """
 

seismic/extract_event_traces.py

@@ -383,7 +383,8 @@ def close(x1, x2):
 # end func
 
 def extract_data(catalog, inventory, waveform_getter, event_trace_datafile,
-                 tt_model='iasp91', wave='P', resample_hz=10, dry_run=True):
+                 tt_model='iasp91', wave='P', resample_hz=10, pad=10,
+                 dry_run=True):
 
     assert wave in ['P', 'S', 'SW'], 'Only P, S and SW (surface wave) is supported. Aborting..'
 
@@ -433,7 +434,9 @@ def extract_data(catalog, inventory, waveform_getter, event_trace_datafile,
     nproc = comm.Get_size()
     rank = comm.Get_rank()
 
+    #################################################
     # data extraction is parallelized over stations
+    #################################################
     nsl_dict = None
     if(rank==0):
         nsl_dict = []
@@ -486,6 +489,7 @@ def extract_data(catalog, inventory, waveform_getter, event_trace_datafile,
                                           phase=phase_map[wave],
                                           tt_model=tt_model, pbar=None,#pbar,
                                           request_window=request_window[wave],
+                                          pad=pad,
                                           dist_range=distance_range[wave]):
                 # Write traces to output file in append mode so that arbitrarily large file
                 # can be processed. If the file already exists, then existing streams will
@@ -634,11 +638,19 @@ def main(inventory_file, network_list, station_list, waveform_database, event_ca
 
     inventory = None
     asdf_dataset = None
+    pad = 10 # nominal padding for waveforms in seconds
     waveform_db_is_web = is_url(waveform_database) or waveform_database in obspy.clients.fdsn.header.URL_MAPPINGS
     if not waveform_db_is_web:
         assert os.path.exists(waveform_database), "Cannot find waveform database file {}".format(waveform_database)
         asdf_dataset = FederatedASDFDataSet(waveform_database)
         inventory = asdf_dataset.get_inventory()
+
+        #################################################
+        # Check if GPS clock-corrections are being applied
+        # A large padding is used to allow for time-shifts
+        # from clock-correction
+        #################################################
+        if (asdf_dataset.corrections_enabled()): pad = 3600
     else:
         assert inventory_file, 'Must provide inventory file if using a URL or an obspy client as waveform source'
         inventory = read_inventory(inventory_file)
@@ -718,7 +730,9 @@ def closure_get_waveforms(network, station, location, channel, starttime, endtim
         if(wave == 'SW'): curr_catalog = sw_catalog(catalog) # remove proximal events for surface-wave output
 
         extract_data(curr_catalog, inventory, waveform_getter, event_trace_datafile,
-                     tt_model=taup_model, wave=wave, resample_hz=resample_hz, dry_run=dry_run)
+                     tt_model=taup_model, wave=wave, resample_hz=resample_hz,
+                     pad=pad,
+                     dry_run=dry_run)
     # end for
     del asdf_dataset
 

seismic/stream_io.py

@@ -18,6 +18,8 @@
 from seismic.units_utils import KM_PER_DEG
 from rf.rfstream import rfstats, obj2stats
 from rf.util import _get_stations
+from obspy.geodetics import gps2dist_azimuth
+from obspy.geodetics import kilometers2degrees
 
 # pylint: disable=invalid-name
 
@@ -92,6 +94,18 @@ def safe_iter_event_data(events, inventory, get_waveforms, use_rfstats=True, pha
                 continue
         else:
             stats = obj2stats(event, coords)
+
+            dist_range = kwargs.get('dist_range')
+            if(dist_range):
+                dist, baz, _ = gps2dist_azimuth(stats.station_latitude,
+                                                stats.station_longitude,
+                                                stats.event_latitude,
+                                                stats.event_longitude)
+                dist = kilometers2degrees(dist / 1e3)
+                if dist_range and not dist_range[0] <= dist <= dist_range[1]:
+                    continue
+                # end if
+            # end if
         # end if
 
         net, sta, loc, cha = seedid.split('.')
@@ -120,11 +134,31 @@ def safe_iter_event_data(events, inventory, get_waveforms, use_rfstats=True, pha
                  'detected for event %s, station %s.'
                  % (len(stream), event.resource_id, seedid))
             continue
+        # end if
+
         if any(isinstance(tr.data, np.ma.masked_array) for tr in stream):
-            from warnings import warn
-            warn('Gaps or overlaps detected for event %s, station %s.'
-                 % (event.resource_id, seedid))
-            continue
+            def has_masked_values(data_stream):
+                rval = False
+                for tr in data_stream:
+                    if(isinstance(tr.data, np.ma.masked_array)):
+                        if (np.any(tr.data.mask)):
+                            rval = True
+                            break
+                        # end if
+                    # end if
+                # end for
+                return rval
+            # end func
+
+            if(has_masked_values(stream)):
+                from warnings import warn
+                warn('Gaps or overlaps detected for event %s, station %s.'
+                     % (event.resource_id, seedid))
+                continue
+            else:
+                for tr in stream: tr.data = np.array(tr.data)
+            # end if
+        # end for
 
         for tr in stream:
             tr.stats.update(stats)