diff --git a/imap_processing/glows/l1b/glows_l1b_data.py b/imap_processing/glows/l1b/glows_l1b_data.py index a5409f187..d0c679f16 100644 --- a/imap_processing/glows/l1b/glows_l1b_data.py +++ b/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 diff --git a/imap_processing/quality_flags.py b/imap_processing/quality_flags.py index 84119f650..cb679d73a 100644 --- a/imap_processing/quality_flags.py +++ b/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): diff --git a/imap_processing/tests/glows/conftest.py b/imap_processing/tests/glows/conftest.py index bb6fa285a..9083a8f18 100644 --- a/imap_processing/tests/glows/conftest.py +++ b/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}, diff --git a/imap_processing/tests/glows/test_glows_l1b.py b/imap_processing/tests/glows/test_glows_l1b.py index 3fec7ba0a..4b621f5c8 100644 --- a/imap_processing/tests/glows/test_glows_l1b.py +++ b/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 diff --git a/imap_processing/tests/glows/test_glows_l1b_data.py b/imap_processing/tests/glows/test_glows_l1b_data.py index d7f9114bb..eaaf43a07 100644 --- a/imap_processing/tests/glows/test_glows_l1b_data.py +++ b/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, diff --git a/imap_processing/tests/glows/test_glows_l2.py b/imap_processing/tests/glows/test_glows_l2.py index e73233a77..d4df5fda3 100644 --- a/imap_processing/tests/glows/test_glows_l2.py +++ b/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,