Merge branch 'dst_changes' into 'master'

Changed dst_tot to calculate average regardless of nan values and added a new count

Closes #78

See merge request !194

geomagio/algorithm/AverageAlgorithm.py

@@ -63,18 +63,14 @@ class AverageAlgorithm(Algorithm):
             if ts.stats.npts != self._npts:
                 raise AlgorithmException("Received timeseries have different lengths")
 
-            if numpy.isnan(ts.data).all():
-                raise AlgorithmException(
-                    "Trace for %s observatory is completely empty." % (ts.stats.station)
-                )
-
             if ts.stats.starttime != self._stt:
                 raise AlgorithmException(
                     "Received timeseries have different starttimes"
                 )
 
     def process(self, timeseries):
-        """averages a channel across multiple stations
+        """averages a channel across multiple stations and counts available
+        observatories per timestep
 
         Parameters
         ----------
@@ -82,7 +78,7 @@ class AverageAlgorithm(Algorithm):
         Returns
         -------
         out_stream:
-            new stream object containing the averaged values.
+            new stream object containing the averaged values and available observatories.
         """
         self.observatories = self.observatories or [t.stats.station for t in timeseries]
 
@@ -115,21 +111,30 @@ class AverageAlgorithm(Algorithm):
             ts = timeseries.select(station=obsy)[0]
             combined.append(ts.data * latcorr)
 
-        # after looping over stations, compute average
-        dst_tot = numpy.mean(combined, axis=0)
-
-        # Create a stream from the trace function
-        new_stats = obspy.core.Stats()
-        new_stats.station = "USGS"
-        new_stats.channel = self.outchannel
-        new_stats.network = "NT"
-        new_stats.location = self.outlocation
-        new_stats.starttime = timeseries[0].stats.starttime
-        new_stats.npts = timeseries[0].stats.npts
-        new_stats.delta = timeseries[0].stats.delta
-        stream = obspy.core.Stream((obspy.core.Trace(dst_tot, new_stats),))
-
-        # return averaged values as a stream
+        # calculate averages
+        average_data = numpy.nanmean(combined, axis=0)
+        average_stats = obspy.core.Stats()
+        average_stats.station = "USGS"
+        average_stats.channel = self.outchannel
+        average_stats.network = "NT"
+        average_stats.location = self.outlocation
+        average_stats.starttime = timeseries[0].stats.starttime
+        average_stats.npts = timeseries[0].stats.npts
+        average_stats.delta = timeseries[0].stats.delta
+
+        # calculate counts
+        ts_counter = numpy.column_stack(timeseries)
+        count_data = numpy.count_nonzero(~numpy.isnan(ts_counter), axis=1)
+        count_stats = average_stats.copy()
+        count_stats.channel = average_stats.channel + "_count"
+
+        # return stream
+        stream = obspy.core.Stream(
+            (
+                obspy.core.Trace(average_data, average_stats),
+                obspy.core.Trace(numpy.asarray(count_data), count_stats),
+            )
+        )
         return stream
 
     @classmethod

pyproject.toml

@@ -53,7 +53,7 @@ mypy = "^0.982"
 poethepoet = "^0.16.4"
 pytest = "^7.1.3"
 pytest-cov = "^4.0.0"
-safety = "^1.10.3"
+safety = "^2.2.0"
 # type hints
 data-science-types = "^0.2.23"
 openpyxl-stubs = "^0.1.24"

test/algorithm_test/AverageAlgorithm_test.py

@@ -8,18 +8,19 @@ from numpy.testing import assert_array_equal, assert_equal
 
 def test_process():
     """AverageAlgorithm_test.test_process()
-    confirms that the output of the algorithm is the average of three
-    different stations.
+    confirms that the output of the algorithm is the average of two
+    stations and ignores the third.
 
     """
-    # Make three test data arrays of 1's, 3's, and 5's
-    # the average of the three arrays should be 3
-    test_data1 = np.ones(5)
+    # Make three test data arrays of NAN's, 3's, and 5's
+    # the average of the three arrays should be 4
+    test_data1 = np.zeros(5)
+    test_data1[:] = np.nan
     test_data2 = np.ones(5) * 5
     test_data3 = np.ones(5) * 3
 
-    # Create expected solution array of 3's
-    expected_solution = np.ones(5) * 3
+    # Create expected solution array of 4's
+    expected_solution = np.ones(5) * 4
 
     # Create timeseries with first trace that uses test_data1
     timeseries = Stream()
@@ -79,8 +80,10 @@ def test_gaps():
     # Run timeseries through the average process
     outstream = alg.process(timeseries)
 
-    # The gaps should not be changed and should remain 'nan'
-    assert_array_equal(outstream[0].data, [1, 1, np.nan, np.nan, 1, 1])
+    # The gaps should not be ignored
+    assert_array_equal(outstream[0].data, [1, 1, 1, 1, 1, 1])
+    # Counts should show where there are fewer input data
+    assert_array_equal(outstream[1].data, [2, 2, 1, 1, 2, 2])
 
 
 def test_metadata():