GLOWS - Propagate bad-angle flags from L1B Histogram to L2 in GLOWS

imap_processing/glows/l2/glows_l2_data.py

@@ -90,7 +90,20 @@ def __post_init__(self, l1b_data: xr.Dataset) -> None:
         # TODO: Average this, or should they all be the same?
         self.spin_angle = np.average(l1b_data["imap_spin_angle_bin_cntr"].data, axis=0)
 
-        self.histogram_flag_array = np.zeros(self.number_of_bins)
+        # Apply 'OR' operation to histogram_flag_array across all
+        # good-time L1B blocks per bin.
+        # Per Section 12.3.4: a flag is True in L2 if it is True in any L1B block.
+        # flags shape: (n_epochs, 4, n_bins)
+        flags = l1b_data["histogram_flag_array"].data
+        if flags.size > 0:
+            # Flatten epochs and flag rows into one axis: (n_epochs * 4, n_bins)
+            flags_2d = flags.reshape(-1, self.number_of_bins)
+            # Apply binary 'OR' operation across all rows per bin: (n_bins,)
+            self.histogram_flag_array = np.bitwise_or.reduce(flags_2d, axis=0).astype(
+                np.uint8
+            )
+        else:
+            self.histogram_flag_array = np.zeros(self.number_of_bins, dtype=np.uint8)
         self.ecliptic_lon = np.zeros(self.number_of_bins)
         self.ecliptic_lat = np.zeros(self.number_of_bins)
 

imap_processing/tests/glows/test_glows_l2_data.py

@@ -54,20 +54,26 @@ def l1b_dataset():
 
     Two timestamps, four bins.
     Bin 3 is masked (-1) at timestamp 0.
+    histogram_flag_array has shape (epoch, bad_angle_flags, bins) with all zeros.
     """
-    n_epochs, n_bins = 2, 4
+    n_epochs, n_bins, n_flags = 2, 4, 4
     epoch = xr.DataArray(np.arange(n_epochs), dims=["epoch"])
     bins = xr.DataArray(np.arange(n_bins), dims=["bins"])
 
     histogram = np.array([[10, 20, 30, -1], [10, 20, 30, 40]], dtype=float)
     spin_angle = np.tile(np.linspace(0, 270, n_bins), (n_epochs, 1))
+    histogram_flag_array = np.zeros((n_epochs, n_flags, n_bins), dtype=np.uint8)
 
     ds = xr.Dataset(
         {
             "histogram": (["epoch", "bins"], histogram),
             "spin_period_average": (["epoch"], [15.0, 15.0]),
             "number_of_spins_per_block": (["epoch"], [5, 5]),
             "imap_spin_angle_bin_cntr": (["epoch", "bins"], spin_angle),
+            "histogram_flag_array": (
+                ["epoch", "bad_angle_flags", "bins"],
+                histogram_flag_array,
+            ),
         },
         coords={"epoch": epoch, "bins": bins},
     )
@@ -108,16 +114,21 @@ def test_zero_exposure_bins():
     when all histogram values are masked (-1). Flux and uncertainty are zero
     because the raw histogram sums are zero.
     """
-    n_epochs, n_bins = 2, 3
+    n_epochs, n_bins, n_flags = 2, 3, 4
     histogram = np.full((n_epochs, n_bins), -1, dtype=float)
     spin_angle = np.tile(np.linspace(0, 240, n_bins), (n_epochs, 1))
+    histogram_flag_array = np.zeros((n_epochs, n_flags, n_bins), dtype=np.uint8)
 
     ds = xr.Dataset(
         {
             "histogram": (["epoch", "bins"], histogram),
             "spin_period_average": (["epoch"], [15.0, 15.0]),
             "number_of_spins_per_block": (["epoch"], [5, 5]),
             "imap_spin_angle_bin_cntr": (["epoch", "bins"], spin_angle),
+            "histogram_flag_array": (
+                ["epoch", "bad_angle_flags", "bins"],
+                histogram_flag_array,
+            ),
         },
         coords={"epoch": xr.DataArray(np.arange(n_epochs), dims=["epoch"])},
     )
@@ -138,6 +149,72 @@ def test_number_of_bins(l1b_dataset):
     assert len(lc.exposure_times) == 4
 
 
+def test_histogram_flag_array_or_propagation():
+    """histogram_flag_array is OR'd across all L1B epochs and flag rows per bin.
+
+    Per Section 12.3.4: a flag is True in L2 if it is True in any L1B block.
+    """
+    n_epochs, n_bins, n_flags = 3, 4, 4
+    histogram = np.ones((n_epochs, n_bins), dtype=float)
+    spin_angle = np.tile(np.linspace(0, 270, n_bins), (n_epochs, 1))
+
+    # epoch 0, flag row 0 (IS_CLOSE_TO_UV_SOURCE=1): bin 0 flagged
+    # epoch 1, flag row 1 (IS_INSIDE_EXCLUDED_REGION=2): bin 2 flagged
+    # epoch 2: no flags set
+    histogram_flag_array = np.zeros((n_epochs, n_flags, n_bins), dtype=np.uint8)
+    histogram_flag_array[0, 0, 0] = 1  # IS_CLOSE_TO_UV_SOURCE on bin 0
+    histogram_flag_array[1, 1, 2] = 2  # IS_INSIDE_EXCLUDED_REGION on bin 0, 2
+    histogram_flag_array[1, 0, 0] = 2
+
+    ds = xr.Dataset(
+        {
+            "histogram": (["epoch", "bins"], histogram),
+            "spin_period_average": (["epoch"], [15.0, 15.0, 15.0]),
+            "number_of_spins_per_block": (["epoch"], [5, 5, 5]),
+            "imap_spin_angle_bin_cntr": (["epoch", "bins"], spin_angle),
+            "histogram_flag_array": (
+                ["epoch", "bad_angle_flags", "bins"],
+                histogram_flag_array,
+            ),
+        },
+        coords={"epoch": xr.DataArray(np.arange(n_epochs), dims=["epoch"])},
+    )
+    lc = DailyLightcurve(ds)
+
+    assert (
+        lc.histogram_flag_array[0] == 3
+    )  # IS_CLOSE_TO_UV_SOURCE and IS_INSIDE_EXCLUDED_REGION
+    assert lc.histogram_flag_array[1] == 0  # no flags on bin 1
+    assert lc.histogram_flag_array[2] == 2  # IS_INSIDE_EXCLUDED_REGION
+    assert lc.histogram_flag_array[3] == 0  # no flags on bin 3
+
+
+def test_histogram_flag_array_zero_epochs():
+    """histogram_flag_array is all zeros when the input dataset is empty."""
+    n_bins, n_flags = 4, 4
+    histogram = np.empty((0, n_bins), dtype=float)
+    spin_angle = np.empty((0, n_bins), dtype=float)
+    histogram_flag_array = np.empty((0, n_flags, n_bins), dtype=np.uint8)
+
+    ds = xr.Dataset(
+        {
+            "histogram": (["epoch", "bins"], histogram),
+            "spin_period_average": (["epoch"], []),
+            "number_of_spins_per_block": (["epoch"], []),
+            "imap_spin_angle_bin_cntr": (["epoch", "bins"], spin_angle),
+            "histogram_flag_array": (
+                ["epoch", "bad_angle_flags", "bins"],
+                histogram_flag_array,
+            ),
+        },
+        coords={"epoch": xr.DataArray(np.arange(0), dims=["epoch"])},
+    )
+    lc = DailyLightcurve(ds)
+
+    assert len(lc.histogram_flag_array) == n_bins
+    assert np.all(lc.histogram_flag_array == 0)
+
+
 def test_filter_good_times():
     """Epochs where any active flag is 0 are excluded; inactive flags are ignored."""
     active_flags = np.ones((17,))