Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
103 changes: 102 additions & 1 deletion imap_processing/glows/l2/glows_l2_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will this always be a value of 6? It might be worth adding this as a constant to the GlowsConstants class

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
Expand Down Expand Up @@ -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
Comment on lines +592 to +595
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since offsets don't get applied to the flags in this case, you could do this check on the flags array before copying it to flags_with_offsets. This way you can avoid an unnecessary copy of the data and make it more clear that the original flags is returned without corrections.


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.
Expand Down
43 changes: 38 additions & 5 deletions imap_processing/tests/glows/test_glows_l2_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]),
],
Comment on lines +346 to +358
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: since the parameters are in order of sunrise_offset then sunset_offset, consider updating your comments to follow this order as well for clarity

)
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 ──────────────────────────────────────────────────────────


Expand Down Expand Up @@ -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
Expand All @@ -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),
Expand All @@ -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]),
Expand Down
Loading