diff --git a/imap_processing/glows/l2/glows_l2_data.py b/imap_processing/glows/l2/glows_l2_data.py index c91b797f9..46a91ad2e 100644 --- a/imap_processing/glows/l2/glows_l2_data.py +++ b/imap_processing/glows/l2/glows_l2_data.py @@ -370,10 +370,20 @@ def __init__( """ active_flags = np.array(pipeline_settings.active_bad_time_flags, dtype=float) + # Apply sunrise/sunset offsets to extend the night region around + # is_night transitions before selecting good blocks. + flags = self.apply_is_night_offsets( + l1b_dataset["flags"].data, + is_night_idx=6, # is_night is the 7th bad-time flag (0-indexed) + sunrise_offset=pipeline_settings.sunrise_offset, + sunset_offset=pipeline_settings.sunset_offset, + ) + flags_da = xr.DataArray(flags, dims=l1b_dataset["flags"].dims) + # Select the good blocks (i.e. epoch values) according to the flags. Drop any # bad blocks before processing. good_data = l1b_dataset.isel( - epoch=self.return_good_times(l1b_dataset["flags"], active_flags) + epoch=self.return_good_times(flags_da, active_flags) ) # TODO: bad angle filter # TODO: filter bad bins out. Needs to happen here while everything is still @@ -533,6 +543,97 @@ def return_good_times(flags: xr.DataArray, active_flags: NDArray) -> NDArray: good_times = np.where(np.all(flags[:, active_flags == 1] == 1, axis=1))[0] return good_times + @staticmethod + def apply_is_night_offsets( + flags: np.ndarray, + is_night_idx: int, + sunrise_offset: float, + sunset_offset: float, + ) -> np.ndarray: + """ + Apply sunrise/sunset offsets to is_night transitions. + + Per algorithm doc v4.4.7, Sec. 3.9.1, item 2 (raw is_night: 1=night, 0=day): + + sunset_offset applies at both transitions: + >0: night shortens by N at each end (first N night epochs at sunset become + day; last N night epochs before sunrise become day) + <0: night extends by |N| at each end + + sunrise_offset is an additional adjustment at sunrise (is_night 1->0) only: + >0: night extends N histograms past the raw sunrise transition + <0: night shortens by |N| before the raw sunrise transition + + In the processed flags array: 0 = bad (night), 1 = good (day). + + Parameters + ---------- + flags : numpy.ndarray + Flags array with shape (n_epochs, FLAG_LENGTH), 0=bad, 1=good. + is_night_idx : int + Column index of the is_night flag in the flags array. + sunrise_offset : float + Additional histogram shift at the sunrise (is_night 1->0) transition. + sunset_offset : float + Histogram shift applied at both the sunset and sunrise transitions. + + Returns + ------- + numpy.ndarray + Copy of flags with the is_night column adjusted per the offsets. + + Notes + ----- + Algorithm doc v4.4.7, Sec. 3.9.1, item 2 + is_night: 1 = daytime (good), 0 = night (bad) + """ + flags_with_offsets = flags.copy() + + # If sunrise_offset=0 and sunset_offset=0 then no corrections are needed + # relative to is_night transition set onboard. + if sunrise_offset == 0 and sunset_offset == 0: + return flags_with_offsets + + is_night_col = flags[:, is_night_idx] + n = flags.shape[0] + diff = np.diff(is_night_col.astype(int)) + sunset_index = np.where(diff == -1)[0] + sunrise_index = np.where(diff == 1)[0] + + if sunrise_offset > 0: + # Night (flag = 0) extends by sunrise_offset relative + # to is_night 0 -> 1 transition. + for i in sunrise_index: + flags_with_offsets[ + i + 1 : min(n, i + 1 + int(sunrise_offset)), is_night_idx + ] = 0 + + if sunrise_offset < 0: + # Night (flag = 0) shortens by sunrise_offset relative + # to is_night 0 -> 1 transition. + for i in sunrise_index: + flags_with_offsets[ + max(0, i + 1 + int(sunrise_offset)) : i + 1, is_night_idx + ] = 1 + + if sunset_offset > 0: + # Night (flag = 0) shortens by sunset_offset relative + # to is_night 1 -> 0 transition. + for i in sunset_index: + flags_with_offsets[ + i + 1 : min(n, i + 1 + int(sunset_offset)), is_night_idx + ] = 1 + + if sunset_offset < 0: + # Night (flag = 0) extends by sunset_offset relative + # to is_night 1 -> 0 transition. + for i in sunset_index: + flags_with_offsets[ + max(0, i + 1 + int(sunset_offset)) : i + 1, is_night_idx + ] = 0 + + return flags_with_offsets + def compute_position_angle(self) -> float: """ Compute the position angle based on the instrument mounting. diff --git a/imap_processing/tests/glows/test_glows_l2_data.py b/imap_processing/tests/glows/test_glows_l2_data.py index cb7a71ac6..028c71e80 100644 --- a/imap_processing/tests/glows/test_glows_l2_data.py +++ b/imap_processing/tests/glows/test_glows_l2_data.py @@ -342,6 +342,38 @@ def test_filter_good_times(): assert np.array_equal(good_times, expected_good_times) +@pytest.mark.parametrize( + "sunrise_offset, sunset_offset, expected_is_night", + [ + # sunset>0 shortens at sunset; sunrise>0 extends at sunrise + (1, 1, [1, 1, 1, 1, 0, 0, 0, 1]), + # sunset>0 shortens at sunset; sunrise<0 shortens at sunrise + (-1, 1, [1, 1, 1, 1, 0, 1, 1, 1]), + # sunset<0 extends at sunset; sunrise>0 extends at sunrise + (1, -1, [1, 1, 0, 0, 0, 0, 0, 1]), + # sunset<0 extends at sunset; sunrise<0 shortens at sunrise + (-1, -1, [1, 1, 0, 0, 0, 1, 1, 1]), + # zero offsets: no change + (0, 0, [1, 1, 1, 0, 0, 0, 1, 1]), + ], +) +def test_apply_is_night_offsets(sunrise_offset, sunset_offset, expected_is_night): + """Test apply_is_night_offsets function.""" + + # Setup: epochs 0-2 day, 3-5 night, 6-7 day (processed flags: 0=night, 1=day). + flags = np.ones((8, 17), dtype=float) + flags[3:6, 6] = 0 # epochs 3-5 are night + + result = HistogramL2.apply_is_night_offsets( + flags, + is_night_idx=6, + sunrise_offset=sunrise_offset, + sunset_offset=sunset_offset, + ) + + assert np.array_equal(result[:, 6], np.array(expected_is_night, dtype=float)) + + # ── spin_angle tests ────────────────────────────────────────────────────────── @@ -392,7 +424,8 @@ def test_compute_position_angle(): def l1b_dataset_full(): """Minimal L1B dataset with all variables required by HistogramL2. - Two epochs, four bins, 17 flags (all good). + Two epochs, four bins, 17 flags. Both epochs are daytime (is_night=1). + All other flags are 1 (good). """ n_epochs, n_bins, n_angle_flags, n_time_flags = 2, 4, 4, 17 fillval = GlowsConstants.HISTOGRAM_FILLVAL @@ -402,6 +435,9 @@ def l1b_dataset_full(): spin_angle = np.tile(np.linspace(0, 270, n_bins), (n_epochs, 1)) histogram_flag_array = np.zeros((n_epochs, n_angle_flags, n_bins), dtype=np.uint8) + # All flags good (1). Index 6 is is_night: 1 = daytime (good). + flags = np.ones((n_epochs, n_time_flags), dtype=float) + return xr.Dataset( { "histogram": (["epoch", "bins"], histogram), @@ -413,10 +449,7 @@ def l1b_dataset_full(): histogram_flag_array, ), "number_of_bins_per_histogram": (["epoch"], [n_bins, n_bins]), - "flags": ( - ["epoch", "flag_index"], - np.ones((n_epochs, n_time_flags)), - ), + "flags": (["epoch", "flag_index"], flags), "filter_temperature_average": (["epoch"], [20.0, 21.0]), "hv_voltage_average": (["epoch"], [1000.0, 1000.0]), "pulse_length_average": (["epoch"], [5.0, 5.0]),