From aa8f6ab2b9c70665a9c04e383aa422a2e979507a Mon Sep 17 00:00:00 2001
From: Alex Wernle <awernle@usgs.gov>
Date: Wed, 2 Nov 2022 11:37:15 -0600
Subject: [PATCH] Renamed variables and put into paragraphs. Asserted that the
 count correctly shows where there are fewer input data.

---
 geomagio/algorithm/AverageAlgorithm.py       | 46 ++++++++++----------
 test/algorithm_test/AverageAlgorithm_test.py |  2 +
 2 files changed, 25 insertions(+), 23 deletions(-)

diff --git a/geomagio/algorithm/AverageAlgorithm.py b/geomagio/algorithm/AverageAlgorithm.py
index 16eafdf98..2b07bda4d 100644
--- a/geomagio/algorithm/AverageAlgorithm.py
+++ b/geomagio/algorithm/AverageAlgorithm.py
@@ -111,30 +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.nanmean(combined, axis=0)
-
-        # Stack timeseries and count non-nan values per timestep
+        # 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)
-        channel_count = numpy.count_nonzero(~numpy.isnan(ts_counter), axis=1)
-
-        # create metadata
-        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
-        new_stats2 = new_stats.copy()
-        new_stats2.channel = new_stats.channel + "_count"
-
-        # create streams from the trace function
-        stream = obspy.core.Stream((obspy.core.Trace(dst_tot, new_stats),))
-        stream += obspy.core.Trace(numpy.asarray(channel_count), new_stats2)
-
-        # return averaged values as a stream
+        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
diff --git a/test/algorithm_test/AverageAlgorithm_test.py b/test/algorithm_test/AverageAlgorithm_test.py
index 88ca37e0d..2480f655b 100644
--- a/test/algorithm_test/AverageAlgorithm_test.py
+++ b/test/algorithm_test/AverageAlgorithm_test.py
@@ -82,6 +82,8 @@ def test_gaps():
 
     # 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():
-- 
GitLab