Glows badtime flag (#2760)

Author: laspsandoval

imap_processing/glows/l1b/glows_l1b_data.py

@@ -9,8 +9,15 @@
 
 from imap_processing.glows import FLAG_LENGTH
 from imap_processing.glows.utils.constants import TimeTuple
+from imap_processing.quality_flags import GLOWSL1bFlags
 from imap_processing.spice import geometry
-from imap_processing.spice.geometry import SpiceBody, SpiceFrame
+from imap_processing.spice.geometry import (
+    SpiceBody,
+    SpiceFrame,
+    frame_transform,
+    get_instrument_mounting_az_el,
+    spherical_to_cartesian,
+)
 from imap_processing.spice.spin import (
     get_instrument_spin_phase,
     get_spin_angle,
@@ -819,8 +826,11 @@ def __post_init__(
 
         # Add SPICE related variables
         self.update_spice_parameters()
-        # Will require some additional inputs
-        self.imap_spin_angle_bin_cntr = np.zeros((3600,))
+        # Calculate the spin angle bin center
+        phi = (
+            np.arange(self.number_of_bins_per_histogram, dtype=np.float64) + 0.5
+        ) / self.number_of_bins_per_histogram
+        self.imap_spin_angle_bin_cntr = phi * 360.0
 
         # TODO: This should probably be an AWS file
         # TODO Pass in AncillaryParameters object instead of reading here.
@@ -970,6 +980,79 @@ def deserialize_flags(raw: int) -> np.ndarray[int]:
 
         return flags
 
+    def flag_uv_source(self, exclusions: AncillaryExclusions) -> np.ndarray:
+        """
+        Create boolean mask where True means bin is within radius of UV source.
+
+        Parameters
+        ----------
+        exclusions : AncillaryExclusions
+            Ancillary exclusions data filtered for the current day.
+
+        Returns
+        -------
+        close_to_uv_source : np.ndarray
+            Boolean mask for uv source.
+        """
+        # Rotate spin-angle bin centers by the instrument position-angle offset
+        # so azimuth=0 aligns with the instrument pointing direction.
+        azimuth = (
+            self.imap_spin_angle_bin_cntr + self.position_angle_offset_average
+        ) % 360.0
+        # Ephemeris start time of the histogram accumulation.
+        data_start_time_et = sct_to_et(met_to_sclkticks(self.imap_start_time))
+
+        # Instrument pointing direction in the DPS frame.
+        az_el = get_instrument_mounting_az_el(SpiceFrame.IMAP_GLOWS)
+        elevation = az_el[1]
+
+        spherical = np.stack(
+            [np.ones_like(azimuth), azimuth, np.full_like(azimuth, elevation)],
+            axis=-1,
+        )  # (nbin, 3)
+
+        # Convert to unit cartesian vectors.
+        look_vecs_dps = spherical_to_cartesian(spherical)  # (nbin, 3)
+
+        # Transform unit cartesian vectors to ECLIPJ2000 frame.
+        look_vecs_ecl = frame_transform(
+            data_start_time_et,
+            look_vecs_dps,
+            SpiceFrame.IMAP_DPS,
+            SpiceFrame.ECLIPJ2000,
+        )
+
+        # UV source vectors.
+        uv_longitude = exclusions.uv_sources[
+            "ecliptic_longitude_deg"
+        ].values  # (n_src,)
+        uv_latitude = exclusions.uv_sources["ecliptic_latitude_deg"].values  # (n_src,)
+        uv_radius = np.deg2rad(
+            exclusions.uv_sources["angular_radius_for_masking"].values
+        )
+
+        uv_spherical = np.stack(
+            [np.ones_like(uv_longitude), uv_longitude, uv_latitude],
+            axis=-1,
+        )  # (n_src, 3): (r, azimuth, elevation) in degrees
+
+        uv_vecs = spherical_to_cartesian(uv_spherical)  # (n_src, 3)
+
+        # Dot product of unit vectors gives cos(separation_angle) for each
+        # histogram bin vs. each UV source -> shape (nbin, n_src).
+        # (nbin, 3) @ (3, n_src) -> (nbin, n_src)
+        # If dot product -> 1 the two vectors point in almost
+        # the same direction and needs mask.
+        # If dot product -> 0 the two directions are perpendicular on the sky.
+        cos_sep = look_vecs_ecl @ uv_vecs.T  # (nbin, n_src)
+
+        # Determine if the pixel is too close to any of the source radii.
+        close_to_uv_source = np.any(
+            cos_sep >= np.cos(uv_radius)[None, :], axis=1
+        )  # (nbin,)
+
+        return close_to_uv_source
+
     def _compute_histogram_flag_array(
         self, exclusions: AncillaryExclusions
     ) -> np.ndarray:
@@ -978,9 +1061,9 @@ def _compute_histogram_flag_array(
 
         Creates a (4, 3600) array where each row represents a different flag type:
         - Row 0: is_close_to_uv_source
-        - Row 1: is_inside_excluded_region
-        - Row 2: is_excluded_by_instr_team
-        - Row 3: is_suspected_transient
+        - Row 1: is_inside_excluded_region (TODO)
+        - Row 2: is_excluded_by_instr_team (TODO)
+        - Row 3: is_suspected_transient (TODO)
 
         Parameters
         ----------
@@ -992,5 +1075,15 @@ def _compute_histogram_flag_array(
         np.ndarray
             Array of shape (4, 3600) with bad-angle flags for each bin.
         """
-        # TODO: fill out once spice data is available
-        return np.zeros((4, 3600), dtype=np.uint8)
+        histogram_flags = np.full(
+            (4, self.number_of_bins_per_histogram),
+            GLOWSL1bFlags.NONE.value,
+            dtype=np.uint8,
+        )
+
+        close_any = self.flag_uv_source(exclusions)
+
+        # close if within radius of any UV source
+        histogram_flags[0][close_any] |= GLOWSL1bFlags.IS_CLOSE_TO_UV_SOURCE.value
+
+        return histogram_flags

imap_processing/quality_flags.py

@@ -148,6 +148,7 @@ class GLOWSL1bFlags(FlagNameMixin):
     """Glows L1b flags."""
 
     NONE = CommonFlags.NONE
+    IS_CLOSE_TO_UV_SOURCE = 2**0  # Is the bin close to a UV source.
 
 
 class ImapHiL1bDeFlags(FlagNameMixin):

imap_processing/tests/glows/conftest.py

@@ -99,17 +99,28 @@ def mock_ancillary_exclusions():
                 ["epoch", "source"],
                 [["star1", "star2", "star3"]] * len(epoch_range),
             ),
+            # degrees in [0, 360)
             "ecliptic_longitude_deg": (
                 ["epoch", "source"],
-                np.random.rand(len(epoch_range), 3),
+                np.tile(
+                    np.array([202.0812, 120.0, 250.0], dtype=np.float64),
+                    (len(epoch_range), 1),
+                ),
             ),
+            # degrees in [-90, 90]
             "ecliptic_latitude_deg": (
                 ["epoch", "source"],
-                np.random.rand(len(epoch_range), 3),
+                np.tile(
+                    np.array([18.4119, 0.0, 35.0], dtype=np.float64),
+                    (len(epoch_range), 1),
+                ),
             ),
+            # masking radius in degrees
             "angular_radius_for_masking": (
                 ["epoch", "source"],
-                np.random.rand(len(epoch_range), 3),
+                np.tile(
+                    np.array([2.0, 0.0, 0.0], dtype=np.float64), (len(epoch_range), 1)
+                ),
             ),
         },
         coords={"epoch": epoch_range},

imap_processing/tests/glows/test_glows_l1b.py

@@ -20,6 +20,7 @@
     HistogramL1B,
     PipelineSettings,
 )
+from imap_processing.spice.time import met_to_datetime64
 from imap_processing.tests.glows.conftest import mock_update_spice_parameters
 
 
@@ -33,7 +34,7 @@ def hist_dataset():
         "flags_set_onboard": np.zeros((20,)),
         "is_generated_on_ground": np.zeros((20,)),
         "number_of_spins_per_block": np.zeros((20,)),
-        "number_of_bins_per_histogram": np.zeros((20,)),
+        "number_of_bins_per_histogram": np.full((20,), 3600),
         "number_of_events": np.zeros((20,)),
         "filter_temperature_average": np.zeros((20,)),
         "filter_temperature_variance": np.zeros((20,)),
@@ -193,9 +194,11 @@ def ancillary_dict():
     return dictionary
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_histogram_mapping(
     mock_spice_function,
+    mock_flag_uv_source,
     mock_ancillary_exclusions,
     mock_ancillary_parameters,
     mock_pipeline_settings,
@@ -228,7 +231,7 @@ def test_histogram_mapping(
                 0,
                 0,
                 0,
-                0,
+                3600,
                 0,
                 encoded_val,
                 encoded_val,
@@ -255,9 +258,11 @@ def test_histogram_mapping(
     assert output[10] - expected_temp < 0.1
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_process_histogram(
     mock_spice_function,
+    mock_flag_uv_source,
     hist_dataset,
     mock_ancillary_exclusions,
     mock_ancillary_parameters,
@@ -289,7 +294,7 @@ def test_process_histogram(
         0,
         0,
         0,
-        0,
+        3600,
         0,
         encoded_val,
         encoded_val,
@@ -335,9 +340,11 @@ def test_process_de(de_dataset, ancillary_dict, mock_ancillary_parameters):
     assert np.isclose(output[8].data[0], expected_temp)
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_glows_l1b(
     mock_spice_function,
+    mock_flag_uv_source,
     de_dataset,
     hist_dataset,
     mock_ancillary_exclusions,
@@ -428,9 +435,11 @@ def test_glows_l1b(
         assert key in de_output
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_generate_histogram_dataset(
     mock_spice_function,
+    mock_flag_uv_source,
     hist_dataset,
     mock_ancillary_exclusions,
     mock_pipeline_settings,
@@ -532,6 +541,11 @@ def test_hist_spice_output(
     with furnish_kernels(kernels):
         hist_data = HistogramL1B(**params)
 
+        day = met_to_datetime64(hist_data.imap_start_time)
+        day_exclusions = mock_ancillary_exclusions.limit_by_day(day)
+
+        mask = hist_data.flag_uv_source(day_exclusions)
+
         # Assert that all these variables are the correct shape:
         assert isinstance(hist_data.spin_period_ground_average, np.float64)
         assert isinstance(hist_data.spin_period_ground_std_dev, np.float64)
@@ -543,5 +557,8 @@ def test_hist_spice_output(
         assert hist_data.spacecraft_location_std_dev.shape == (3,)
         assert hist_data.spacecraft_velocity_average.shape == (3,)
         assert hist_data.spacecraft_velocity_std_dev.shape == (3,)
+        assert mask.shape == (3600,)
+        # For 2 degree radius: 20 + 20 + 1(center) ≈ 41 bins.
+        assert np.count_nonzero(mask) == 41
 
         # TODO: Maxine will validate actual data with GLOWS team

imap_processing/tests/glows/test_glows_l1b_data.py

@@ -83,9 +83,11 @@ def test_glows_l1b_de():
     assert np.allclose(pulse_len, expected_pulse)
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_validation_data_histogram(
     mock_spice_function,
+    mock_flag_uv_source,
     l1a_dataset,
     mock_ancillary_exclusions,
     mock_pipeline_settings,

imap_processing/tests/glows/test_glows_l2.py

@@ -35,9 +35,11 @@ def l1b_hists():
     return input
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_glows_l2(
     mock_spice_function,
+    mock_flag_uv_source,
     l1a_dataset,
     mock_ancillary_exclusions,
     mock_pipeline_settings,
@@ -60,9 +62,11 @@ def test_glows_l2(
     assert np.allclose(l2["filter_temperature_average"].values, [57.6], rtol=0.1)
 
 
+@patch.object(HistogramL1B, "flag_uv_source", return_value=np.zeros(3600, dtype=bool))
 @patch.object(HistogramL1B, "update_spice_parameters", autospec=True)
 def test_generate_l2(
     mock_spice_function,
+    mock_flag_uv_source,
     l1a_dataset,
     mock_ancillary_exclusions,
     mock_pipeline_settings,