[new] mass channel surveillance added

parameters.yaml

@@ -3,7 +3,8 @@ datapath: "/data/SDS/"                    # SC3 Datapath
 networks: ["1Y", "HA"]                    # select networks, list or str
 stations: "*"                             # select stations, list or str
 locations: "*"                            # select locations, list or str
-channels: ["EX1", "EX2", "EX3", "VEI"]    # Specify SOH channels, currently supported EX[1-3] and VEI
+channels: ["EX1", "EX2", "EX3", "VEI",
+           "VM1", "VM2", "VM3"]    # Specify SOH channels, currently supported EX[1-3], VEI and VM[1-3]
 stations_blacklist: ["TEST", "EREA"]      # exclude these stations
 networks_blacklist: []                    # exclude these networks
 interval: 60               # Perform checks every x seconds
@@ -43,6 +44,7 @@ THRESHOLDS:
   low_volt: 12                 # min voltage for low voltage warning
   high_volt: 14.8               # max voltage for over voltage warning
   unclassified: 5              # min voltage samples not classified for warning
+  max_vm: [1.5, 2.5]           # thresholds for mass offset (warn, fail)
 
 # ---------------------------------------- OPTIONAL PARAMETERS ---------------------------------------------------------
 
@@ -62,13 +64,23 @@ EMAIL:
   networks_blacklist: []                           # do not send emails for specific network
 
 # names for plotting of the above defined parameter "channels" in the same order
-channel_names: ["Temperature (°C)", "230V/12V Status (V)", "Router/Charger State (V)", "Logger Voltage (V)"]
+channel_names: ["Temperature (°C)",
+                "230V/12V (V)",
+                "Rout/Charge (V)",
+                "Logger (V)",
+                "Mass 1 (V)",
+                "Mass 2 (V)",
+                "Mass 3 (V)"]
+
 # specify y-ticks (and ylims) giving, (ymin, ymax, step) for each of the above channels (0: default)
 CHANNEL_TICKS:
   - [-10, 50, 10]
   - [1, 5, 1]
   - [1, 5, 1]
   - [9, 15, 1]
+  - [-2, 2, 1]
+  - [-2, 2, 1]
+  - [-2, 2, 1]
 
 # Factor for channel to SI-units (for plotting)
 CHANNEL_UNITS:
@@ -76,6 +88,9 @@ CHANNEL_UNITS:
   EX2: 1e-6
   EX3: 1e-6
   VEI: 1e-3
+  VM1: 1e-6
+  VM2: 1e-6
+  VM3: 1e-6
 
 # Transform channel for plotting, perform arithmetic operations in given order, e.g.: PBox EX1 V to deg C: 20 * x -20
 CHANNEL_TRANSFORM:

survBot.py

@@ -61,7 +61,7 @@ def fancy_timestr(dt, thresh=600, modif='+'):
 
 class SurveillanceBot(object):
     def __init__(self, parameter_path, outpath_html=None):
-        self.keys = ['last active', '230V', '12V', 'router', 'charger', 'voltage', 'temp', 'other']
+        self.keys = ['last active', '230V', '12V', 'router', 'charger', 'voltage', 'mass', 'temp', 'other']
         self.parameter_path = parameter_path
         self.update_parameters()
         self.starttime = UTCDateTime()
@@ -242,7 +242,7 @@ class SurveillanceBot(object):
     def get_station_delay(self, nwst_id):
         """ try to get station delay from SDS archive using client"""
         locations = ['', '0', '00']
-        channels = ['HHZ', 'HHE', 'HHN', 'VEI', 'EX1', 'EX2', 'EX3']
+        channels = ['HHZ', 'HHE', 'HHN'] + self.parameters.get('channels')
         network, station = nwst_id.split('.')[:2]
 
         times = []
@@ -685,6 +685,7 @@ class StationQC(object):
         self.pb_temp_analysis()
         self.pb_power_analysis()
         self.pb_rout_charge_analysis()
+        self.mass_analysis()
 
     def return_print_analysis(self):
         items = [self.nwst_id]
@@ -718,7 +719,8 @@ class StationQC(object):
         key = 'voltage'
         st = self.stream.select(channel=channel)
         trace = self.get_trace(st, key)
-        if not trace: return
+        if not trace:
+            return
         voltage = trace.data * 1e-3
         low_volt = self.parameters.get('THRESHOLDS').get('low_volt')
         high_volt = self.parameters.get('THRESHOLDS').get('high_volt')
@@ -756,14 +758,15 @@ class StationQC(object):
         key = 'temp'
         st = self.stream.select(channel=channel)
         trace = self.get_trace(st, key)
-        if not trace: return
+        if not trace:
+            return
         voltage = trace.data * 1e-6
         thresholds = self.parameters.get('THRESHOLDS')
         temp = 20. * voltage - 20
         # average temp
         timespan = min([self.parameters.get('timespan') * 24 * 3600, int(len(temp) / trace.stats.sampling_rate)])
         nsamp_av = int(trace.stats.sampling_rate) * timespan
-        av_temp_str = str(round(np.mean(temp[-nsamp_av:]), 1)) + deg_str
+        av_temp_str = str(round(np.nanmean(temp[-nsamp_av:]), 1)) + deg_str
         # dt of average
         dt_t_str = str(timedelta(seconds=int(timespan))).replace(', 0:00:00', '')
         # current temp
@@ -771,7 +774,7 @@ class StationQC(object):
         if self.verbosity > 1:
             self.print(40 * '-')
             self.print('Performing PowBox temperature check (EX1)', flush=False)
-            self.print(f'Average temperature at {np.mean(temp)}\N{DEGREE SIGN}', flush=False)
+            self.print(f'Average temperature at {np.nanmean(temp)}\N{DEGREE SIGN}', flush=False)
             self.print(f'Peak temperature at {max(temp)}\N{DEGREE SIGN}', flush=False)
             self.print(f'Min temperature at {min(temp)}\N{DEGREE SIGN}', flush=False)
         max_temp = thresholds.get('max_temp')
@@ -787,6 +790,52 @@ class StationQC(object):
                            status_message=cur_temp,
                            detailed_message=f'Average temperature of last {dt_t_str}: {av_temp_str}')
 
+    def mass_analysis(self, channels=('VM1', 'VM2', 'VM3'), n_samp_mean=10):
+        """ Analyse datalogger mass channels. """
+        key = 'mass'
+
+        # build stream with all channels
+        st = Stream()
+        for channel in channels:
+            st += self.stream.select(channel=channel).copy()
+        st.merge()
+
+        # return if there are no three components
+        if not len(st) == 3:
+            return
+
+        # correct for channel unit
+        for trace in st:
+            trace.data = trace.data * 1e-6  # hardcoded, change this?
+
+        # calculate average of absolute maximum of mass offset of last n_samp_mean
+        last_values = np.array([trace.data[-n_samp_mean:] for trace in st])
+        last_val_mean = np.nanmean(last_values, axis=1)
+        common_highest_val = np.nanmax(abs(last_val_mean))
+        common_highest_val = round(common_highest_val, 1)
+
+        # get thresholds for WARN (max_vm1) and FAIL (max_vm2)
+        thresholds = self.parameters.get('THRESHOLDS')
+        max_vm = thresholds.get('max_vm')
+        if not max_vm:
+            return
+        max_vm1, max_vm2 = max_vm
+
+        # change status depending on common_highest_val
+        if common_highest_val < max_vm1:
+            self.status_ok(key, detailed_message=f'{common_highest_val}V')
+        elif max_vm1 <= common_highest_val < max_vm2:
+            self.warn(key=key,
+                      detailed_message=f'Warning raised for mass centering. Highest val {common_highest_val}V', )
+        else:
+            self.error(key=key,
+                      detailed_message=f'Fail status for mass centering. Highest val {common_highest_val}V',)
+
+        if self.verbosity > 1:
+            self.print(40 * '-')
+            self.print('Performing mass position check', flush=False)
+            self.print(f'Average mass position at {common_highest_val}', flush=False)
+
     def pb_power_analysis(self, channel='EX2', pb_dict_key='pb_SOH2'):
         """ Analyse EX2 channel of PowBox """
         keys = ['230V', '12V']
@@ -1058,6 +1107,23 @@ class StatusOther(Status):
         return message, detailed_message
 
 
+def common_mass_trace(st):
+    traces = st.traces
+    if not len(traces) == 3:
+        return
+    check_keys = ['sampling_rate', 'starttime', 'endtime', 'npts']
+    # check if values of the above keys are identical for all traces
+    for c_key in check_keys:
+        if not traces[0].stats.get(c_key) == traces[1].stats.get(c_key) == traces[2].stats.get(c_key):
+            return
+    max_1_2 = np.fmax(abs(traces[0]), abs(traces[1]))
+    abs_max = np.fmax(max_1_2, abs(traces[2]))
+    return_trace = traces[0].copy()
+    return_trace.data = abs_max
+    return return_trace
+
+
+
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='Call survBot')
     parser.add_argument('-html', dest='html_path', default=None, help='filepath for HTML output')

utils.py

@@ -10,9 +10,11 @@ def get_bg_color(check_key, status, dt_thresh=None, hex=False):
         bg_color = get_time_delay_color(message, dt_thresh)
     elif check_key == 'temp':
         bg_color = get_temp_color(message)
+    elif check_key == 'mass':
+        bg_color = get_mass_color(message)
     else:
         if status.is_warn:
-            bg_color = get_color('WARNX')(status.count)
+            bg_color = get_warn_color(status.count)
         elif status.is_error:
             bg_color = get_color('FAIL')
         else:
@@ -30,7 +32,6 @@ def get_color(key):
     colors_dict = {'FAIL': (255, 50, 0, 255),
                    'NO DATA': (255, 255, 125, 255),
                    'WARN': (255, 255, 80, 255),
-                   'WARNX': lambda x: (min([255, 200 + x ** 2]), 255, 80, 255),
                    'OK': (125, 255, 125, 255),
                    'undefined': (230, 230, 230, 255)}
     return colors_dict.get(key)
@@ -45,6 +46,18 @@ def get_time_delay_color(dt, dt_thresh):
     return get_color('FAIL')
 
 
+def get_warn_color(count):
+    color = (min([255, 200 + count ** 2]), 255, 80, 255)
+    return color
+
+
+def get_mass_color(message):
+    # can change this to something else if wanted. This way it always returns get_color (without warn count)
+    if isinstance(message, (float, int)):
+        return get_color('OK')
+    return get_color(message)
+
+
 def get_temp_color(temp, vmin=-10, vmax=60, cmap='coolwarm'):
     """ Get an rgba temperature value back from specified cmap, linearly interpolated between vmin and vmax. """
     if type(temp) in [str]:
@@ -60,9 +73,8 @@ def modify_stream_for_plot(st, parameters):
     ch_units = parameters.get('CHANNEL_UNITS')
     ch_transf = parameters.get('CHANNEL_TRANSFORM')
 
-    # if either of both are defined make copy
+    # make a copy
-    if ch_units or ch_transf:
+    st = st.copy()
-        st = st.copy()
 
     # modify trace for plotting by multiplying unit factor (e.g. 1e-3 mV to V)
     if ch_units:
@@ -71,6 +83,7 @@ def modify_stream_for_plot(st, parameters):
             unit_factor = ch_units.get(channel)
             if unit_factor:
                 tr.data = tr.data * float(unit_factor)
+
     # modify trace for plotting by other arithmetic expressions
     if ch_transf:
         for tr in st:
@@ -79,6 +92,10 @@ def modify_stream_for_plot(st, parameters):
             if transf:
                 tr.data = transform_trace(tr.data, transf)
 
+    # change channel IDs to prevent re-sorting in obspy routine
+    for index, trace in enumerate(st):
+        trace.id = f'trace {index + 1}: {trace.id}'
+
     return st
 
 
@@ -142,4 +159,4 @@ def trace_yticks(fig, parameters, verbosity=0):
 
         yticks = list(range(ymin, ymax + step, step))
         ax.set_yticks(yticks)
-        ax.set_ylim(ymin - step, ymax + step)
+        ax.set_ylim(ymin - 0.33 * step, ymax + 0.33 * step)