Rf quality filter (#241)

* Minor bug fixes * Minor improvement to rf plots
GeoscienceAustralia · Apr 8, 2022 · ae5db59 · ae5db59
1 parent 45e0ef7
commit ae5db59
Show file tree

Hide file tree

Showing 4 changed files with 17 additions and 8 deletions.
diff --git a/seismic/extract_event_traces.py b/seismic/extract_event_traces.py
@@ -31,7 +31,7 @@
 from obspy.taup import TauPyModel
 
 from PhasePApy.phasepapy.phasepicker import aicdpicker
-from seismic.pick_harvester.utils import Event, Origin
+from seismic.pick_harvester.utils import Event, Origin, Magnitude
 from seismic.pick_harvester.pick import extract_p, extract_s
 from seismic.stream_processing import zerophase_resample
 
@@ -273,6 +273,8 @@ def trim_inventory(inventory, network_list, station_list):
 #end func
 
 class Picker():
+    counter = 0 # static variable to count number of calls made
+
     def __init__(self, taup_model_name):
         self._taup_model = TauPyModel(model=taup_model_name)
         self._picker_list_p = []
@@ -299,6 +301,8 @@ def pick(self, ztrace, ntrace, etrace, phase='P'):
                         ztrace.stats.event_depth)
         event = Event()
         event.preferred_origin = origin
+        event.public_id = Picker.counter; Picker.counter += 1
+        event.preferred_magnitude = Magnitude(mag=ztrace.stats.event_magnitude, mag_type='M')
 
         result = None
         if(phase == 'P'):

diff --git a/seismic/pick_harvester/pick.py b/seismic/pick_harvester/pick.py
@@ -64,8 +64,8 @@ def extract_p(taupy_model, pickerlist, event, station_longitude, station_latitud
     try:
         snrst = st.slice(po.utctime + tat + win_start + buffer_start, po.utctime + tat + win_end + buffer_end)
         snrst = snrst.copy()
-        zerophase_resample(snrst, resample_hz)
         snrst.detrend('linear')
+        zerophase_resample(snrst, resample_hz)
     except:
         return None
     # end try
@@ -178,14 +178,14 @@ def extract_s(taupy_model, pickerlist, event, station_longitude, station_latitud
     try:
         stn = stn.slice(po.utctime + tat + win_start + buffer_start, po.utctime + tat + win_end + buffer_end)
         stn = stn.copy()
-        zerophase_resample(stn, resample_hz)
         stn.detrend('linear')
+        zerophase_resample(stn, resample_hz)
 
         if (ste):
             ste = ste.slice(po.utctime + tat + win_start + buffer_start, po.utctime + tat + win_end + buffer_end)
             ste = ste.copy()
-            zerophase_resample(ste, resample_hz)
             ste.detrend('linear')
+            zerophase_resample(ste, resample_hz)
         # end if
 
         if (ste):

diff --git a/seismic/receiver_fn/generate_rf.py b/seismic/receiver_fn/generate_rf.py
@@ -176,9 +176,8 @@ def event_waveforms_to_rf(input_file, output_file, config, network_list='*', sta
                     logger.info('Rank {}: {}: Applying a baz correction '
                                 'of {}'.format(rank,
                                                nsl,
-                                               result['.'.join([net, sta])]['azimuth_correction']))
-                    ned.apply(lambda st: correct_back_azimuth(None, st,
-                                                              result['.'.join([net, sta])]['azimuth_correction']))
+                                               result[nsl]['azimuth_correction']))
+                    ned.apply(lambda st: correct_back_azimuth(None, st, result[nsl]['azimuth_correction']))
                 except:
                     logger.warning('Channel rotation failed for {}. Moving along..'.format(nsl))
                 # end try
@@ -242,7 +241,11 @@ def event_waveforms_to_rf(input_file, output_file, config, network_list='*', sta
             if(len(proc_rf_stream)):
                 for hdf_key in proc_hdfkeys[rank]:
                     # remove existing traces if there are any
-                    remove_group(output_file, hdf_key, logger=logger)
+                    try:
+                        remove_group(output_file, hdf_key, logger=logger)
+                    except Exception as e:
+                        logger.warning(str(e))
+                    # end try
 
                     logger.info("Writing RF stream(s) for {} on rank {}...".format(hdf_key, rank))
                 # end for

diff --git a/seismic/receiver_fn/rf_plot_utils.py b/seismic/receiver_fn/rf_plot_utils.py
@@ -60,6 +60,8 @@ def plot_rf_stack(rf_stream, time_window=(-10.0, 25.0), trace_height=0.2, stack_
 
     rf_stream = revert_baz(rf_stream)
 
+    if(time_window): rf_stream = rf_stream.copy().slice2(*time_window, reftime='onset')
+
     all_trace_lens = np.array([len(tr) for tr in rf_stream])
     most_common_len, _ = stats.mode(all_trace_lens, axis=None)
     stackable_stream = rf.RFStream([tr for tr in rf_stream if len(tr) == most_common_len])