Glows obs day updates (IMAP-Science-Operations-Center#1301)

* Adding flags and obs_day updates

* Updating flags, adding tests

* Fixing tests, and epoch time type

* fixing mypy

* Fixing mypy bug

* fixing merge

Author: maxinelasp

imap_processing/glows/l1a/glows_l1a.py

@@ -8,9 +8,13 @@
 
 from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes
 from imap_processing.glows.l0.decom_glows import decom_packets
-from imap_processing.glows.l0.glows_l0_data import DirectEventL0
+from imap_processing.glows.l0.glows_l0_data import DirectEventL0, HistogramL0
 from imap_processing.glows.l1a.glows_l1a_data import DirectEventL1A, HistogramL1A
-from imap_processing.spice.time import TTJ2000_EPOCH, met_to_ttj2000ns
+from imap_processing.glows.l1b.glows_l1b_data import HistogramL1B
+from imap_processing.spice.time import (
+    met_to_datetime64,
+    met_to_ttj2000ns,
+)
 
 
 def create_glows_attr_obj(data_version: str) -> ImapCdfAttributes:
@@ -66,13 +70,19 @@ def glows_l1a(packet_filepath: Path, data_version: str) -> list[xr.Dataset]:
     # Create dictionaries to group data by day
     de_by_day = process_de_l0(de_l0)
     hists_by_day = defaultdict(list)
+    # Assume the observational day starts with the first packet, then find any new
+    # observation days.
+    # TODO: replace determine_observational_day with spin table API
+    obs_days = [hist_l0[0].SEC]
+    obs_days += determine_observational_day(hist_l0)
 
-    # TODO: Make this its own function?
     for hist in hist_l0:
         hist_l1a = HistogramL1A(hist)
-        # Split by IMAP start time
-        # TODO: Should this be MET?
-        hist_day = (TTJ2000_EPOCH + met_to_ttj2000ns(hist.SEC)).astype("datetime64[D]")
+        # Determine the day the histogram belongs to. This finds the observation
+        # day in obs_day that is nearest the histogram timestamp without going over.
+        hist_day = next(
+            (day for day in reversed(obs_days) if day <= hist.SEC), obs_days[-1]
+        )
         hists_by_day[hist_day].append(hist_l1a)
 
     # Generate CDF files for each day
@@ -88,6 +98,36 @@ def glows_l1a(packet_filepath: Path, data_version: str) -> list[xr.Dataset]:
     return output_datasets
 
 
+def determine_observational_day(hist_l0: list[HistogramL0]) -> list:
+    """
+    Find the timestamps for each observational day.
+
+    This function temporarily uses the is_night flag to determine the start of a new
+    observational day, but should eventually use the spin table APIs.
+
+    Parameters
+    ----------
+    hist_l0 : list[HistogramL0]
+        List of HistogramL0 objects.
+
+    Returns
+    -------
+    list
+        List of start times for each observational day.
+    """
+    prev_is_night = -1
+    obs_day_change = []
+    for hist in hist_l0:
+        flags = HistogramL1B.deserialize_flags(hist.FLAGS)
+        is_night: int = int(flags[6])
+        if prev_is_night and not is_night:
+            obs_day_change.append(hist.SEC)
+
+        prev_is_night = is_night
+
+    return obs_day_change
+
+
 def process_de_l0(
     de_l0: list[DirectEventL0],
 ) -> dict[np.datetime64, list[DirectEventL1A]]:
@@ -111,7 +151,7 @@ def process_de_l0(
     de_by_day = dict()
 
     for de in de_l0:
-        de_day = (TTJ2000_EPOCH + met_to_ttj2000ns(de.MET)).astype("datetime64[D]")
+        de_day = (met_to_datetime64(de.MET)).astype("datetime64[D]")
         if de_day not in de_by_day:
             de_by_day[de_day] = [DirectEventL1A(de)]
         # Putting not first data int o last direct event list.
@@ -315,7 +355,7 @@ def generate_histogram_dataset(
         Dataset containing the GLOWS L1A histogram CDF output.
     """
     # Store timestamps for each HistogramL1A object.
-    time_data = np.zeros(len(hist_l1a_list), dtype="datetime64[ns]")
+    time_data = np.zeros(len(hist_l1a_list), dtype="int64")
     # TODO Add daily average of histogram counts
     # TODO compute average temperature etc
     # Data in lists, for each of the 25 time varying datapoints in HistogramL1A

imap_processing/glows/l1b/glows_l1b_data.py

@@ -611,3 +611,26 @@ def output_data(self) -> tuple:
             A tuple containing each attribute value in the class.
         """
         return tuple(getattr(self, out.name) for out in dataclasses.fields(self))
+
+    @staticmethod
+    def deserialize_flags(raw: int) -> np.ndarray[int]:
+        """
+        Deserialize the flags into a list.
+
+        Parameters
+        ----------
+        raw : int
+            16 bit integer containing the on-board flags to deserialize.
+
+        Returns
+        -------
+        flags : np.ndarray
+            Array of flags as a boolean.
+        """
+        # there are only 10 flags in the on-board flag array, additional flags are added
+        # later.
+        flags: np.ndarray[bool] = np.array(
+            [bool((raw >> i) & 1) for i in range(10)], dtype=bool
+        )
+
+        return flags

imap_processing/tests/glows/conftest.py

@@ -46,7 +46,7 @@ def l1a_dataset(packet_path):
 
 @pytest.fixture()
 def l1b_hist_dataset(l1a_dataset):
-    return glows_l1b(l1a_dataset[0], "v001")
+    return glows_l1b(l1a_dataset[1], "v001")
 
 
 @pytest.fixture()

imap_processing/tests/glows/test_glows_l1a_cdf.py

@@ -1,6 +1,5 @@
 import dataclasses
 from functools import reduce
-from pathlib import Path
 
 import numpy as np
 import pytest
@@ -16,13 +15,10 @@
 from imap_processing.glows.utils.constants import TimeTuple
 
 
-@pytest.fixture(scope="module")
-def l1a_data():
+@pytest.fixture()
+def l1a_data(packet_path):
     """Read test data from file"""
-    current_directory = Path(__file__).parent
-    packet_path = (
-        current_directory / "validation_data" / "glows_test_packet_20110921_v01.pkts"
-    )
+
     histogram_l0, de_l0 = decom_glows.decom_packets(packet_path)
 
     histogram_l1a = [HistogramL1A(hist) for hist in histogram_l0]

imap_processing/tests/glows/test_glows_l1a_data.py

@@ -8,12 +8,14 @@
 import pytest
 
 from imap_processing.glows import __version__
+from imap_processing.glows.l1a.glows_l1a import glows_l1a
 from imap_processing.glows.l1a.glows_l1a_data import (
     DirectEventL1A,
     HistogramL1A,
     StatusData,
 )
 from imap_processing.glows.utils.constants import DirectEvent, GlowsConstants, TimeTuple
+from imap_processing.spice.time import met_to_ttj2000ns
 
 
 @pytest.fixture()
@@ -38,6 +40,22 @@ def test_histogram_list(histogram_test_data, decom_test_data):
     assert sum(histogram_test_data.histogram) == histl0.EVENTS
 
 
+def test_histogram_obs_day(packet_path):
+    l1a = glows_l1a(packet_path, "v001")
+
+    assert len(l1a) == 3
+
+    assert "hist" in l1a[0].attrs["Logical_source"]
+    assert "hist" in l1a[1].attrs["Logical_source"]
+
+    # Numbers pulled from the validation data.
+    # this test assumes that the "is_night" flag switching from true to false is the
+    # start of the observation day.
+
+    assert np.array_equal(l1a[0]["imap_start_time"].data[0], 54232215.0)
+    assert np.array_equal(l1a[1]["imap_start_time"].data[0], 54232455.0)
+
+
 def test_histogram_attributes(histogram_test_data):
     """Test other data in histogram packet"""
 
@@ -488,7 +506,8 @@ def test_expected_de_results(l1a_test_data):
 
 
 def test_expected_hist_results(l1a_dataset):
-    hist = l1a_dataset[0]
+    end_time = l1a_dataset[0]["epoch"].data[-1]
+
     validation_data = (
         Path(__file__).parent / "validation_data" / "glows_l1a_hist_validation.json"
     )
@@ -521,15 +540,24 @@ def test_expected_hist_results(l1a_dataset):
         "pulse_length_variance",
     ]
 
-    for index, data in enumerate(out["output"]):
+    for data in out["output"]:
+        epoch_val = met_to_ttj2000ns(
+            TimeTuple(
+                data["imap_start_time"]["seconds"],
+                data["imap_start_time"]["subseconds"],
+            ).to_seconds()
+        )
+
+        # Validation data spans the two obs days, so this selects the correct output
+        dataset_index = 1 if epoch_val > end_time else 0
+        datapoint = l1a_dataset[dataset_index].sel(epoch=epoch_val)
+
         for field in time_fields.keys():
             expected_time = (
                 data[field]["seconds"]
                 + data[field]["subseconds"] / GlowsConstants.SUBSECOND_LIMIT
             )
-            assert np.array_equal(
-                expected_time, hist.isel(epoch=index)[time_fields[field]].data
-            )
+            assert np.array_equal(expected_time, datapoint[time_fields[field]].data)
 
         for field in compare_fields:
-            assert np.array_equal(data[field], hist.isel(epoch=index)[field].data)
+            assert np.array_equal(data[field], datapoint[field].data)

imap_processing/tests/glows/test_glows_l1b_data.py

@@ -5,7 +5,12 @@
 import pytest
 
 from imap_processing.glows.l1b.glows_l1b import glows_l1b
-from imap_processing.glows.l1b.glows_l1b_data import AncillaryParameters, DirectEventL1B
+from imap_processing.glows.l1b.glows_l1b_data import (
+    AncillaryParameters,
+    DirectEventL1B,
+    HistogramL1B,
+)
+from imap_processing.spice.time import met_to_ttj2000ns
 
 
 def test_glows_l1b_ancillary_file():
@@ -75,11 +80,8 @@ def test_glows_l1b_de():
 
 
 def test_validation_data_histogram(l1a_dataset):
-    hist_day_one = l1a_dataset[0]
-    hist_day_two = l1a_dataset[1]
-
-    l1b_day_one = glows_l1b(hist_day_one, "v001")
-    l1b_day_two = glows_l1b(hist_day_two, "v001")
+    l1b = [glows_l1b(l1a_dataset[0], "v001"), glows_l1b(l1a_dataset[1], "v001")]
+    end_time = l1b[0]["epoch"].data[-1]
 
     validation_data = (
         Path(__file__).parent
@@ -122,26 +124,23 @@ def test_validation_data_histogram(l1a_dataset):
         # "spacecraft_velocity_std_dev": "spacecraft_velocity_std_dev",
     }
 
-    for index, validation_output in enumerate(out["output"]):
-        if validation_output["imap_start_time"] < 54259215:
-            # day of 2011-09-20
-            l1b = l1b_day_one
-            l1b_index = index
-        else:
-            l1b_index = index - l1b_day_one.epoch.size
-            l1b = l1b_day_two
+    for validation_output in out["output"]:
+        epoch_val = met_to_ttj2000ns(validation_output["imap_start_time"])
+
+        # Validation data spans the two obs days, so this selects the correct output
+        dataset_index = 1 if epoch_val > end_time else 0
+        datapoint = l1b[dataset_index].sel(epoch=epoch_val)
 
         assert np.equal(
             validation_output["imap_start_time"],
-            l1b.isel(epoch=l1b_index).imap_start_time.data,
+            datapoint.imap_start_time.data,
         )
 
         for key in validation_output:
             if key not in expected_matching_columns.keys():
                 continue
-
             np.testing.assert_array_almost_equal(
-                l1b[expected_matching_columns[key]].isel(epoch=l1b_index).data,
+                datapoint[expected_matching_columns[key]].data,
                 validation_output[key],
                 decimal=1,
             )
@@ -185,3 +184,16 @@ def test_validation_data_de(l1a_dataset):
                 np.testing.assert_array_almost_equal(
                     l1b[key].isel(epoch=index).data, validation_output[key], decimal=1
                 )
+
+
+@pytest.mark.parametrize(
+    "flags, expected",
+    [
+        (0, np.zeros(10)),
+        (64, np.array([0, 0, 0, 0, 0, 0, 1, 0, 0, 0])),
+        (65, np.array([1, 0, 0, 0, 0, 0, 1, 0, 0, 0])),
+    ],
+)
+def test_deserialize_flags(flags, expected):
+    output = HistogramL1B.deserialize_flags(flags)
+    assert np.array_equal(output, expected)

imap_processing/tests/glows/test_glows_l2.py

@@ -60,7 +60,7 @@ def test_generate_l2(l1b_hist_dataset):
 
     expected_values = {
         "filter_temperature_average": [57.59],
-        "filter_temperature_std_dev": [0.23],
+        "filter_temperature_std_dev": [0.21],
         "hv_voltage_average": [1715.4],
         "hv_voltage_std_dev": [0.0],
     }