From 216dc8c4f935d801da3663ac430ff1dc57b97a4a Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 23 Jan 2026 08:34:45 -0700 Subject: [PATCH 01/11] Add star sensor code to lo l1b --- imap_processing/lo/l1b/lo_l1b.py | 354 +++++++++++++++++- imap_processing/tests/lo/test_lo_l1b.py | 466 ++++++++++++++++++++++++ 2 files changed, 819 insertions(+), 1 deletion(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index f8cfe4afb3..0497c23cd2 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -22,7 +22,7 @@ frame_transform, lo_instrument_pointing, ) -from imap_processing.spice.repoint import get_pointing_times +from imap_processing.spice.repoint import get_pointing_times, interpolate_repoint_data from imap_processing.spice.spin import get_spin_data, get_spin_number from imap_processing.spice.time import ( et_to_utc, @@ -90,6 +90,25 @@ def lo_l1b( ds = calculate_de_rates(sci_dependencies, anc_dependencies, attr_mgr_l1b) datasets_to_return.append(ds) + # If dependencies are used to create Star Sensor profile + if ( + "imap_lo_l1a_star" in sci_dependencies + and "imap_lo_l1b_nhk" in sci_dependencies + and "imap_lo_l1a_spin" in sci_dependencies + ): + logger.info("\nProcessing IMAP-Lo L1B Star Sensor Profile...") + logical_source = "imap_lo_l1b_star" + + l1a_star = sci_dependencies["imap_lo_l1a_star"] + l1b_nhk = sci_dependencies["imap_lo_l1b_nhk"] + spin_data = sci_dependencies["imap_lo_l1a_spin"] + + l1b_star = initialize_l1b_star( + l1a_star, l1b_nhk, spin_data, attr_mgr_l1b, logical_source + ) + + datasets_to_return.append(l1b_star) + return datasets_to_return @@ -1894,3 +1913,336 @@ def _get_esa_level_indices(epochs: np.ndarray, anc_dependencies: list) -> np.nda energy_step_mapping[esa_idx] = true_esa_step return energy_step_mapping + + +# ============================================================================ +# Star Sensor L1B Processing Functions +# ============================================================================ + + +def filter_valid_star_records( + l1a_star: xr.Dataset, + min_count: int = 700, + time_window_offset: float = 0.0, + time_window_duration: float | None = None, +) -> np.ndarray: + """ + Create boolean mask for valid star sensor records. + + Records are valid if: + 1. COUNT >= min_count (default 700, per algorithm Section 5) + 2. Within specified time window (if provided) + 3. Not during a repoint maneuver + + Parameters + ---------- + l1a_star : xr.Dataset + L1A star sensor dataset containing 'shcoarse' (MET seconds) and 'count'. + min_count : int + Minimum acceptable COUNT value (default: 700). + time_window_offset : float + Time offset in seconds from first record (default: 0.0). + time_window_duration : float | None + Duration of valid time window in seconds (None = no filter, default). + + Returns + ------- + valid_mask : np.ndarray + Boolean array indicating valid records. + """ + # Section 5: Acceptance Criteria - COUNT >= 700 + count_mask = l1a_star["count"].values >= min_count + + # shcoarse is already in MET seconds + shcoarse_sec = l1a_star["shcoarse"].values.astype(np.float64) + + # Section 2.2: Time window filter (if specified) + if time_window_duration is not None: + t0 = shcoarse_sec[0] + time_mask = (shcoarse_sec >= (t0 + time_window_offset)) & ( + shcoarse_sec <= (t0 + time_window_offset + time_window_duration) + ) + valid_mask = count_mask & time_mask + else: + valid_mask = count_mask + + # Filter out repoint maneuvers + repoint_df = interpolate_repoint_data(shcoarse_sec) + # Exclude times where repoint_in_progress is True + repoint_mask = ~repoint_df["repoint_in_progress"].values + valid_mask = valid_mask & repoint_mask + + n_valid = valid_mask.sum() + n_total = len(valid_mask) + logger.info( + f"Star sensor valid records: {n_valid}/{n_total} " + f"({100 * n_valid / n_total:.1f}%)" + ) + + return valid_mask + + +def calculate_star_sensor_profile( + l1a_star: xr.Dataset, + sampling_cadence: float, + spin_period: float, + time_window_offset: float = 0.0, + time_window_duration: float | None = None, + start_angle_offset: float = 62.0, + edge_bins_to_exclude: int = 2, + min_count_threshold: int = 700, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Calculate averaged star sensor amplitude profile vs spin angle. + + Implements the star sensor L1B algorithm. + + Parameters + ---------- + l1a_star : xr.Dataset + L1A star sensor data. + sampling_cadence : float + Sampling period in milliseconds (IFB_DATA_INTERVAL). + spin_period : float + Spin period in seconds. + time_window_offset : float + Time offset for window filtering in seconds (default: 0.0). + time_window_duration : float | None + Duration of time window in seconds (default: None = use all data). + start_angle_offset : float + Starting angle offset in degrees (default: 62.0 = 90° - 28°). + edge_bins_to_exclude : int + Number of edge bins to exclude from each end of the data (default: 2). + min_count_threshold : int + Minimum COUNT value for valid record (default: 700). + + Returns + ------- + spin_angle : np.ndarray + Spin angles in degrees [0-360], shape (720,). + avg_amplitude : np.ndarray + Average amplitude in mV per bin, shape (720,). + count_per_bin : np.ndarray + Number of samples accumulated per bin, shape (720,). + """ + # Section 4, Step 1: Initialize 720-bin sum and count arrays + sum_array = np.zeros(720, dtype=np.float64) + count_array = np.zeros(720, dtype=np.int32) + + # Section 4, Step 2: Get valid record mask + valid_mask = filter_valid_star_records( + l1a_star, min_count_threshold, time_window_offset, time_window_duration + ) + + valid_indices = np.where(valid_mask)[0] + + if len(valid_indices) == 0: + logger.warning( + "No valid star sensor records found. Returning empty profile with FILLVAL." + ) + # Return arrays with FILLVAL for amplitude + spin_angle = np.arange(720) * 0.5 # nominal 0.5 deg bins + avg_amplitude = np.full(720, -1.0e31, dtype=np.float64) + return spin_angle, avg_amplitude, count_array + + # Section 4, Step 3: Accumulate data from valid records (vectorized) + # Get all valid data at once - shape: (n_valid_records, 720) + valid_data = l1a_star["data"].values[valid_indices] + valid_counts = l1a_star["count"].values[valid_indices] + + # Section 4, Step 4: Determine valid bin ranges for each record + # Apply edge exclusion only when count > 2 * edge_bins_to_exclude + use_edge_exclusion = (edge_bins_to_exclude > 0) & ( + valid_counts > 2 * edge_bins_to_exclude + ) + start_bins = np.where(use_edge_exclusion, edge_bins_to_exclude, 0) + end_bins = np.where( + use_edge_exclusion, + np.minimum(valid_counts - edge_bins_to_exclude, 720), + np.minimum(valid_counts, 720), + ) + + # Create mask for valid bins: shape (n_valid_records, 720) + bin_indices = np.arange(720) + valid_bin_mask = (bin_indices[None, :] >= start_bins[:, None]) & ( + bin_indices[None, :] < end_bins[:, None] + ) + + # Apply mask and sum across all valid records + masked_data = np.where(valid_bin_mask, valid_data, 0) + sum_array = masked_data.sum(axis=0).astype(np.float64) + count_array = valid_bin_mask.sum(axis=0).astype(np.int32) + + # Section 4, Step 5: Compute average amplitude per bin + avg_amplitude = np.full(720, -1.0e31, dtype=np.float64) # Initialize with FILLVAL + mask = count_array > 0 + avg_amplitude[mask] = sum_array[mask] / count_array[mask] + + # Section 4, Step 6: Convert bin indices to spin angles + # Section 2.3: DEG_PER_BIN = 360 * (sampling_cadence/1000) / spin_period + deg_per_bin = 360.0 * (sampling_cadence / 1000.0) / spin_period + + # Sample centers at bin center (index + 0.5) * DEG_PER_BIN + bin_indices = np.arange(720) + sample_centers = (bin_indices + 0.5) * deg_per_bin + + # Apply start_angle offset and wrap to [0, 360) + spin_angle = (start_angle_offset + sample_centers) % 360.0 + + logger.info( + f"Star sensor profile calculated: {mask.sum()}/720 bins with valid data" + ) + + return spin_angle, avg_amplitude, count_array + + +def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: + """ + Extract IFB_DATA_INTERVAL from NHK dataset. + + The sampling cadence is already in engineering units after L1B processing. + Formula applied in XTCE: IFB_DATA_INTERVAL = 13.3344 + 0.06945 * DN + + Parameters + ---------- + l1b_nhk : xr.Dataset + L1B NHK dataset with derived values (engineering units). + + Returns + ------- + sampling_cadence : float + Average sampling cadence in milliseconds. + """ + if "IFB_DATA_INTERVAL" not in l1b_nhk: + raise ValueError( + "IFB_DATA_INTERVAL field not found in L1B NHK dataset. " + "Cannot calculate sampling cadence." + ) + + # Get mean value across all epochs (should be relatively constant) + sampling_cadence = float(l1b_nhk["IFB_DATA_INTERVAL"].values.mean()) + + logger.info(f"Sampling cadence from NHK: {sampling_cadence:.3f} ms") + return sampling_cadence + + +def initialize_l1b_star( + l1a_star: xr.Dataset, + l1b_nhk: xr.Dataset, + spin_data: xr.Dataset, + attr_mgr_l1b: ImapCdfAttributes, + logical_source: str, +) -> xr.Dataset: + """ + Initialize and process L1B star sensor dataset. + + Creates an averaged spin profile from L1A star sensor data, computing + the average amplitude per spin angle bin across all valid records. + + Parameters + ---------- + l1a_star : xr.Dataset + The L1A star sensor dataset containing SHCOARSE, COUNT, and DATA fields. + l1b_nhk : xr.Dataset + The L1B NHK dataset containing IFB_DATA_INTERVAL field for sampling cadence. + spin_data : xr.Dataset + The L1A spin dataset used to calculate spin duration. + attr_mgr_l1b : ImapCdfAttributes + Attribute manager for L1B dataset metadata. + logical_source : str + The logical source identifier (e.g., "imap_lo_l1b_star"). + + Returns + ------- + l1b_star : xr.Dataset + L1B star sensor dataset with spin_angle, avg_amplitude, count_per_bin, + and time range metadata. + """ + # Get sampling cadence from NHK + sampling_cadence = get_sampling_cadence_from_nhk(l1b_nhk) + + # Get spin duration from spin data + acq_start, acq_end = convert_start_end_acq_times(spin_data) + avg_spin_durations = get_avg_spin_durations_per_cycle(acq_start, acq_end) + spin_duration = float(avg_spin_durations.mean().values) + logger.info(f"Using spin duration from spin data: {spin_duration:.6f} s") + + # TODO: Read from ancillary config file when available + time_window_offset = 0.0 + time_window_duration = None # None = process all data + start_angle_offset = 62.0 # 90° - 28° + edge_bins_to_exclude = 2 + + # Calculate profile + spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( + l1a_star, + sampling_cadence, + spin_duration, + time_window_offset, + time_window_duration, + start_angle_offset, + edge_bins_to_exclude, + ) + + # Get epoch times from L1A data + start_epoch = l1a_star["epoch"].values[0] + end_epoch = l1a_star["epoch"].values[-1] + epoch_delta = end_epoch - start_epoch + + # Create dataset with global attributes + l1b_star = xr.Dataset( + coords={ + "epoch": xr.DataArray( + [start_epoch], + dims=["epoch"], + attrs=attr_mgr_l1b.get_variable_attributes("epoch"), + ), + "spin_angle_bin": xr.DataArray( + np.arange(720, dtype=np.uint16), + dims=["spin_angle_bin"], + attrs=attr_mgr_l1b.get_variable_attributes("spin_angle_bin"), + ), + }, + attrs=attr_mgr_l1b.get_global_attributes(logical_source), + ) + + # Add data variables + l1b_star["spin_angle"] = xr.DataArray( + spin_angle, + dims=["spin_angle_bin"], + attrs=attr_mgr_l1b.get_variable_attributes("spin_angle"), + ) + + l1b_star["avg_amplitude"] = xr.DataArray( + avg_amplitude, + dims=["spin_angle_bin"], + attrs=attr_mgr_l1b.get_variable_attributes("avg_amplitude"), + ) + + l1b_star["count_per_bin"] = xr.DataArray( + count_per_bin, + dims=["spin_angle_bin"], + attrs=attr_mgr_l1b.get_variable_attributes("count_per_bin"), + ) + + # Add epoch delta (duration in nanoseconds) + l1b_star["epoch_delta"] = xr.DataArray( + [epoch_delta], + dims=["epoch"], + attrs=attr_mgr_l1b.get_variable_attributes("epoch_delta"), + ) + + # Add processing parameters as metadata + l1b_star.attrs["sampling_cadence_ms"] = sampling_cadence + l1b_star.attrs["spin_duration_sec"] = spin_duration + l1b_star.attrs["start_angle_offset_deg"] = start_angle_offset + l1b_star.attrs["edge_bins_excluded"] = edge_bins_to_exclude + l1b_star.attrs["min_count_threshold"] = 700 + l1b_star.attrs["time_window_offset_sec"] = time_window_offset + l1b_star.attrs["time_window_duration_sec"] = ( + "all_data" if time_window_duration is None else time_window_duration + ) + + logger.info("L1B star sensor dataset created successfully") + + return l1b_star diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 5734615ca9..a3c9d2ec6c 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -2,6 +2,7 @@ from unittest.mock import patch import numpy as np +import pandas as pd import pytest import xarray as xr @@ -11,15 +12,19 @@ from imap_processing.lo.l1b.lo_l1b import ( calculate_de_rates, calculate_histogram_rates, + calculate_star_sensor_profile, calculate_tof1_for_golden_triples, convert_start_end_acq_times, convert_tofs_to_eu, create_badtimes_dataset, create_datasets, + filter_valid_star_records, get_avg_spin_durations_per_cycle, + get_sampling_cadence_from_nhk, get_spin_start_times, identify_species, initialize_l1b_de, + initialize_l1b_star, lo_l1b, resweep_histogram_data, set_avg_spin_durations_per_event, @@ -1411,3 +1416,464 @@ def test_calculate_de_rates( output_datasets = lo_l1b(sci_dependencies, anc_dependencies, descriptor="derates") assert len(output_datasets) == 1 assert output_datasets[0].attrs["Logical_source"] == "imap_lo_l1b_derates" + +# ============================================================================ +# Star Sensor L1B Tests +# ============================================================================ +class TestGetSamplingCadenceFromNhk: + """Tests for get_sampling_cadence_from_nhk function.""" + + def test_extracts_mean_cadence(self): + """Test extracting sampling cadence from NHK dataset.""" + # Arrange + l1b_nhk = xr.Dataset( + { + "IFB_DATA_INTERVAL": ("epoch", [20.0, 20.5, 21.0]), + }, + coords={"epoch": [0, 1, 2]}, + ) + expected_cadence = 20.5 # Mean of [20.0, 20.5, 21.0] + + # Act + sampling_cadence = get_sampling_cadence_from_nhk(l1b_nhk) + + # Assert + assert sampling_cadence == expected_cadence + + def test_raises_error_when_field_missing(self): + """Test error when IFB_DATA_INTERVAL field is missing.""" + # Arrange + l1b_nhk = xr.Dataset( + { + "other_field": ("epoch", [1, 2, 3]), + }, + coords={"epoch": [0, 1, 2]}, + ) + + # Act / Assert + with pytest.raises( + ValueError, + match="IFB_DATA_INTERVAL field not found in L1B NHK dataset", + ): + get_sampling_cadence_from_nhk(l1b_nhk) + + +class TestFilterValidStarRecords: + """Tests for filter_valid_star_records function.""" + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_filters_by_count_threshold(self, mock_repoint): + """Test filtering star records by COUNT >= 700.""" + # Arrange - Mock repoint data (no repoints in progress) + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False, False, False]} + ) + + l1a_star = xr.Dataset( + { + "count": ("epoch", [650, 700, 720, 699, 715]), + "shcoarse": ( + "epoch", + np.arange(5, dtype=np.float64), + ), # Already in seconds + }, + coords={"epoch": [0, 1, 2, 3, 4]}, + ) + expected_mask = np.array([False, True, True, False, True]) + + # Act + valid_mask = filter_valid_star_records(l1a_star, min_count=700) + + # Assert + np.testing.assert_array_equal(valid_mask, expected_mask) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_filters_by_count_and_time_window(self, mock_repoint): + """Test filtering star records by both COUNT and time window.""" + # Arrange - Mock repoint data (no repoints in progress) + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False, False, False]} + ) + + # Create times: 0s, 10s, 20s, 30s, 40s (already in seconds) + l1a_star = xr.Dataset( + { + "count": ("epoch", [700, 710, 720, 715, 720]), + "shcoarse": ("epoch", np.array([0, 10, 20, 30, 40], dtype=np.float64)), + }, + coords={"epoch": [0, 1, 2, 3, 4]}, + ) + # Time window: [5s, 25s] - should include epochs 1 and 2 + expected_mask = np.array([False, True, True, False, False]) + + # Act + valid_mask = filter_valid_star_records( + l1a_star, + min_count=700, + time_window_offset=5.0, + time_window_duration=20.0, + ) + + # Assert + np.testing.assert_array_equal(valid_mask, expected_mask) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_processes_all_data_without_time_window(self, mock_repoint): + """Test filtering without time window (process all data).""" + # Arrange - Mock repoint data (no repoints in progress) + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False]} + ) + + l1a_star = xr.Dataset( + { + "count": ("epoch", [700, 710, 720]), + "shcoarse": ("epoch", np.array([0, 10, 20], dtype=np.float64)), + }, + coords={"epoch": [0, 1, 2]}, + ) + expected_mask = np.array([True, True, True]) + + # Act + valid_mask = filter_valid_star_records( + l1a_star, min_count=700, time_window_duration=None + ) + + # Assert + np.testing.assert_array_equal(valid_mask, expected_mask) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_excludes_records_during_repoint(self, mock_repoint): + """Test filtering records during repoint maneuvers.""" + # Arrange - Mock repoint data with some repoints in progress + # Epochs 1 and 3 are during repoint maneuvers + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, True, False, True, False]} + ) + + l1a_star = xr.Dataset( + { + "count": ("epoch", [700, 710, 720, 715, 720]), + "shcoarse": ("epoch", np.arange(5, dtype=np.float64)), + }, + coords={"epoch": [0, 1, 2, 3, 4]}, + ) + # Expected: epochs 0, 2, 4 pass (COUNT >= 700 AND not during repoint) + # Epochs 1 and 3 fail because they are during repoint + expected_mask = np.array([True, False, True, False, True]) + + # Act + valid_mask = filter_valid_star_records(l1a_star, min_count=700) + + # Assert + np.testing.assert_array_equal(valid_mask, expected_mask) + + +class TestCalculateStarSensorProfile: + """Tests for calculate_star_sensor_profile function.""" + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_basic_profile_calculation(self, mock_repoint): + """Test basic star sensor profile calculation.""" + # Arrange + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False]} + ) + # Create simple mock data with 3 records, each with 720 samples + np.random.seed(42) + l1a_star = xr.Dataset( + { + "count": ("epoch", [720, 720, 720]), + "shcoarse": ("epoch", np.array([0.0, 15.0, 30.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.random.randint(100, 200, size=(3, 720), dtype=np.uint16), + ), + }, + coords={"epoch": [0, 1, 2], "samples": np.arange(720)}, + ) + sampling_cadence = 21.0 # ms + spin_duration = 15.0 # seconds + + # Act + spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( + l1a_star, + sampling_cadence, + spin_duration, + time_window_offset=0.0, + time_window_duration=None, + start_angle_offset=62.0, + edge_bins_to_exclude=0, # No edge bins excluded for simplicity + min_count_threshold=700, + ) + + # Assert + assert len(spin_angle) == 720 + assert len(avg_amplitude) == 720 + assert len(count_per_bin) == 720 + # All bins should have 3 samples (3 valid records) + np.testing.assert_array_equal(count_per_bin, np.full(720, 3)) + # Spin angles should be in [0, 360) + assert np.all(spin_angle >= 0) + assert np.all(spin_angle < 360) + # Averages should be reasonable (between 100 and 200 from our mock data) + assert np.all(avg_amplitude >= 100) + assert np.all(avg_amplitude <= 200) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_edge_bins_excluded(self, mock_repoint): + """Test that edge bins are properly excluded.""" + # Arrange + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False]} + ) + l1a_star = xr.Dataset( + { + "count": ("epoch", [720, 720]), + "shcoarse": ("epoch", np.array([0.0, 15.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.ones((2, 720), dtype=np.uint16) * 100, + ), + }, + coords={"epoch": [0, 1], "samples": np.arange(720)}, + ) + sampling_cadence = 21.0 + spin_duration = 15.0 + + # Act + spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( + l1a_star, + sampling_cadence, + spin_duration, + edge_bins_to_exclude=2, # Exclude 2 bins from each end + ) + + # Assert + # First 2 bins and last 2 bins should have count=0 + assert count_per_bin[0] == 0 + assert count_per_bin[1] == 0 + assert count_per_bin[718] == 0 + assert count_per_bin[719] == 0 + # Middle bins should have count=2 (2 valid records) + assert np.all(count_per_bin[2:718] == 2) + # Edge bins should have FILLVAL + assert avg_amplitude[0] == -1.0e31 + assert avg_amplitude[1] == -1.0e31 + assert avg_amplitude[718] == -1.0e31 + assert avg_amplitude[719] == -1.0e31 + # Middle bins should have average value + assert np.all(avg_amplitude[2:718] == 100.0) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_handles_no_valid_records(self, mock_repoint): + """Test handling when no records pass the COUNT threshold.""" + # Arrange + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False]} + ) + l1a_star = xr.Dataset( + { + "count": ("epoch", [650, 600, 699]), # All below 700 + "shcoarse": ("epoch", np.array([0.0, 15.0, 30.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.ones((3, 720), dtype=np.uint16) * 100, + ), + }, + coords={"epoch": [0, 1, 2], "samples": np.arange(720)}, + ) + sampling_cadence = 21.0 + spin_duration = 15.0 + + # Act + spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( + l1a_star, + sampling_cadence, + spin_duration, + min_count_threshold=700, + ) + + # Assert + # All bins should have count=0 + np.testing.assert_array_equal(count_per_bin, np.zeros(720)) + # All averages should be FILLVAL + np.testing.assert_array_equal(avg_amplitude, np.full(720, -1.0e31)) + # Spin angles should still be calculated correctly + assert len(spin_angle) == 720 + assert np.all(spin_angle >= 0) + assert np.all(spin_angle < 360) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_angle_wrapping(self, mock_repoint): + """Test that spin angles wrap correctly to [0, 360) range.""" + # Arrange + mock_repoint.return_value = pd.DataFrame({"repoint_in_progress": [False]}) + l1a_star = xr.Dataset( + { + "count": ("epoch", [720]), + "shcoarse": ("epoch", np.array([0.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.ones((1, 720), dtype=np.uint16) * 100, + ), + }, + coords={"epoch": [0], "samples": np.arange(720)}, + ) + sampling_cadence = 21.0 + spin_duration = 15.0 + start_angle_offset = 350.0 # Large offset to test wrapping + + # Act + spin_angle, _, _ = calculate_star_sensor_profile( + l1a_star, + sampling_cadence, + spin_duration, + start_angle_offset=start_angle_offset, + ) + + # Assert + # All angles should be in [0, 360) + assert np.all(spin_angle >= 0) + assert np.all(spin_angle < 360) + # Check that angles are properly wrapped (not just clamped) + # With offset=350°, first bin should be around 350° + # DEG_PER_BIN = 360 * 0.021 / 15 = 0.504 degrees + # So first bin (index 0.5) should be at 350° + 0.252° = 350.252° + assert 350.0 < spin_angle[0] < 351.0 # First bin near 350° + # Some bins will wrap to the lower range (angles < 100°) + # Check that we have angles both above 300° and below 100° (proof of wrapping) + assert np.any(spin_angle > 300) # Some angles in upper range + assert np.any(spin_angle < 100) # Some angles wrapped to lower range + + +class TestInitializeL1bStar: + """Tests for initialize_l1b_star function.""" + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): + """Test successful initialization of L1B star dataset with spin data.""" + # Arrange + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False]} + ) + np.random.seed(42) + l1a_star = xr.Dataset( + { + "count": ("epoch", [720, 720, 720]), + "shcoarse": ("epoch", np.array([0.0, 15.0, 30.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.random.randint(100, 200, size=(3, 720), dtype=np.uint16), + ), + }, + coords={ + "epoch": met_to_ttj2000ns([0.0, 15.0, 30.0]), + "samples": np.arange(720), + }, + ) + l1b_nhk = xr.Dataset( + { + "IFB_DATA_INTERVAL": ("epoch", [21.0, 21.0, 21.0]), + }, + coords={"epoch": [0, 1, 2]}, + ) + # Create spin data with known spin durations + spin_data = xr.Dataset( + { + "acq_start_sec": ("epoch", [0, 15]), + "acq_start_subsec": ("epoch", [0, 0]), + "acq_end_sec": ("epoch", [420, 435]), # 420s = 28 spins * 15s + "acq_end_subsec": ("epoch", [0, 0]), + "num_completed": ("epoch", [28, 28]), + }, + coords={"epoch": [0, 1]}, + ) + logical_source = "imap_lo_l1b_star" + + # Act + l1b_star = initialize_l1b_star( + l1a_star, l1b_nhk, spin_data, attr_mgr_l1b, logical_source + ) + + # Assert + assert l1b_star.attrs["Logical_source"] == logical_source + assert "epoch" in l1b_star.coords + assert len(l1b_star.coords["epoch"]) == 1 + assert "spin_angle_bin" in l1b_star.coords + assert len(l1b_star.coords["spin_angle_bin"]) == 720 + assert "spin_angle" in l1b_star.data_vars + assert "avg_amplitude" in l1b_star.data_vars + assert "count_per_bin" in l1b_star.data_vars + assert "epoch_delta" in l1b_star.data_vars + # Check attributes + assert "sampling_cadence_ms" in l1b_star.attrs + assert "spin_duration_sec" in l1b_star.attrs + assert l1b_star.attrs["sampling_cadence_ms"] == 21.0 + assert l1b_star.attrs["spin_duration_sec"] == 15.0 + # Check data shapes + assert l1b_star["spin_angle"].shape == (720,) + assert l1b_star["avg_amplitude"].shape == (720,) + assert l1b_star["count_per_bin"].shape == (720,) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): + """Test that L1B star dataset has correct structure and attributes.""" + # Arrange + mock_repoint.return_value = pd.DataFrame({"repoint_in_progress": [False]}) + l1a_star = xr.Dataset( + { + "count": ("epoch", [720]), + "shcoarse": ("epoch", np.array([0.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.ones((1, 720), dtype=np.uint16) * 150, + ), + }, + coords={"epoch": met_to_ttj2000ns([0.0]), "samples": np.arange(720)}, + ) + l1b_nhk = xr.Dataset( + { + "IFB_DATA_INTERVAL": ("epoch", [21.0]), + }, + coords={"epoch": [0]}, + ) + spin_data = xr.Dataset( + { + "acq_start_sec": ("epoch", [0]), + "acq_start_subsec": ("epoch", [0]), + "acq_end_sec": ("epoch", [420]), # 420s = 28 spins * 15s + "acq_end_subsec": ("epoch", [0]), + }, + coords={"epoch": [0]}, + ) + logical_source = "imap_lo_l1b_star" + + # Act + l1b_star = initialize_l1b_star( + l1a_star, l1b_nhk, spin_data, attr_mgr_l1b, logical_source + ) + + # Assert - Check coordinate attributes + assert l1b_star["spin_angle_bin"].attrs["CATDESC"] == "Spin angle bin index" + assert l1b_star["spin_angle_bin"].attrs["VALIDMIN"] == 0 + assert l1b_star["spin_angle_bin"].attrs["VALIDMAX"] == 719 + + # Assert - Check variable attributes + assert l1b_star["spin_angle"].attrs["UNITS"] == "degrees" + assert l1b_star["spin_angle"].attrs["VALIDMIN"] == 0.0 + assert l1b_star["spin_angle"].attrs["VALIDMAX"] == 360.0 + + assert l1b_star["avg_amplitude"].attrs["UNITS"] == "mV" + assert l1b_star["avg_amplitude"].attrs["FILLVAL"] == -1.0e31 + + assert l1b_star["count_per_bin"].attrs["VALIDMIN"] == 0 + assert l1b_star["count_per_bin"].attrs["VALIDMAX"] == 100000 + + # Assert - Check processing parameter attributes + assert "start_angle_offset_deg" in l1b_star.attrs + assert "edge_bins_excluded" in l1b_star.attrs + assert "min_count_threshold" in l1b_star.attrs + assert l1b_star.attrs["start_angle_offset_deg"] == 62.0 + assert l1b_star.attrs["edge_bins_excluded"] == 2 + assert l1b_star.attrs["min_count_threshold"] == 700 From 923cc53968e2a1b7420e31163293a9251af9db27 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 23 Jan 2026 08:39:37 -0700 Subject: [PATCH 02/11] Add star sensor attributes to lo YAMLs --- .../cdf/config/imap_lo_global_cdf_attrs.yaml | 6 +++ .../config/imap_lo_l1b_variable_attrs.yaml | 54 ++++++++++++++++++- 2 files changed, 59 insertions(+), 1 deletion(-) diff --git a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml index 942a317ba9..4527e7c5eb 100644 --- a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml @@ -83,6 +83,12 @@ imap_lo_l1b_prostar: Logical_source: imap_lo_l1b_prostar Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1B Data +imap_lo_l1b_star: + <<: *instrument_base + Data_type: L1B_star>Level-1B Star Sensor Spin-Averaged Profile + Logical_source: imap_lo_l1b_star + Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1B Star Sensor Spin-Averaged Profile Data + imap_lo_l1b_nhk: <<: *instrument_base Data_type: L1B_star>Level-1B Nominal Housekeeping diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index 0e57459581..49d4c95246 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -118,4 +118,56 @@ direction: DEPEND_1: direction_vec FORMAT: I1 LABLAXIS: Direction - LABL_PTR_1: direction_vec_label \ No newline at end of file + LABL_PTR_1: direction_vec_label + +# Star Sensor L1B Attributes +spin_angle_bin: + CATDESC: Spin angle bin index + FIELDNAM: Spin Angle Bin + FORMAT: I3 + VALIDMIN: 0 + VALIDMAX: 719 + VAR_TYPE: support_data + UNITS: ' ' + +spin_angle: + <<: *default + CATDESC: Spin angle in degrees + FIELDNAM: Spin Angle + FORMAT: F8.3 + VALIDMIN: 0.0 + VALIDMAX: 360.0 + UNITS: degrees + LABLAXIS: Spin Angle + DEPEND_0: spin_angle_bin + +avg_amplitude: + <<: *default + CATDESC: Average star sensor amplitude per spin angle bin + FIELDNAM: Average Amplitude + FORMAT: F12.4 + FILLVAL: -1.0000000E+31 + UNITS: mV + LABLAXIS: Amplitude + DEPEND_0: spin_angle_bin + +count_per_bin: + <<: *default + CATDESC: Number of samples per spin angle bin + FIELDNAM: Count Per Bin + FORMAT: I6 + VALIDMIN: 0 + VALIDMAX: 100000 + UNITS: ' ' + LABLAXIS: Sample Count + DEPEND_0: spin_angle_bin + +epoch_delta: + <<: *default + CATDESC: Duration of data range in nanoseconds from start epoch + FIELDNAM: Epoch Delta + FORMAT: I19 + UNITS: ns + VAR_TYPE: support_data + DISPLAY_TYPE: no_plot + TIME_SCALE: Terrestrial Time \ No newline at end of file From 6027d79bed7ae204c6eeee3b7024cb89ea936b86 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 23 Jan 2026 08:56:23 -0700 Subject: [PATCH 03/11] Updates to follow changes in dev after rebase --- imap_processing/lo/l1b/lo_l1b.py | 10 ++-------- imap_processing/tests/lo/test_lo_l1b.py | 2 ++ 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 0497c23cd2..826546c904 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -90,12 +90,7 @@ def lo_l1b( ds = calculate_de_rates(sci_dependencies, anc_dependencies, attr_mgr_l1b) datasets_to_return.append(ds) - # If dependencies are used to create Star Sensor profile - if ( - "imap_lo_l1a_star" in sci_dependencies - and "imap_lo_l1b_nhk" in sci_dependencies - and "imap_lo_l1a_spin" in sci_dependencies - ): + if descriptor == "star": logger.info("\nProcessing IMAP-Lo L1B Star Sensor Profile...") logical_source = "imap_lo_l1b_star" @@ -2162,8 +2157,7 @@ def initialize_l1b_star( sampling_cadence = get_sampling_cadence_from_nhk(l1b_nhk) # Get spin duration from spin data - acq_start, acq_end = convert_start_end_acq_times(spin_data) - avg_spin_durations = get_avg_spin_durations_per_cycle(acq_start, acq_end) + avg_spin_durations = get_avg_spin_durations_per_cycle(spin_data) spin_duration = float(avg_spin_durations.mean().values) logger.info(f"Using spin duration from spin data: {spin_duration:.6f} s") diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index a3c9d2ec6c..13f42c93c0 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1417,6 +1417,7 @@ def test_calculate_de_rates( assert len(output_datasets) == 1 assert output_datasets[0].attrs["Logical_source"] == "imap_lo_l1b_derates" + # ============================================================================ # Star Sensor L1B Tests # ============================================================================ @@ -1844,6 +1845,7 @@ def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): "acq_start_subsec": ("epoch", [0]), "acq_end_sec": ("epoch", [420]), # 420s = 28 spins * 15s "acq_end_subsec": ("epoch", [0]), + "num_completed": ("epoch", [28]), }, coords={"epoch": [0]}, ) From 69ed0d6fa42806d424898b4fcef9052724596473 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 23 Jan 2026 10:14:54 -0700 Subject: [PATCH 04/11] Use spin_angle as coordiante so that data goes from 0->360 --- .../config/imap_lo_l1b_variable_attrs.yaml | 29 ++-- imap_processing/lo/l1b/lo_l1b.py | 117 +++++++------- imap_processing/tests/lo/test_lo_l1b.py | 149 ++++++++++++------ 3 files changed, 170 insertions(+), 125 deletions(-) diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index 49d4c95246..d3369b4061 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -121,25 +121,26 @@ direction: LABL_PTR_1: direction_vec_label # Star Sensor L1B Attributes +spin_angle: + CATDESC: Spin angle in degrees + FIELDNAM: Spin Angle + FORMAT: F8.3 + FILLVAL: -1.0e31 + VALIDMIN: 0.0 + VALIDMAX: 360.0 + VAR_TYPE: support_data + UNITS: deg + spin_angle_bin: - CATDESC: Spin angle bin index + <<: *default + CATDESC: Original spin angle bin index (before sorting) FIELDNAM: Spin Angle Bin FORMAT: I3 VALIDMIN: 0 VALIDMAX: 719 VAR_TYPE: support_data UNITS: ' ' - -spin_angle: - <<: *default - CATDESC: Spin angle in degrees - FIELDNAM: Spin Angle - FORMAT: F8.3 - VALIDMIN: 0.0 - VALIDMAX: 360.0 - UNITS: degrees - LABLAXIS: Spin Angle - DEPEND_0: spin_angle_bin + DEPEND_1: spin_angle avg_amplitude: <<: *default @@ -149,7 +150,7 @@ avg_amplitude: FILLVAL: -1.0000000E+31 UNITS: mV LABLAXIS: Amplitude - DEPEND_0: spin_angle_bin + DEPEND_1: spin_angle count_per_bin: <<: *default @@ -160,7 +161,7 @@ count_per_bin: VALIDMAX: 100000 UNITS: ' ' LABLAXIS: Sample Count - DEPEND_0: spin_angle_bin + DEPEND_1: spin_angle epoch_delta: <<: *default diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 826546c904..75cf1d5e79 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -92,17 +92,8 @@ def lo_l1b( if descriptor == "star": logger.info("\nProcessing IMAP-Lo L1B Star Sensor Profile...") - logical_source = "imap_lo_l1b_star" - - l1a_star = sci_dependencies["imap_lo_l1a_star"] - l1b_nhk = sci_dependencies["imap_lo_l1b_nhk"] - spin_data = sci_dependencies["imap_lo_l1a_spin"] - - l1b_star = initialize_l1b_star( - l1a_star, l1b_nhk, spin_data, attr_mgr_l1b, logical_source - ) - - datasets_to_return.append(l1b_star) + ds = l1b_star(sci_dependencies, attr_mgr_l1b) + datasets_to_return.append(ds) return datasets_to_return @@ -1997,7 +1988,7 @@ def calculate_star_sensor_profile( l1a_star : xr.Dataset L1A star sensor data. sampling_cadence : float - Sampling period in milliseconds (IFB_DATA_INTERVAL). + Sampling period in milliseconds (ifb_data_interval). spin_period : float Spin period in seconds. time_window_offset : float @@ -2093,10 +2084,10 @@ def calculate_star_sensor_profile( def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: """ - Extract IFB_DATA_INTERVAL from NHK dataset. + Extract ifb_data_interval from NHK dataset. The sampling cadence is already in engineering units after L1B processing. - Formula applied in XTCE: IFB_DATA_INTERVAL = 13.3344 + 0.06945 * DN + Formula applied in XTCE: ifb_data_interval = 13.3344 + 0.06945 * DN Parameters ---------- @@ -2108,51 +2099,47 @@ def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: sampling_cadence : float Average sampling cadence in milliseconds. """ - if "IFB_DATA_INTERVAL" not in l1b_nhk: - raise ValueError( - "IFB_DATA_INTERVAL field not found in L1B NHK dataset. " + if "ifb_data_interval" not in l1b_nhk: + raise KeyError( + "ifb_data_interval field not found in L1B NHK dataset. " "Cannot calculate sampling cadence." ) # Get mean value across all epochs (should be relatively constant) - sampling_cadence = float(l1b_nhk["IFB_DATA_INTERVAL"].values.mean()) + sampling_cadence = float(l1b_nhk["ifb_data_interval"].values.mean()) logger.info(f"Sampling cadence from NHK: {sampling_cadence:.3f} ms") return sampling_cadence -def initialize_l1b_star( - l1a_star: xr.Dataset, - l1b_nhk: xr.Dataset, - spin_data: xr.Dataset, +def l1b_star( + sci_dependencies: dict, attr_mgr_l1b: ImapCdfAttributes, - logical_source: str, ) -> xr.Dataset: """ - Initialize and process L1B star sensor dataset. + Create the IMAP-Lo L1B Star Sensor dataset. Creates an averaged spin profile from L1A star sensor data, computing the average amplitude per spin angle bin across all valid records. Parameters ---------- - l1a_star : xr.Dataset - The L1A star sensor dataset containing SHCOARSE, COUNT, and DATA fields. - l1b_nhk : xr.Dataset - The L1B NHK dataset containing IFB_DATA_INTERVAL field for sampling cadence. - spin_data : xr.Dataset - The L1A spin dataset used to calculate spin duration. + sci_dependencies : dict + Dictionary of datasets needed for L1B data product creation in xarray Datasets. attr_mgr_l1b : ImapCdfAttributes Attribute manager for L1B dataset metadata. - logical_source : str - The logical source identifier (e.g., "imap_lo_l1b_star"). Returns ------- - l1b_star : xr.Dataset + l1b_star_ds : xr.Dataset L1B star sensor dataset with spin_angle, avg_amplitude, count_per_bin, and time range metadata. """ + logical_source = "imap_lo_l1b_star" + l1a_star = sci_dependencies["imap_lo_l1a_star"] + l1b_nhk = sci_dependencies["imap_lo_l1b_nhk"] + spin_data = sci_dependencies["imap_lo_l1a_spin"] + # Get sampling cadence from NHK sampling_cadence = get_sampling_cadence_from_nhk(l1b_nhk) @@ -2178,65 +2165,77 @@ def initialize_l1b_star( edge_bins_to_exclude, ) + # Sort data so spin_angle is monotonically increasing from 0 to 360 + # Use argsort to get indices that would sort spin_angle + sort_indices = np.argsort(spin_angle) + spin_angle_sorted = spin_angle[sort_indices] + avg_amplitude_sorted = avg_amplitude[sort_indices] + count_per_bin_sorted = count_per_bin[sort_indices] + # Original bin indices, reordered to match the sorted spin_angle + original_bin_indices = sort_indices.astype(np.uint16) + # Get epoch times from L1A data start_epoch = l1a_star["epoch"].values[0] end_epoch = l1a_star["epoch"].values[-1] epoch_delta = end_epoch - start_epoch - # Create dataset with global attributes - l1b_star = xr.Dataset( + # Create dataset with spin_angle as the coordinate + l1b_star_ds = xr.Dataset( coords={ "epoch": xr.DataArray( [start_epoch], dims=["epoch"], attrs=attr_mgr_l1b.get_variable_attributes("epoch"), ), - "spin_angle_bin": xr.DataArray( - np.arange(720, dtype=np.uint16), - dims=["spin_angle_bin"], - attrs=attr_mgr_l1b.get_variable_attributes("spin_angle_bin"), + "spin_angle": xr.DataArray( + spin_angle_sorted, + dims=["spin_angle"], + attrs=attr_mgr_l1b.get_variable_attributes( + "spin_angle", check_schema=False + ), ), }, attrs=attr_mgr_l1b.get_global_attributes(logical_source), ) - # Add data variables - l1b_star["spin_angle"] = xr.DataArray( - spin_angle, - dims=["spin_angle_bin"], - attrs=attr_mgr_l1b.get_variable_attributes("spin_angle"), + # Add spin_angle_bin as a variable (original bin indices) + # All variables must have epoch as first dimension for SPDF CDF compliance + l1b_star_ds["spin_angle_bin"] = xr.DataArray( + original_bin_indices[np.newaxis, :], + dims=["epoch", "spin_angle"], + attrs=attr_mgr_l1b.get_variable_attributes("spin_angle_bin"), ) - l1b_star["avg_amplitude"] = xr.DataArray( - avg_amplitude, - dims=["spin_angle_bin"], + l1b_star_ds["avg_amplitude"] = xr.DataArray( + avg_amplitude_sorted[np.newaxis, :], + dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("avg_amplitude"), ) - l1b_star["count_per_bin"] = xr.DataArray( - count_per_bin, - dims=["spin_angle_bin"], + l1b_star_ds["count_per_bin"] = xr.DataArray( + count_per_bin_sorted[np.newaxis, :], + dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("count_per_bin"), ) # Add epoch delta (duration in nanoseconds) - l1b_star["epoch_delta"] = xr.DataArray( + l1b_star_ds["epoch_delta"] = xr.DataArray( [epoch_delta], dims=["epoch"], attrs=attr_mgr_l1b.get_variable_attributes("epoch_delta"), ) # Add processing parameters as metadata - l1b_star.attrs["sampling_cadence_ms"] = sampling_cadence - l1b_star.attrs["spin_duration_sec"] = spin_duration - l1b_star.attrs["start_angle_offset_deg"] = start_angle_offset - l1b_star.attrs["edge_bins_excluded"] = edge_bins_to_exclude - l1b_star.attrs["min_count_threshold"] = 700 - l1b_star.attrs["time_window_offset_sec"] = time_window_offset - l1b_star.attrs["time_window_duration_sec"] = ( + l1b_star_ds.attrs["sampling_cadence_ms"] = sampling_cadence + l1b_star_ds.attrs["spin_duration_sec"] = spin_duration + l1b_star_ds.attrs["start_angle_offset_deg"] = start_angle_offset + l1b_star_ds.attrs["edge_bins_excluded"] = edge_bins_to_exclude + l1b_star_ds.attrs["min_count_threshold"] = 700 + l1b_star_ds.attrs["time_window_offset_sec"] = time_window_offset + l1b_star_ds.attrs["time_window_duration_sec"] = ( "all_data" if time_window_duration is None else time_window_duration ) logger.info("L1B star sensor dataset created successfully") - return l1b_star + return l1b_star_ds diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 13f42c93c0..d350b3b0a0 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1,4 +1,5 @@ from collections import namedtuple +from pathlib import Path from unittest.mock import patch import numpy as np @@ -8,7 +9,7 @@ from imap_processing import imap_module_directory from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes -from imap_processing.cdf.utils import load_cdf +from imap_processing.cdf.utils import load_cdf, write_cdf from imap_processing.lo.l1b.lo_l1b import ( calculate_de_rates, calculate_histogram_rates, @@ -24,7 +25,7 @@ get_spin_start_times, identify_species, initialize_l1b_de, - initialize_l1b_star, + l1b_star, lo_l1b, resweep_histogram_data, set_avg_spin_durations_per_event, @@ -1429,7 +1430,7 @@ def test_extracts_mean_cadence(self): # Arrange l1b_nhk = xr.Dataset( { - "IFB_DATA_INTERVAL": ("epoch", [20.0, 20.5, 21.0]), + "ifb_data_interval": ("epoch", [20.0, 20.5, 21.0]), }, coords={"epoch": [0, 1, 2]}, ) @@ -1442,7 +1443,7 @@ def test_extracts_mean_cadence(self): assert sampling_cadence == expected_cadence def test_raises_error_when_field_missing(self): - """Test error when IFB_DATA_INTERVAL field is missing.""" + """Test error when ifb_data_interval field is missing.""" # Arrange l1b_nhk = xr.Dataset( { @@ -1454,7 +1455,7 @@ def test_raises_error_when_field_missing(self): # Act / Assert with pytest.raises( ValueError, - match="IFB_DATA_INTERVAL field not found in L1B NHK dataset", + match="ifb_data_interval field not found in L1B NHK dataset", ): get_sampling_cadence_from_nhk(l1b_nhk) @@ -1748,8 +1749,8 @@ def test_angle_wrapping(self, mock_repoint): assert np.any(spin_angle < 100) # Some angles wrapped to lower range -class TestInitializeL1bStar: - """Tests for initialize_l1b_star function.""" +class TestL1bStar: + """Tests for l1b_star function.""" @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): @@ -1775,7 +1776,7 @@ def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): ) l1b_nhk = xr.Dataset( { - "IFB_DATA_INTERVAL": ("epoch", [21.0, 21.0, 21.0]), + "ifb_data_interval": ("epoch", [21.0, 21.0, 21.0]), }, coords={"epoch": [0, 1, 2]}, ) @@ -1790,32 +1791,43 @@ def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): }, coords={"epoch": [0, 1]}, ) - logical_source = "imap_lo_l1b_star" + sci_dependencies = { + "imap_lo_l1a_star": l1a_star, + "imap_lo_l1b_nhk": l1b_nhk, + "imap_lo_l1a_spin": spin_data, + } # Act - l1b_star = initialize_l1b_star( - l1a_star, l1b_nhk, spin_data, attr_mgr_l1b, logical_source - ) + l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) # Assert - assert l1b_star.attrs["Logical_source"] == logical_source - assert "epoch" in l1b_star.coords - assert len(l1b_star.coords["epoch"]) == 1 - assert "spin_angle_bin" in l1b_star.coords - assert len(l1b_star.coords["spin_angle_bin"]) == 720 - assert "spin_angle" in l1b_star.data_vars - assert "avg_amplitude" in l1b_star.data_vars - assert "count_per_bin" in l1b_star.data_vars - assert "epoch_delta" in l1b_star.data_vars + assert l1b_star_ds.attrs["Logical_source"] == "imap_lo_l1b_star" + assert "epoch" in l1b_star_ds.coords + assert len(l1b_star_ds.coords["epoch"]) == 1 + # spin_angle is now the coordinate (monotonically increasing) + assert "spin_angle" in l1b_star_ds.coords + assert len(l1b_star_ds.coords["spin_angle"]) == 720 + # spin_angle_bin is now a data variable + assert "spin_angle_bin" in l1b_star_ds.data_vars + assert "avg_amplitude" in l1b_star_ds.data_vars + assert "count_per_bin" in l1b_star_ds.data_vars + assert "epoch_delta" in l1b_star_ds.data_vars + # Check that spin_angle is monotonically increasing + spin_angles = l1b_star_ds.coords["spin_angle"].values + assert np.all(np.diff(spin_angles) > 0), ( + "spin_angle should be monotonically increasing" + ) + assert spin_angles[0] >= 0.0 + assert spin_angles[-1] < 360.0 # Check attributes - assert "sampling_cadence_ms" in l1b_star.attrs - assert "spin_duration_sec" in l1b_star.attrs - assert l1b_star.attrs["sampling_cadence_ms"] == 21.0 - assert l1b_star.attrs["spin_duration_sec"] == 15.0 - # Check data shapes - assert l1b_star["spin_angle"].shape == (720,) - assert l1b_star["avg_amplitude"].shape == (720,) - assert l1b_star["count_per_bin"].shape == (720,) + assert "sampling_cadence_ms" in l1b_star_ds.attrs + assert "spin_duration_sec" in l1b_star_ds.attrs + assert l1b_star_ds.attrs["sampling_cadence_ms"] == 21.0 + assert l1b_star_ds.attrs["spin_duration_sec"] == 15.0 + # Check data shapes - all variables have epoch as first dimension + assert l1b_star_ds["spin_angle_bin"].shape == (1, 720) + assert l1b_star_ds["avg_amplitude"].shape == (1, 720) + assert l1b_star_ds["count_per_bin"].shape == (1, 720) @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): @@ -1835,7 +1847,7 @@ def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): ) l1b_nhk = xr.Dataset( { - "IFB_DATA_INTERVAL": ("epoch", [21.0]), + "ifb_data_interval": ("epoch", [21.0]), }, coords={"epoch": [0]}, ) @@ -1849,33 +1861,66 @@ def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): }, coords={"epoch": [0]}, ) - logical_source = "imap_lo_l1b_star" + sci_dependencies = { + "imap_lo_l1a_star": l1a_star, + "imap_lo_l1b_nhk": l1b_nhk, + "imap_lo_l1a_spin": spin_data, + } # Act - l1b_star = initialize_l1b_star( - l1a_star, l1b_nhk, spin_data, attr_mgr_l1b, logical_source - ) + l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) - # Assert - Check coordinate attributes - assert l1b_star["spin_angle_bin"].attrs["CATDESC"] == "Spin angle bin index" - assert l1b_star["spin_angle_bin"].attrs["VALIDMIN"] == 0 - assert l1b_star["spin_angle_bin"].attrs["VALIDMAX"] == 719 + # Assert - Check spin_angle coordinate attributes + assert l1b_star_ds.coords["spin_angle"].attrs["UNITS"] == "degrees" + assert l1b_star_ds.coords["spin_angle"].attrs["VALIDMIN"] == 0.0 + assert l1b_star_ds.coords["spin_angle"].attrs["VALIDMAX"] == 360.0 - # Assert - Check variable attributes - assert l1b_star["spin_angle"].attrs["UNITS"] == "degrees" - assert l1b_star["spin_angle"].attrs["VALIDMIN"] == 0.0 - assert l1b_star["spin_angle"].attrs["VALIDMAX"] == 360.0 + # Assert - Check spin_angle_bin variable attributes (now a data variable) + assert ( + "Original spin angle bin index" + in l1b_star_ds["spin_angle_bin"].attrs["CATDESC"] + ) + assert l1b_star_ds["spin_angle_bin"].attrs["VALIDMIN"] == 0 + assert l1b_star_ds["spin_angle_bin"].attrs["VALIDMAX"] == 719 - assert l1b_star["avg_amplitude"].attrs["UNITS"] == "mV" - assert l1b_star["avg_amplitude"].attrs["FILLVAL"] == -1.0e31 + assert l1b_star_ds["avg_amplitude"].attrs["UNITS"] == "mV" + assert l1b_star_ds["avg_amplitude"].attrs["FILLVAL"] == -1.0e31 - assert l1b_star["count_per_bin"].attrs["VALIDMIN"] == 0 - assert l1b_star["count_per_bin"].attrs["VALIDMAX"] == 100000 + assert l1b_star_ds["count_per_bin"].attrs["VALIDMIN"] == 0 + assert l1b_star_ds["count_per_bin"].attrs["VALIDMAX"] == 100000 # Assert - Check processing parameter attributes - assert "start_angle_offset_deg" in l1b_star.attrs - assert "edge_bins_excluded" in l1b_star.attrs - assert "min_count_threshold" in l1b_star.attrs - assert l1b_star.attrs["start_angle_offset_deg"] == 62.0 - assert l1b_star.attrs["edge_bins_excluded"] == 2 - assert l1b_star.attrs["min_count_threshold"] == 700 + assert "start_angle_offset_deg" in l1b_star_ds.attrs + assert "edge_bins_excluded" in l1b_star_ds.attrs + assert "min_count_threshold" in l1b_star_ds.attrs + assert l1b_star_ds.attrs["start_angle_offset_deg"] == 62.0 + assert l1b_star_ds.attrs["edge_bins_excluded"] == 2 + assert l1b_star_ds.attrs["min_count_threshold"] == 700 + + +def test_star_integration(use_test_repoint_data_csv): + """Temporary integration test for star data.""" + use_test_repoint_data_csv( + Path( + "/Users/plummert/Projects/imap/data/prod/imap/spice/repoint/imap_2026_022_01.repoint" + ) + ) + star_path = Path( + "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2025/11/imap_lo_l1a_star_20251110-repoint00044_v001.cdf" + ) + spin_path = Path( + "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2025/11/imap_lo_l1a_spin_20251110-repoint00044_v001.cdf" + ) + nhk_path = Path( + "/Users/plummert/Projects/imap/data/prod/imap/lo/l1b/2025/11/imap_lo_l1b_nhk_20251110-repoint00044_v001.cdf" + ) + sci_dependencies = { + "imap_lo_l1a_star": load_cdf(star_path), + "imap_lo_l1a_spin": load_cdf(spin_path), + "imap_lo_l1b_nhk": load_cdf(nhk_path), + } + anc_dependencies = [] + descriptor = "star" + result = lo_l1b(sci_dependencies, anc_dependencies, descriptor) + assert len(result) == 1 + print(write_cdf(result[0])) From 2a86469896b6ce0eb9b182db734e0baffdfa5b0b Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Fri, 23 Jan 2026 14:21:11 -0700 Subject: [PATCH 05/11] Move epoch to fractional DOY to time module --- .../config/imap_lo_l1b_variable_attrs.yaml | 26 ++++++- imap_processing/lo/l1b/lo_l1b.py | 17 +++++ imap_processing/spice/time.py | 30 ++++++++ imap_processing/tests/lo/test_lo_l1b.py | 71 +++++++++++++++++++ 4 files changed, 143 insertions(+), 1 deletion(-) diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index d3369b4061..ce795592be 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -171,4 +171,28 @@ epoch_delta: UNITS: ns VAR_TYPE: support_data DISPLAY_TYPE: no_plot - TIME_SCALE: Terrestrial Time \ No newline at end of file + TIME_SCALE: Terrestrial Time + +start_doy: + <<: *default + CATDESC: Fractional day of year for start of data range + FIELDNAM: Start Day of Year + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 1.0 + VALIDMAX: 367.0 + UNITS: day + VAR_TYPE: support_data + LABLAXIS: Start DOY + +end_doy: + <<: *default + CATDESC: Fractional day of year for end of data range + FIELDNAM: End Day of Year + FORMAT: F12.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 1.0 + VALIDMAX: 367.0 + UNITS: day + VAR_TYPE: support_data + LABLAXIS: End DOY \ No newline at end of file diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 75cf1d5e79..e54ceb1d09 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -25,6 +25,7 @@ from imap_processing.spice.repoint import get_pointing_times, interpolate_repoint_data from imap_processing.spice.spin import get_spin_data, get_spin_number from imap_processing.spice.time import ( + epoch_to_fractional_doy, et_to_utc, met_to_ttj2000ns, ttj2000ns_to_et, @@ -2225,6 +2226,22 @@ def l1b_star( attrs=attr_mgr_l1b.get_variable_attributes("epoch_delta"), ) + # Add start and end day of year as floating point values + start_doy = epoch_to_fractional_doy(start_epoch) + end_doy = epoch_to_fractional_doy(start_epoch + epoch_delta) + + l1b_star_ds["start_doy"] = xr.DataArray( + [start_doy], + dims=["epoch"], + attrs=attr_mgr_l1b.get_variable_attributes("start_doy"), + ) + + l1b_star_ds["end_doy"] = xr.DataArray( + [end_doy], + dims=["epoch"], + attrs=attr_mgr_l1b.get_variable_attributes("end_doy"), + ) + # Add processing parameters as metadata l1b_star_ds.attrs["sampling_cadence_ms"] = sampling_cadence l1b_star_ds.attrs["spin_duration_sec"] = spin_duration diff --git a/imap_processing/spice/time.py b/imap_processing/spice/time.py index 8129cb2fe9..1b91972232 100644 --- a/imap_processing/spice/time.py +++ b/imap_processing/spice/time.py @@ -407,3 +407,33 @@ def epoch_to_doy(epoch: np.ndarray) -> npt.NDArray: return np.array( [datetime.fromisoformat(date).timetuple().tm_yday for date in time_strings] ) + + +def epoch_to_fractional_doy(epoch_ttj2000ns: int) -> float: + """ + Convert epoch in TTJ2000ns to floating point day of year. + + Parameters + ---------- + epoch_ttj2000ns : int + Epoch in TTJ2000ns format (nanoseconds since J2000). + + Returns + ------- + doy : float + Floating point day of year (1.0 = Jan 1 00:00:00). + """ + # Convert to ephemeris time, then to UTC string + et = ttj2000ns_to_et(epoch_ttj2000ns) + utc_str = et_to_utc(et) # Returns ISO format: "YYYY-MM-DDTHH:MM:SS.sss" + + # Parse the datetime (remove trailing 'Z' if present) + dt = datetime.fromisoformat(utc_str.rstrip("Z")) + + # Calculate day of year as floating point + # Day of year starts at 1, so Jan 1 00:00:00 = 1.0 + start_of_year = datetime(dt.year, 1, 1) + delta = dt - start_of_year + doy = 1.0 + delta.total_seconds() / 86400.0 + + return doy diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index d350b3b0a0..76daf4c70e 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1897,6 +1897,77 @@ def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): assert l1b_star_ds.attrs["edge_bins_excluded"] == 2 assert l1b_star_ds.attrs["min_count_threshold"] == 700 + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_start_and_end_doy_variables(self, mock_repoint, attr_mgr_l1b): + """Test that start_doy and end_doy variables are computed correctly.""" + # Arrange + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False, False, False]} + ) + np.random.seed(42) + # Create epochs spanning 30 seconds + l1a_star = xr.Dataset( + { + "count": ("epoch", [720, 720, 720]), + "shcoarse": ("epoch", np.array([0.0, 15.0, 30.0], dtype=np.float64)), + "data": ( + ("epoch", "samples"), + np.random.randint(100, 200, size=(3, 720), dtype=np.uint16), + ), + }, + coords={ + "epoch": met_to_ttj2000ns([0.0, 15.0, 30.0]), + "samples": np.arange(720), + }, + ) + l1b_nhk = xr.Dataset( + { + "ifb_data_interval": ("epoch", [21.0, 21.0, 21.0]), + }, + coords={"epoch": [0, 1, 2]}, + ) + spin_data = xr.Dataset( + { + "acq_start_sec": ("epoch", [0, 15]), + "acq_start_subsec": ("epoch", [0, 0]), + "acq_end_sec": ("epoch", [420, 435]), + "acq_end_subsec": ("epoch", [0, 0]), + "num_completed": ("epoch", [28, 28]), + }, + coords={"epoch": [0, 1]}, + ) + sci_dependencies = { + "imap_lo_l1a_star": l1a_star, + "imap_lo_l1b_nhk": l1b_nhk, + "imap_lo_l1a_spin": spin_data, + } + + # Act + l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) + + # Assert - Check that start_doy and end_doy exist + assert "start_doy" in l1b_star_ds.data_vars + assert "end_doy" in l1b_star_ds.data_vars + + # Assert - Check dimensions + assert l1b_star_ds["start_doy"].dims == ("epoch",) + assert l1b_star_ds["end_doy"].dims == ("epoch",) + + # Assert - Check values are valid day of year (1.0 to 366.x for leap years) + start_doy = l1b_star_ds["start_doy"].values[0] + end_doy = l1b_star_ds["end_doy"].values[0] + assert 1.0 <= start_doy <= 367.0 + assert 1.0 <= end_doy <= 367.0 + + # Assert - end_doy should be >= start_doy (data spans 30 seconds) + assert end_doy >= start_doy + + # Assert - Check attributes + assert l1b_star_ds["start_doy"].attrs["UNITS"] == "day" + assert l1b_star_ds["end_doy"].attrs["UNITS"] == "day" + assert "Fractional day of year" in l1b_star_ds["start_doy"].attrs["CATDESC"] + assert "Fractional day of year" in l1b_star_ds["end_doy"].attrs["CATDESC"] + def test_star_integration(use_test_repoint_data_csv): """Temporary integration test for star data.""" From 8778dfbfd2ec97f360df27e188ed8bc60dab310c Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 26 Jan 2026 10:23:02 -0700 Subject: [PATCH 06/11] Average star sensor data every 64 spin set --- .../config/imap_lo_l1b_variable_attrs.yaml | 24 +- imap_processing/lo/l1b/lo_l1b.py | 270 ++++++++------ imap_processing/tests/lo/test_lo_l1b.py | 337 +++++++++++------- 3 files changed, 395 insertions(+), 236 deletions(-) diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index ce795592be..3e01676182 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -163,16 +163,6 @@ count_per_bin: LABLAXIS: Sample Count DEPEND_1: spin_angle -epoch_delta: - <<: *default - CATDESC: Duration of data range in nanoseconds from start epoch - FIELDNAM: Epoch Delta - FORMAT: I19 - UNITS: ns - VAR_TYPE: support_data - DISPLAY_TYPE: no_plot - TIME_SCALE: Terrestrial Time - start_doy: <<: *default CATDESC: Fractional day of year for start of data range @@ -195,4 +185,16 @@ end_doy: VALIDMAX: 367.0 UNITS: day VAR_TYPE: support_data - LABLAXIS: End DOY \ No newline at end of file + LABLAXIS: End DOY + +pointing_mid_met: + <<: *default + CATDESC: Mission Elapsed Time at the midpoint of the pointing + FIELDNAM: Pointing Mid MET + FORMAT: F15.6 + FILLVAL: -1.0000000E+31 + VALIDMIN: 0.0 + VALIDMAX: 1.0E+10 + UNITS: s + VAR_TYPE: support_data + LABLAXIS: Pointing Mid MET \ No newline at end of file diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index e54ceb1d09..b0d3ef32cd 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -22,7 +22,11 @@ frame_transform, lo_instrument_pointing, ) -from imap_processing.spice.repoint import get_pointing_times, interpolate_repoint_data +from imap_processing.spice.repoint import ( + get_pointing_mid_time, + get_pointing_times, + interpolate_repoint_data, +) from imap_processing.spice.spin import get_spin_data, get_spin_number from imap_processing.spice.time import ( epoch_to_fractional_doy, @@ -1969,20 +1973,77 @@ def filter_valid_star_records( return valid_mask -def calculate_star_sensor_profile( +def calculate_star_sensor_profile_for_group( + data: np.ndarray, + counts: np.ndarray, + edge_bins_to_exclude: int = 2, +) -> tuple[np.ndarray, np.ndarray]: + """ + Calculate averaged star sensor amplitude profile for a group of records. + + Parameters + ---------- + data : np.ndarray + Star sensor data array, shape (n_records, 720). + counts : np.ndarray + Count values for each record, shape (n_records,). + edge_bins_to_exclude : int + Number of edge bins to exclude from each end of the data (default: 2). + + Returns + ------- + avg_amplitude : np.ndarray + Average amplitude in mV per bin, shape (720,). + count_per_bin : np.ndarray + Number of samples accumulated per bin, shape (720,). + """ + if len(data) == 0: + return np.full(720, -1.0e31, dtype=np.float64), np.zeros(720, dtype=np.int32) + + # Determine valid bin ranges for each record + use_edge_exclusion = (edge_bins_to_exclude > 0) & ( + counts > 2 * edge_bins_to_exclude + ) + start_bins = np.where(use_edge_exclusion, edge_bins_to_exclude, 0) + end_bins = np.where( + use_edge_exclusion, + np.minimum(counts - edge_bins_to_exclude, 720), + np.minimum(counts, 720), + ) + + # Create mask for valid bins: shape (n_records, 720) + bin_indices = np.arange(720) + valid_bin_mask = (bin_indices[None, :] >= start_bins[:, None]) & ( + bin_indices[None, :] < end_bins[:, None] + ) + + # Apply mask and sum across all records + masked_data = np.where(valid_bin_mask, data, 0) + sum_array = masked_data.sum(axis=0).astype(np.float64) + count_array = valid_bin_mask.sum(axis=0).astype(np.int32) + + # Compute average amplitude per bin + avg_amplitude = np.full(720, -1.0e31, dtype=np.float64) # Initialize with FILLVAL + mask = count_array > 0 + avg_amplitude[mask] = sum_array[mask] / count_array[mask] + + return avg_amplitude, count_array + + +def calculate_star_sensor_profiles_by_group( l1a_star: xr.Dataset, sampling_cadence: float, spin_period: float, - time_window_offset: float = 0.0, - time_window_duration: float | None = None, + group_size: int = 64, start_angle_offset: float = 62.0, edge_bins_to_exclude: int = 2, min_count_threshold: int = 700, -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ - Calculate averaged star sensor amplitude profile vs spin angle. + Calculate averaged star sensor amplitude profiles for groups of records. - Implements the star sensor L1B algorithm. + Groups L1A star sensor records into chunks of `group_size` and calculates + an averaged profile for each group. Parameters ---------- @@ -1992,10 +2053,8 @@ def calculate_star_sensor_profile( Sampling period in milliseconds (ifb_data_interval). spin_period : float Spin period in seconds. - time_window_offset : float - Time offset for window filtering in seconds (default: 0.0). - time_window_duration : float | None - Duration of time window in seconds (default: None = use all data). + group_size : int + Number of records per group (default: 64). start_angle_offset : float Starting angle offset in degrees (default: 62.0 = 90° - 28°). edge_bins_to_exclude : int @@ -2007,80 +2066,73 @@ def calculate_star_sensor_profile( ------- spin_angle : np.ndarray Spin angles in degrees [0-360], shape (720,). - avg_amplitude : np.ndarray - Average amplitude in mV per bin, shape (720,). - count_per_bin : np.ndarray - Number of samples accumulated per bin, shape (720,). - """ - # Section 4, Step 1: Initialize 720-bin sum and count arrays - sum_array = np.zeros(720, dtype=np.float64) - count_array = np.zeros(720, dtype=np.int32) - - # Section 4, Step 2: Get valid record mask + group_epochs : np.ndarray + Start epoch for each group, shape (n_groups,). + avg_amplitudes : np.ndarray + Average amplitude in mV per bin per group, shape (n_groups, 720). + counts_per_bin : np.ndarray + Number of samples accumulated per bin per group, shape (n_groups, 720). + """ + # Get valid record mask valid_mask = filter_valid_star_records( - l1a_star, min_count_threshold, time_window_offset, time_window_duration + l1a_star, min_count_threshold, time_window_offset=0.0, time_window_duration=None ) valid_indices = np.where(valid_mask)[0] + n_valid = len(valid_indices) + + # Calculate spin angles (same for all groups) + deg_per_bin = 360.0 * (sampling_cadence / 1000.0) / spin_period + bin_indices = np.arange(720) + sample_centers = (bin_indices + 0.5) * deg_per_bin + spin_angle = (start_angle_offset + sample_centers) % 360.0 - if len(valid_indices) == 0: + if n_valid == 0: logger.warning( "No valid star sensor records found. Returning empty profile with FILLVAL." ) - # Return arrays with FILLVAL for amplitude - spin_angle = np.arange(720) * 0.5 # nominal 0.5 deg bins - avg_amplitude = np.full(720, -1.0e31, dtype=np.float64) - return spin_angle, avg_amplitude, count_array + return ( + spin_angle, + np.array([], dtype=np.int64), + np.empty((0, 720), dtype=np.float64), + np.empty((0, 720), dtype=np.int32), + ) - # Section 4, Step 3: Accumulate data from valid records (vectorized) - # Get all valid data at once - shape: (n_valid_records, 720) + # Get valid data valid_data = l1a_star["data"].values[valid_indices] valid_counts = l1a_star["count"].values[valid_indices] + valid_epochs = l1a_star["epoch"].values[valid_indices] - # Section 4, Step 4: Determine valid bin ranges for each record - # Apply edge exclusion only when count > 2 * edge_bins_to_exclude - use_edge_exclusion = (edge_bins_to_exclude > 0) & ( - valid_counts > 2 * edge_bins_to_exclude - ) - start_bins = np.where(use_edge_exclusion, edge_bins_to_exclude, 0) - end_bins = np.where( - use_edge_exclusion, - np.minimum(valid_counts - edge_bins_to_exclude, 720), - np.minimum(valid_counts, 720), - ) + # Calculate number of groups (include partial groups) + n_groups = (n_valid + group_size - 1) // group_size - # Create mask for valid bins: shape (n_valid_records, 720) - bin_indices = np.arange(720) - valid_bin_mask = (bin_indices[None, :] >= start_bins[:, None]) & ( - bin_indices[None, :] < end_bins[:, None] + logger.info( + f"Processing {n_valid} valid records into {n_groups} groups of {group_size}" ) - # Apply mask and sum across all valid records - masked_data = np.where(valid_bin_mask, valid_data, 0) - sum_array = masked_data.sum(axis=0).astype(np.float64) - count_array = valid_bin_mask.sum(axis=0).astype(np.int32) + # Initialize output arrays + avg_amplitudes = np.zeros((n_groups, 720), dtype=np.float64) + counts_per_bin = np.zeros((n_groups, 720), dtype=np.int32) + group_epochs = np.zeros(n_groups, dtype=np.int64) - # Section 4, Step 5: Compute average amplitude per bin - avg_amplitude = np.full(720, -1.0e31, dtype=np.float64) # Initialize with FILLVAL - mask = count_array > 0 - avg_amplitude[mask] = sum_array[mask] / count_array[mask] + # Process each group + for group_idx in range(n_groups): + start_idx = group_idx * group_size + end_idx = min(start_idx + group_size, n_valid) - # Section 4, Step 6: Convert bin indices to spin angles - # Section 2.3: DEG_PER_BIN = 360 * (sampling_cadence/1000) / spin_period - deg_per_bin = 360.0 * (sampling_cadence / 1000.0) / spin_period + group_data = valid_data[start_idx:end_idx] + group_counts = valid_counts[start_idx:end_idx] - # Sample centers at bin center (index + 0.5) * DEG_PER_BIN - bin_indices = np.arange(720) - sample_centers = (bin_indices + 0.5) * deg_per_bin - - # Apply start_angle offset and wrap to [0, 360) - spin_angle = (start_angle_offset + sample_centers) % 360.0 + # Calculate profile for this group + avg_amp, count_arr = calculate_star_sensor_profile_for_group( + group_data, group_counts, edge_bins_to_exclude + ) - logger.info( - f"Star sensor profile calculated: {mask.sum()}/720 bins with valid data" - ) + avg_amplitudes[group_idx] = avg_amp + counts_per_bin[group_idx] = count_arr + group_epochs[group_idx] = valid_epochs[start_idx] - return spin_angle, avg_amplitude, count_array + return spin_angle, group_epochs, avg_amplitudes, counts_per_bin def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: @@ -2116,12 +2168,14 @@ def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: def l1b_star( sci_dependencies: dict, attr_mgr_l1b: ImapCdfAttributes, + group_size: int = 64, ) -> xr.Dataset: """ Create the IMAP-Lo L1B Star Sensor dataset. - Creates an averaged spin profile from L1A star sensor data, computing - the average amplitude per spin angle bin across all valid records. + Creates averaged spin profiles from L1A star sensor data, computing + the average amplitude per spin angle bin for each group of records. + Each group contains `group_size` consecutive valid records. Parameters ---------- @@ -2129,12 +2183,14 @@ def l1b_star( Dictionary of datasets needed for L1B data product creation in xarray Datasets. attr_mgr_l1b : ImapCdfAttributes Attribute manager for L1B dataset metadata. + group_size : int + Number of records to average per group (default: 64). Returns ------- l1b_star_ds : xr.Dataset L1B star sensor dataset with spin_angle, avg_amplitude, count_per_bin, - and time range metadata. + and time range metadata. Each epoch corresponds to a group of records. """ logical_source = "imap_lo_l1b_star" l1a_star = sci_dependencies["imap_lo_l1a_star"] @@ -2150,41 +2206,42 @@ def l1b_star( logger.info(f"Using spin duration from spin data: {spin_duration:.6f} s") # TODO: Read from ancillary config file when available - time_window_offset = 0.0 - time_window_duration = None # None = process all data start_angle_offset = 62.0 # 90° - 28° edge_bins_to_exclude = 2 - # Calculate profile - spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( + # Calculate profiles for each 64-spin group + ( + spin_angle, + group_epochs, + avg_amplitudes, + counts_per_bin, + ) = calculate_star_sensor_profiles_by_group( l1a_star, sampling_cadence, spin_duration, - time_window_offset, - time_window_duration, - start_angle_offset, - edge_bins_to_exclude, + group_size=group_size, + start_angle_offset=start_angle_offset, + edge_bins_to_exclude=edge_bins_to_exclude, ) # Sort data so spin_angle is monotonically increasing from 0 to 360 - # Use argsort to get indices that would sort spin_angle sort_indices = np.argsort(spin_angle) spin_angle_sorted = spin_angle[sort_indices] - avg_amplitude_sorted = avg_amplitude[sort_indices] - count_per_bin_sorted = count_per_bin[sort_indices] + # Apply sorting to all groups' data + avg_amplitudes_sorted = avg_amplitudes[:, sort_indices] + counts_per_bin_sorted = counts_per_bin[:, sort_indices] # Original bin indices, reordered to match the sorted spin_angle original_bin_indices = sort_indices.astype(np.uint16) - # Get epoch times from L1A data - start_epoch = l1a_star["epoch"].values[0] - end_epoch = l1a_star["epoch"].values[-1] - epoch_delta = end_epoch - start_epoch + # Get global epoch times from L1A data for start_doy and end_doy + global_start_epoch = l1a_star["epoch"].values[0] + global_end_epoch = l1a_star["epoch"].values[-1] - # Create dataset with spin_angle as the coordinate + # Create dataset with spin_angle as coordinate and multiple epochs l1b_star_ds = xr.Dataset( coords={ "epoch": xr.DataArray( - [start_epoch], + group_epochs, dims=["epoch"], attrs=attr_mgr_l1b.get_variable_attributes("epoch"), ), @@ -2200,45 +2257,47 @@ def l1b_star( ) # Add spin_angle_bin as a variable (original bin indices) - # All variables must have epoch as first dimension for SPDF CDF compliance + # Broadcast to all epochs since bin mapping is the same for all groups l1b_star_ds["spin_angle_bin"] = xr.DataArray( - original_bin_indices[np.newaxis, :], + np.broadcast_to(original_bin_indices, (len(group_epochs), 720)), dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("spin_angle_bin"), ) l1b_star_ds["avg_amplitude"] = xr.DataArray( - avg_amplitude_sorted[np.newaxis, :], + avg_amplitudes_sorted, dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("avg_amplitude"), ) l1b_star_ds["count_per_bin"] = xr.DataArray( - count_per_bin_sorted[np.newaxis, :], + counts_per_bin_sorted, dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("count_per_bin"), ) - # Add epoch delta (duration in nanoseconds) - l1b_star_ds["epoch_delta"] = xr.DataArray( - [epoch_delta], - dims=["epoch"], - attrs=attr_mgr_l1b.get_variable_attributes("epoch_delta"), + # Add pointing mid time (MET) as a scalar value + # Use the first epoch to determine which pointing we're in + first_met = ttj2000ns_to_met(global_start_epoch) + pointing_mid_met = get_pointing_mid_time(first_met) + l1b_star_ds["pointing_mid_met"] = xr.DataArray( + pointing_mid_met, + attrs=attr_mgr_l1b.get_variable_attributes( + "pointing_mid_met", check_schema=False + ), ) - # Add start and end day of year as floating point values - start_doy = epoch_to_fractional_doy(start_epoch) - end_doy = epoch_to_fractional_doy(start_epoch + epoch_delta) + # Add global start and end day of year as scalar values + start_doy = epoch_to_fractional_doy(global_start_epoch) + end_doy = epoch_to_fractional_doy(global_end_epoch) l1b_star_ds["start_doy"] = xr.DataArray( - [start_doy], - dims=["epoch"], + start_doy, attrs=attr_mgr_l1b.get_variable_attributes("start_doy"), ) l1b_star_ds["end_doy"] = xr.DataArray( - [end_doy], - dims=["epoch"], + end_doy, attrs=attr_mgr_l1b.get_variable_attributes("end_doy"), ) @@ -2248,11 +2307,10 @@ def l1b_star( l1b_star_ds.attrs["start_angle_offset_deg"] = start_angle_offset l1b_star_ds.attrs["edge_bins_excluded"] = edge_bins_to_exclude l1b_star_ds.attrs["min_count_threshold"] = 700 - l1b_star_ds.attrs["time_window_offset_sec"] = time_window_offset - l1b_star_ds.attrs["time_window_duration_sec"] = ( - "all_data" if time_window_duration is None else time_window_duration - ) + l1b_star_ds.attrs["group_size"] = group_size - logger.info("L1B star sensor dataset created successfully") + logger.info( + f"L1B star sensor dataset created successfully with {len(group_epochs)} groups" + ) return l1b_star_ds diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 76daf4c70e..57d04c8c2c 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -13,7 +13,8 @@ from imap_processing.lo.l1b.lo_l1b import ( calculate_de_rates, calculate_histogram_rates, - calculate_star_sensor_profile, + calculate_star_sensor_profile_for_group, + calculate_star_sensor_profiles_by_group, calculate_tof1_for_golden_triples, convert_start_end_acq_times, convert_tofs_to_eu, @@ -1572,83 +1573,38 @@ def test_excludes_records_during_repoint(self, mock_repoint): class TestCalculateStarSensorProfile: - """Tests for calculate_star_sensor_profile function.""" + """Tests for star sensor profile calculation functions.""" - @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_basic_profile_calculation(self, mock_repoint): - """Test basic star sensor profile calculation.""" - # Arrange - mock_repoint.return_value = pd.DataFrame( - {"repoint_in_progress": [False, False, False]} - ) - # Create simple mock data with 3 records, each with 720 samples + def test_profile_for_group_basic(self): + """Test basic star sensor profile calculation for a group.""" + # Arrange - 3 records with uniform data np.random.seed(42) - l1a_star = xr.Dataset( - { - "count": ("epoch", [720, 720, 720]), - "shcoarse": ("epoch", np.array([0.0, 15.0, 30.0], dtype=np.float64)), - "data": ( - ("epoch", "samples"), - np.random.randint(100, 200, size=(3, 720), dtype=np.uint16), - ), - }, - coords={"epoch": [0, 1, 2], "samples": np.arange(720)}, - ) - sampling_cadence = 21.0 # ms - spin_duration = 15.0 # seconds + data = np.random.randint(100, 200, size=(3, 720)).astype(np.uint16) + counts = np.array([720, 720, 720]) # Act - spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( - l1a_star, - sampling_cadence, - spin_duration, - time_window_offset=0.0, - time_window_duration=None, - start_angle_offset=62.0, - edge_bins_to_exclude=0, # No edge bins excluded for simplicity - min_count_threshold=700, + avg_amplitude, count_per_bin = calculate_star_sensor_profile_for_group( + data, counts, edge_bins_to_exclude=0 ) # Assert - assert len(spin_angle) == 720 assert len(avg_amplitude) == 720 assert len(count_per_bin) == 720 - # All bins should have 3 samples (3 valid records) + # All bins should have 3 samples np.testing.assert_array_equal(count_per_bin, np.full(720, 3)) - # Spin angles should be in [0, 360) - assert np.all(spin_angle >= 0) - assert np.all(spin_angle < 360) - # Averages should be reasonable (between 100 and 200 from our mock data) + # Averages should be between 100 and 200 assert np.all(avg_amplitude >= 100) assert np.all(avg_amplitude <= 200) - @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_edge_bins_excluded(self, mock_repoint): + def test_profile_for_group_edge_bins_excluded(self): """Test that edge bins are properly excluded.""" - # Arrange - mock_repoint.return_value = pd.DataFrame( - {"repoint_in_progress": [False, False]} - ) - l1a_star = xr.Dataset( - { - "count": ("epoch", [720, 720]), - "shcoarse": ("epoch", np.array([0.0, 15.0], dtype=np.float64)), - "data": ( - ("epoch", "samples"), - np.ones((2, 720), dtype=np.uint16) * 100, - ), - }, - coords={"epoch": [0, 1], "samples": np.arange(720)}, - ) - sampling_cadence = 21.0 - spin_duration = 15.0 + # Arrange - 2 records with uniform data + data = np.ones((2, 720), dtype=np.uint16) * 100 + counts = np.array([720, 720]) # Act - spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( - l1a_star, - sampling_cadence, - spin_duration, - edge_bins_to_exclude=2, # Exclude 2 bins from each end + avg_amplitude, count_per_bin = calculate_star_sensor_profile_for_group( + data, counts, edge_bins_to_exclude=2 ) # Assert @@ -1657,7 +1613,7 @@ def test_edge_bins_excluded(self, mock_repoint): assert count_per_bin[1] == 0 assert count_per_bin[718] == 0 assert count_per_bin[719] == 0 - # Middle bins should have count=2 (2 valid records) + # Middle bins should have count=2 assert np.all(count_per_bin[2:718] == 2) # Edge bins should have FILLVAL assert avg_amplitude[0] == -1.0e31 @@ -1667,8 +1623,72 @@ def test_edge_bins_excluded(self, mock_repoint): # Middle bins should have average value assert np.all(avg_amplitude[2:718] == 100.0) + def test_profile_for_group_empty_data(self): + """Test handling of empty data array.""" + # Arrange + data = np.empty((0, 720), dtype=np.uint16) + counts = np.array([], dtype=np.int32) + + # Act + avg_amplitude, count_per_bin = calculate_star_sensor_profile_for_group( + data, counts + ) + + # Assert + np.testing.assert_array_equal(count_per_bin, np.zeros(720)) + np.testing.assert_array_equal(avg_amplitude, np.full(720, -1.0e31)) + + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_profiles_by_group_creates_correct_groups(self, mock_repoint): + """Test that profiles are grouped correctly into 64-record groups.""" + # Arrange - Create 150 records (should produce 3 groups: 64, 64, 22) + n_records = 150 + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False] * n_records} + ) + l1a_star = xr.Dataset( + { + "count": ("epoch", [720] * n_records), + "shcoarse": ( + "epoch", + np.arange(n_records, dtype=np.float64) * 15.0, + ), + "data": ( + ("epoch", "samples"), + np.ones((n_records, 720), dtype=np.uint16) * 100, + ), + }, + coords={ + "epoch": met_to_ttj2000ns(np.arange(n_records) * 15.0), + "samples": np.arange(720), + }, + ) + + # Act + ( + spin_angle, + group_epochs, + avg_amplitudes, + counts_per_bin, + ) = calculate_star_sensor_profiles_by_group( + l1a_star, + sampling_cadence=21.0, + spin_period=15.0, + group_size=64, + ) + + # Assert + assert len(spin_angle) == 720 + assert len(group_epochs) == 3 # 150 records -> 3 groups + assert avg_amplitudes.shape == (3, 720) + assert counts_per_bin.shape == (3, 720) + # First two groups should have 64 samples per bin, last group 22 + assert np.all(counts_per_bin[0, 2:718] == 64) + assert np.all(counts_per_bin[1, 2:718] == 64) + assert np.all(counts_per_bin[2, 2:718] == 22) + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_handles_no_valid_records(self, mock_repoint): + def test_profiles_by_group_handles_no_valid_records(self, mock_repoint): """Test handling when no records pass the COUNT threshold.""" # Arrange mock_repoint.return_value = pd.DataFrame( @@ -1683,31 +1703,32 @@ def test_handles_no_valid_records(self, mock_repoint): np.ones((3, 720), dtype=np.uint16) * 100, ), }, - coords={"epoch": [0, 1, 2], "samples": np.arange(720)}, + coords={ + "epoch": met_to_ttj2000ns([0.0, 15.0, 30.0]), + "samples": np.arange(720), + }, ) - sampling_cadence = 21.0 - spin_duration = 15.0 # Act - spin_angle, avg_amplitude, count_per_bin = calculate_star_sensor_profile( + ( + spin_angle, + group_epochs, + avg_amplitudes, + counts_per_bin, + ) = calculate_star_sensor_profiles_by_group( l1a_star, - sampling_cadence, - spin_duration, + sampling_cadence=21.0, + spin_period=15.0, min_count_threshold=700, ) # Assert - # All bins should have count=0 - np.testing.assert_array_equal(count_per_bin, np.zeros(720)) - # All averages should be FILLVAL - np.testing.assert_array_equal(avg_amplitude, np.full(720, -1.0e31)) - # Spin angles should still be calculated correctly assert len(spin_angle) == 720 - assert np.all(spin_angle >= 0) - assert np.all(spin_angle < 360) + assert len(group_epochs) == 0 # No valid records + assert avg_amplitudes.shape == (0, 720) @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_angle_wrapping(self, mock_repoint): + def test_profiles_by_group_angle_wrapping(self, mock_repoint): """Test that spin angles wrap correctly to [0, 360) range.""" # Arrange mock_repoint.return_value = pd.DataFrame({"repoint_in_progress": [False]}) @@ -1720,65 +1741,63 @@ def test_angle_wrapping(self, mock_repoint): np.ones((1, 720), dtype=np.uint16) * 100, ), }, - coords={"epoch": [0], "samples": np.arange(720)}, + coords={"epoch": met_to_ttj2000ns([0.0]), "samples": np.arange(720)}, ) - sampling_cadence = 21.0 - spin_duration = 15.0 - start_angle_offset = 350.0 # Large offset to test wrapping # Act - spin_angle, _, _ = calculate_star_sensor_profile( + spin_angle, _, _, _ = calculate_star_sensor_profiles_by_group( l1a_star, - sampling_cadence, - spin_duration, - start_angle_offset=start_angle_offset, + sampling_cadence=21.0, + spin_period=15.0, + start_angle_offset=350.0, # Large offset to test wrapping ) # Assert - # All angles should be in [0, 360) assert np.all(spin_angle >= 0) assert np.all(spin_angle < 360) - # Check that angles are properly wrapped (not just clamped) # With offset=350°, first bin should be around 350° - # DEG_PER_BIN = 360 * 0.021 / 15 = 0.504 degrees - # So first bin (index 0.5) should be at 350° + 0.252° = 350.252° - assert 350.0 < spin_angle[0] < 351.0 # First bin near 350° - # Some bins will wrap to the lower range (angles < 100°) - # Check that we have angles both above 300° and below 100° (proof of wrapping) - assert np.any(spin_angle > 300) # Some angles in upper range + assert 350.0 < spin_angle[0] < 351.0 + # Some bins will wrap to the lower range + assert np.any(spin_angle > 300) assert np.any(spin_angle < 100) # Some angles wrapped to lower range class TestL1bStar: """Tests for l1b_star function.""" + @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): + def test_initializes_with_spin_data( + self, mock_repoint, mock_pointing_mid, attr_mgr_l1b + ): """Test successful initialization of L1B star dataset with spin data.""" - # Arrange + # Arrange - Create 150 records to produce multiple groups + n_records = 150 mock_repoint.return_value = pd.DataFrame( - {"repoint_in_progress": [False, False, False]} + {"repoint_in_progress": [False] * n_records} ) + mock_pointing_mid.return_value = 1000.0 # Mock pointing mid time in MET np.random.seed(42) + met_times = np.arange(n_records, dtype=np.float64) * 15.0 l1a_star = xr.Dataset( { - "count": ("epoch", [720, 720, 720]), - "shcoarse": ("epoch", np.array([0.0, 15.0, 30.0], dtype=np.float64)), + "count": ("epoch", [720] * n_records), + "shcoarse": ("epoch", met_times), "data": ( ("epoch", "samples"), - np.random.randint(100, 200, size=(3, 720), dtype=np.uint16), + np.random.randint(100, 200, size=(n_records, 720), dtype=np.uint16), ), }, coords={ - "epoch": met_to_ttj2000ns([0.0, 15.0, 30.0]), + "epoch": met_to_ttj2000ns(met_times), "samples": np.arange(720), }, ) l1b_nhk = xr.Dataset( { - "ifb_data_interval": ("epoch", [21.0, 21.0, 21.0]), + "ifb_data_interval": ("epoch", [21.0] * n_records), }, - coords={"epoch": [0, 1, 2]}, + coords={"epoch": list(range(n_records))}, ) # Create spin data with known spin durations spin_data = xr.Dataset( @@ -1798,12 +1817,13 @@ def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): } # Act - l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) + l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b, group_size=64) # Assert assert l1b_star_ds.attrs["Logical_source"] == "imap_lo_l1b_star" assert "epoch" in l1b_star_ds.coords - assert len(l1b_star_ds.coords["epoch"]) == 1 + # 150 records / 64 group_size = 3 groups (64 + 64 + 22) + assert len(l1b_star_ds.coords["epoch"]) == 3 # spin_angle is now the coordinate (monotonically increasing) assert "spin_angle" in l1b_star_ds.coords assert len(l1b_star_ds.coords["spin_angle"]) == 720 @@ -1811,7 +1831,7 @@ def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): assert "spin_angle_bin" in l1b_star_ds.data_vars assert "avg_amplitude" in l1b_star_ds.data_vars assert "count_per_bin" in l1b_star_ds.data_vars - assert "epoch_delta" in l1b_star_ds.data_vars + assert "pointing_mid_met" in l1b_star_ds.data_vars # Check that spin_angle is monotonically increasing spin_angles = l1b_star_ds.coords["spin_angle"].values assert np.all(np.diff(spin_angles) > 0), ( @@ -1822,18 +1842,27 @@ def test_initializes_with_spin_data(self, mock_repoint, attr_mgr_l1b): # Check attributes assert "sampling_cadence_ms" in l1b_star_ds.attrs assert "spin_duration_sec" in l1b_star_ds.attrs + assert "group_size" in l1b_star_ds.attrs assert l1b_star_ds.attrs["sampling_cadence_ms"] == 21.0 assert l1b_star_ds.attrs["spin_duration_sec"] == 15.0 + assert l1b_star_ds.attrs["group_size"] == 64 # Check data shapes - all variables have epoch as first dimension - assert l1b_star_ds["spin_angle_bin"].shape == (1, 720) - assert l1b_star_ds["avg_amplitude"].shape == (1, 720) - assert l1b_star_ds["count_per_bin"].shape == (1, 720) - + assert l1b_star_ds["spin_angle_bin"].shape == (3, 720) + assert l1b_star_ds["avg_amplitude"].shape == (3, 720) + assert l1b_star_ds["count_per_bin"].shape == (3, 720) + # Check pointing_mid_met is a scalar with expected value + assert l1b_star_ds["pointing_mid_met"].dims == () + assert float(l1b_star_ds["pointing_mid_met"].values) == 1000.0 + + @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): + def test_dataset_structure_and_attributes( + self, mock_repoint, mock_pointing_mid, attr_mgr_l1b + ): """Test that L1B star dataset has correct structure and attributes.""" # Arrange mock_repoint.return_value = pd.DataFrame({"repoint_in_progress": [False]}) + mock_pointing_mid.return_value = 1000.0 l1a_star = xr.Dataset( { "count": ("epoch", [720]), @@ -1871,7 +1900,7 @@ def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) # Assert - Check spin_angle coordinate attributes - assert l1b_star_ds.coords["spin_angle"].attrs["UNITS"] == "degrees" + assert l1b_star_ds.coords["spin_angle"].attrs["UNITS"] == "deg" assert l1b_star_ds.coords["spin_angle"].attrs["VALIDMIN"] == 0.0 assert l1b_star_ds.coords["spin_angle"].attrs["VALIDMAX"] == 360.0 @@ -1897,10 +1926,14 @@ def test_dataset_structure_and_attributes(self, mock_repoint, attr_mgr_l1b): assert l1b_star_ds.attrs["edge_bins_excluded"] == 2 assert l1b_star_ds.attrs["min_count_threshold"] == 700 + @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") - def test_start_and_end_doy_variables(self, mock_repoint, attr_mgr_l1b): + def test_start_and_end_doy_variables( + self, mock_repoint, mock_pointing_mid, attr_mgr_l1b + ): """Test that start_doy and end_doy variables are computed correctly.""" # Arrange + mock_pointing_mid.return_value = 1000.0 # Mock pointing mid time in MET mock_repoint.return_value = pd.DataFrame( {"repoint_in_progress": [False, False, False]} ) @@ -1945,17 +1978,17 @@ def test_start_and_end_doy_variables(self, mock_repoint, attr_mgr_l1b): # Act l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) - # Assert - Check that start_doy and end_doy exist + # Assert - Check that start_doy and end_doy exist as scalars (global values) assert "start_doy" in l1b_star_ds.data_vars assert "end_doy" in l1b_star_ds.data_vars - # Assert - Check dimensions - assert l1b_star_ds["start_doy"].dims == ("epoch",) - assert l1b_star_ds["end_doy"].dims == ("epoch",) + # Assert - Check they are scalar values (no dimensions) + assert l1b_star_ds["start_doy"].dims == () + assert l1b_star_ds["end_doy"].dims == () # Assert - Check values are valid day of year (1.0 to 366.x for leap years) - start_doy = l1b_star_ds["start_doy"].values[0] - end_doy = l1b_star_ds["end_doy"].values[0] + start_doy = float(l1b_star_ds["start_doy"].values) + end_doy = float(l1b_star_ds["end_doy"].values) assert 1.0 <= start_doy <= 367.0 assert 1.0 <= end_doy <= 367.0 @@ -1968,6 +2001,72 @@ def test_start_and_end_doy_variables(self, mock_repoint, attr_mgr_l1b): assert "Fractional day of year" in l1b_star_ds["start_doy"].attrs["CATDESC"] assert "Fractional day of year" in l1b_star_ds["end_doy"].attrs["CATDESC"] + @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") + @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") + def test_multiple_groups_created( + self, mock_repoint, mock_pointing_mid, attr_mgr_l1b + ): + """Test that multiple 64-spin groups are created correctly.""" + # Arrange - Create 150 records to produce 3 groups (64 + 64 + 22) + n_records = 150 + mock_pointing_mid.return_value = 1000.0 # Mock pointing mid time in MET + mock_repoint.return_value = pd.DataFrame( + {"repoint_in_progress": [False] * n_records} + ) + met_times = np.arange(n_records, dtype=np.float64) * 15.0 + l1a_star = xr.Dataset( + { + "count": ("epoch", [720] * n_records), + "shcoarse": ("epoch", met_times), + "data": ( + ("epoch", "samples"), + np.ones((n_records, 720), dtype=np.uint16) * 100, + ), + }, + coords={ + "epoch": met_to_ttj2000ns(met_times), + "samples": np.arange(720), + }, + ) + l1b_nhk = xr.Dataset( + { + "ifb_data_interval": ("epoch", [21.0] * n_records), + }, + coords={"epoch": list(range(n_records))}, + ) + spin_data = xr.Dataset( + { + "acq_start_sec": ("epoch", [0]), + "acq_start_subsec": ("epoch", [0]), + "acq_end_sec": ("epoch", [420]), + "acq_end_subsec": ("epoch", [0]), + "num_completed": ("epoch", [28]), + }, + coords={"epoch": [0]}, + ) + sci_dependencies = { + "imap_lo_l1a_star": l1a_star, + "imap_lo_l1b_nhk": l1b_nhk, + "imap_lo_l1a_spin": spin_data, + } + + # Act + l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b, group_size=64) + + # Assert + assert len(l1b_star_ds.coords["epoch"]) == 3 + # Check pointing_mid_met is present (scalar value) + assert "pointing_mid_met" in l1b_star_ds.data_vars + assert l1b_star_ds["pointing_mid_met"].dims == () + # First group epoch should be the first L1A epoch + assert l1b_star_ds.coords["epoch"].values[0] == met_to_ttj2000ns([0.0])[0] + # Second group epoch should be record 64 + assert l1b_star_ds.coords["epoch"].values[1] == met_to_ttj2000ns([64 * 15.0])[0] + # Third group epoch should be record 128 + assert ( + l1b_star_ds.coords["epoch"].values[2] == met_to_ttj2000ns([128 * 15.0])[0] + ) + def test_star_integration(use_test_repoint_data_csv): """Temporary integration test for star data.""" From 5f8e28cc9efe56f4816ae2be3ff30891a68f423c Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 26 Jan 2026 13:47:12 -0700 Subject: [PATCH 07/11] Leverage xarray instead of using numpy sorting and selection --- .../config/imap_lo_l1b_variable_attrs.yaml | 45 ++----- imap_processing/lo/l1b/lo_l1b.py | 121 ++++++++---------- imap_processing/tests/lo/test_lo_l1b.py | 67 ++++------ 3 files changed, 87 insertions(+), 146 deletions(-) diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index 3e01676182..cb673f9d54 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -121,6 +121,13 @@ direction: LABL_PTR_1: direction_vec_label # Star Sensor L1B Attributes +met: + <<: *default + CATDESC: Start Mission Elapsed Time of averaging period + FIELDNAM: MET + FORMAT: I16 + VAR_TYPE: support_data + spin_angle: CATDESC: Spin angle in degrees FIELDNAM: Spin Angle @@ -140,7 +147,7 @@ spin_angle_bin: VALIDMAX: 719 VAR_TYPE: support_data UNITS: ' ' - DEPEND_1: spin_angle + DEPEND_0: spin_angle avg_amplitude: <<: *default @@ -162,39 +169,3 @@ count_per_bin: UNITS: ' ' LABLAXIS: Sample Count DEPEND_1: spin_angle - -start_doy: - <<: *default - CATDESC: Fractional day of year for start of data range - FIELDNAM: Start Day of Year - FORMAT: F12.6 - FILLVAL: -1.0000000E+31 - VALIDMIN: 1.0 - VALIDMAX: 367.0 - UNITS: day - VAR_TYPE: support_data - LABLAXIS: Start DOY - -end_doy: - <<: *default - CATDESC: Fractional day of year for end of data range - FIELDNAM: End Day of Year - FORMAT: F12.6 - FILLVAL: -1.0000000E+31 - VALIDMIN: 1.0 - VALIDMAX: 367.0 - UNITS: day - VAR_TYPE: support_data - LABLAXIS: End DOY - -pointing_mid_met: - <<: *default - CATDESC: Mission Elapsed Time at the midpoint of the pointing - FIELDNAM: Pointing Mid MET - FORMAT: F15.6 - FILLVAL: -1.0000000E+31 - VALIDMIN: 0.0 - VALIDMAX: 1.0E+10 - UNITS: s - VAR_TYPE: support_data - LABLAXIS: Pointing Mid MET \ No newline at end of file diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index b0d3ef32cd..ab82d6c92e 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -21,6 +21,8 @@ cartesian_to_latitudinal, frame_transform, lo_instrument_pointing, + get_spacecraft_to_instrument_spin_phase_offset, + instrument_pointing, ) from imap_processing.spice.repoint import ( get_pointing_mid_time, @@ -1976,7 +1978,7 @@ def filter_valid_star_records( def calculate_star_sensor_profile_for_group( data: np.ndarray, counts: np.ndarray, - edge_bins_to_exclude: int = 2, + end_bins_to_exclude: int = 2, ) -> tuple[np.ndarray, np.ndarray]: """ Calculate averaged star sensor amplitude profile for a group of records. @@ -1987,8 +1989,8 @@ def calculate_star_sensor_profile_for_group( Star sensor data array, shape (n_records, 720). counts : np.ndarray Count values for each record, shape (n_records,). - edge_bins_to_exclude : int - Number of edge bins to exclude from each end of the data (default: 2). + end_bins_to_exclude : int + Number of bins to exclude from end of each row of data (default: 2). Returns ------- @@ -2001,21 +2003,16 @@ def calculate_star_sensor_profile_for_group( return np.full(720, -1.0e31, dtype=np.float64), np.zeros(720, dtype=np.int32) # Determine valid bin ranges for each record - use_edge_exclusion = (edge_bins_to_exclude > 0) & ( - counts > 2 * edge_bins_to_exclude - ) - start_bins = np.where(use_edge_exclusion, edge_bins_to_exclude, 0) + use_edge_exclusion = (end_bins_to_exclude > 0) & (counts > end_bins_to_exclude) end_bins = np.where( use_edge_exclusion, - np.minimum(counts - edge_bins_to_exclude, 720), + np.minimum(counts - end_bins_to_exclude, 720), np.minimum(counts, 720), ) # Create mask for valid bins: shape (n_records, 720) bin_indices = np.arange(720) - valid_bin_mask = (bin_indices[None, :] >= start_bins[:, None]) & ( - bin_indices[None, :] < end_bins[:, None] - ) + valid_bin_mask = bin_indices[None, :] < end_bins[:, None] # Apply mask and sum across all records masked_data = np.where(valid_bin_mask, data, 0) @@ -2023,7 +2020,7 @@ def calculate_star_sensor_profile_for_group( count_array = valid_bin_mask.sum(axis=0).astype(np.int32) # Compute average amplitude per bin - avg_amplitude = np.full(720, -1.0e31, dtype=np.float64) # Initialize with FILLVAL + avg_amplitude = np.full(720, np.nan, dtype=np.float64) mask = count_array > 0 avg_amplitude[mask] = sum_array[mask] / count_array[mask] @@ -2036,7 +2033,7 @@ def calculate_star_sensor_profiles_by_group( spin_period: float, group_size: int = 64, start_angle_offset: float = 62.0, - edge_bins_to_exclude: int = 2, + end_bins_to_exclude: int = 2, min_count_threshold: int = 700, ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ @@ -2057,8 +2054,8 @@ def calculate_star_sensor_profiles_by_group( Number of records per group (default: 64). start_angle_offset : float Starting angle offset in degrees (default: 62.0 = 90° - 28°). - edge_bins_to_exclude : int - Number of edge bins to exclude from each end of the data (default: 2). + end_bins_to_exclude : int + Number of ending bins to exclude from each average (default: 2). min_count_threshold : int Minimum COUNT value for valid record (default: 700). @@ -2066,8 +2063,8 @@ def calculate_star_sensor_profiles_by_group( ------- spin_angle : np.ndarray Spin angles in degrees [0-360], shape (720,). - group_epochs : np.ndarray - Start epoch for each group, shape (n_groups,). + group_mets : np.ndarray + Start MET for each group, shape (n_groups,). avg_amplitudes : np.ndarray Average amplitude in mV per bin per group, shape (n_groups, 720). counts_per_bin : np.ndarray @@ -2098,10 +2095,8 @@ def calculate_star_sensor_profiles_by_group( np.empty((0, 720), dtype=np.int32), ) - # Get valid data - valid_data = l1a_star["data"].values[valid_indices] - valid_counts = l1a_star["count"].values[valid_indices] - valid_epochs = l1a_star["epoch"].values[valid_indices] + # Keep valid data + l1a_star = l1a_star.isel(epoch=valid_indices) # Calculate number of groups (include partial groups) n_groups = (n_valid + group_size - 1) // group_size @@ -2113,26 +2108,25 @@ def calculate_star_sensor_profiles_by_group( # Initialize output arrays avg_amplitudes = np.zeros((n_groups, 720), dtype=np.float64) counts_per_bin = np.zeros((n_groups, 720), dtype=np.int32) - group_epochs = np.zeros(n_groups, dtype=np.int64) + group_mets = np.zeros((n_groups,), dtype=np.int64) # Process each group for group_idx in range(n_groups): start_idx = group_idx * group_size end_idx = min(start_idx + group_size, n_valid) - group_data = valid_data[start_idx:end_idx] - group_counts = valid_counts[start_idx:end_idx] + group_data = l1a_star.isel(epoch=slice(start_idx, end_idx)) # Calculate profile for this group avg_amp, count_arr = calculate_star_sensor_profile_for_group( - group_data, group_counts, edge_bins_to_exclude + group_data["data"], group_data["count"], end_bins_to_exclude ) avg_amplitudes[group_idx] = avg_amp counts_per_bin[group_idx] = count_arr - group_epochs[group_idx] = valid_epochs[start_idx] + group_mets[group_idx] = group_data["shcoarse"].values[0] - return spin_angle, group_epochs, avg_amplitudes, counts_per_bin + return spin_angle, group_mets, avg_amplitudes, counts_per_bin def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: @@ -2206,13 +2200,17 @@ def l1b_star( logger.info(f"Using spin duration from spin data: {spin_duration:.6f} s") # TODO: Read from ancillary config file when available - start_angle_offset = 62.0 # 90° - 28° - edge_bins_to_exclude = 2 + lo_angle_offset = 2 + sc_to_inst_angle_offset = ( + 360 * get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) + + lo_angle_offset + ) + end_bins_to_exclude = 2 # Calculate profiles for each 64-spin group ( spin_angle, - group_epochs, + group_mets, avg_amplitudes, counts_per_bin, ) = calculate_star_sensor_profiles_by_group( @@ -2220,24 +2218,16 @@ def l1b_star( sampling_cadence, spin_duration, group_size=group_size, - start_angle_offset=start_angle_offset, - edge_bins_to_exclude=edge_bins_to_exclude, + start_angle_offset=sc_to_inst_angle_offset, + end_bins_to_exclude=end_bins_to_exclude, ) - # Sort data so spin_angle is monotonically increasing from 0 to 360 - sort_indices = np.argsort(spin_angle) - spin_angle_sorted = spin_angle[sort_indices] - # Apply sorting to all groups' data - avg_amplitudes_sorted = avg_amplitudes[:, sort_indices] - counts_per_bin_sorted = counts_per_bin[:, sort_indices] - # Original bin indices, reordered to match the sorted spin_angle - original_bin_indices = sort_indices.astype(np.uint16) - # Get global epoch times from L1A data for start_doy and end_doy global_start_epoch = l1a_star["epoch"].values[0] global_end_epoch = l1a_star["epoch"].values[-1] # Create dataset with spin_angle as coordinate and multiple epochs + group_epochs = met_to_ttj2000ns(group_mets) l1b_star_ds = xr.Dataset( coords={ "epoch": xr.DataArray( @@ -2246,7 +2236,7 @@ def l1b_star( attrs=attr_mgr_l1b.get_variable_attributes("epoch"), ), "spin_angle": xr.DataArray( - spin_angle_sorted, + spin_angle, dims=["spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes( "spin_angle", check_schema=False @@ -2257,55 +2247,52 @@ def l1b_star( ) # Add spin_angle_bin as a variable (original bin indices) - # Broadcast to all epochs since bin mapping is the same for all groups l1b_star_ds["spin_angle_bin"] = xr.DataArray( - np.broadcast_to(original_bin_indices, (len(group_epochs), 720)), - dims=["epoch", "spin_angle"], - attrs=attr_mgr_l1b.get_variable_attributes("spin_angle_bin"), + np.arange(720, dtype=np.uint16), + dims=["spin_angle"], + attrs=attr_mgr_l1b.get_variable_attributes( + "spin_angle_bin", check_schema=False + ), + ) + + l1b_star_ds["met"] = xr.DataArray( + group_mets, + dims=["epoch"], + attrs=attr_mgr_l1b.get_variable_attributes("met"), ) l1b_star_ds["avg_amplitude"] = xr.DataArray( - avg_amplitudes_sorted, + avg_amplitudes, dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("avg_amplitude"), ) l1b_star_ds["count_per_bin"] = xr.DataArray( - counts_per_bin_sorted, + counts_per_bin, dims=["epoch", "spin_angle"], attrs=attr_mgr_l1b.get_variable_attributes("count_per_bin"), ) + # Sort the dataset by spin_angle + l1b_star_ds = l1b_star_ds.sortby("spin_angle") + # Add pointing mid time (MET) as a scalar value # Use the first epoch to determine which pointing we're in - first_met = ttj2000ns_to_met(global_start_epoch) + first_met = l1a_star["shcoarse"].values[0] pointing_mid_met = get_pointing_mid_time(first_met) - l1b_star_ds["pointing_mid_met"] = xr.DataArray( - pointing_mid_met, - attrs=attr_mgr_l1b.get_variable_attributes( - "pointing_mid_met", check_schema=False - ), - ) # Add global start and end day of year as scalar values start_doy = epoch_to_fractional_doy(global_start_epoch) end_doy = epoch_to_fractional_doy(global_end_epoch) - l1b_star_ds["start_doy"] = xr.DataArray( - start_doy, - attrs=attr_mgr_l1b.get_variable_attributes("start_doy"), - ) - - l1b_star_ds["end_doy"] = xr.DataArray( - end_doy, - attrs=attr_mgr_l1b.get_variable_attributes("end_doy"), - ) - # Add processing parameters as metadata + l1b_star_ds.attrs["start_doy"] = start_doy + l1b_star_ds.attrs["end_doy"] = end_doy + l1b_star_ds.attrs["pointing_mid_met"] = pointing_mid_met l1b_star_ds.attrs["sampling_cadence_ms"] = sampling_cadence l1b_star_ds.attrs["spin_duration_sec"] = spin_duration - l1b_star_ds.attrs["start_angle_offset_deg"] = start_angle_offset - l1b_star_ds.attrs["edge_bins_excluded"] = edge_bins_to_exclude + l1b_star_ds.attrs["lo_angle_offset_deg"] = lo_angle_offset + l1b_star_ds.attrs["end_bins_excluded"] = end_bins_to_exclude l1b_star_ds.attrs["min_count_threshold"] = 700 l1b_star_ds.attrs["group_size"] = group_size diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index 57d04c8c2c..bd33e57ed1 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1455,7 +1455,7 @@ def test_raises_error_when_field_missing(self): # Act / Assert with pytest.raises( - ValueError, + KeyError, match="ifb_data_interval field not found in L1B NHK dataset", ): get_sampling_cadence_from_nhk(l1b_nhk) @@ -1584,7 +1584,7 @@ def test_profile_for_group_basic(self): # Act avg_amplitude, count_per_bin = calculate_star_sensor_profile_for_group( - data, counts, edge_bins_to_exclude=0 + data, counts, end_bins_to_exclude=0 ) # Assert @@ -1596,7 +1596,7 @@ def test_profile_for_group_basic(self): assert np.all(avg_amplitude >= 100) assert np.all(avg_amplitude <= 200) - def test_profile_for_group_edge_bins_excluded(self): + def test_profile_for_group_end_bins_excluded(self): """Test that edge bins are properly excluded.""" # Arrange - 2 records with uniform data data = np.ones((2, 720), dtype=np.uint16) * 100 @@ -1604,24 +1604,19 @@ def test_profile_for_group_edge_bins_excluded(self): # Act avg_amplitude, count_per_bin = calculate_star_sensor_profile_for_group( - data, counts, edge_bins_to_exclude=2 + data, counts, end_bins_to_exclude=2 ) # Assert - # First 2 bins and last 2 bins should have count=0 - assert count_per_bin[0] == 0 - assert count_per_bin[1] == 0 + # Last 2 bins should have count=0 assert count_per_bin[718] == 0 assert count_per_bin[719] == 0 - # Middle bins should have count=2 - assert np.all(count_per_bin[2:718] == 2) - # Edge bins should have FILLVAL - assert avg_amplitude[0] == -1.0e31 - assert avg_amplitude[1] == -1.0e31 - assert avg_amplitude[718] == -1.0e31 - assert avg_amplitude[719] == -1.0e31 + # All other bins should have count=2 + assert np.all(count_per_bin[:718] == 2) + # End bins should have FILLVAL + assert np.all(np.isnan(avg_amplitude[718:])) # Middle bins should have average value - assert np.all(avg_amplitude[2:718] == 100.0) + assert np.all(avg_amplitude[:718] == 100.0) def test_profile_for_group_empty_data(self): """Test handling of empty data array.""" @@ -1831,7 +1826,7 @@ def test_initializes_with_spin_data( assert "spin_angle_bin" in l1b_star_ds.data_vars assert "avg_amplitude" in l1b_star_ds.data_vars assert "count_per_bin" in l1b_star_ds.data_vars - assert "pointing_mid_met" in l1b_star_ds.data_vars + assert "pointing_mid_met" in l1b_star_ds.attrs # Check that spin_angle is monotonically increasing spin_angles = l1b_star_ds.coords["spin_angle"].values assert np.all(np.diff(spin_angles) > 0), ( @@ -1847,12 +1842,11 @@ def test_initializes_with_spin_data( assert l1b_star_ds.attrs["spin_duration_sec"] == 15.0 assert l1b_star_ds.attrs["group_size"] == 64 # Check data shapes - all variables have epoch as first dimension - assert l1b_star_ds["spin_angle_bin"].shape == (3, 720) + assert l1b_star_ds["spin_angle_bin"].shape == (720,) assert l1b_star_ds["avg_amplitude"].shape == (3, 720) assert l1b_star_ds["count_per_bin"].shape == (3, 720) # Check pointing_mid_met is a scalar with expected value - assert l1b_star_ds["pointing_mid_met"].dims == () - assert float(l1b_star_ds["pointing_mid_met"].values) == 1000.0 + assert float(l1b_star_ds.attrs["pointing_mid_met"]) == 1000.0 @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") @@ -1919,11 +1913,11 @@ def test_dataset_structure_and_attributes( assert l1b_star_ds["count_per_bin"].attrs["VALIDMAX"] == 100000 # Assert - Check processing parameter attributes - assert "start_angle_offset_deg" in l1b_star_ds.attrs - assert "edge_bins_excluded" in l1b_star_ds.attrs + assert "lo_angle_offset_deg" in l1b_star_ds.attrs + assert "end_bins_excluded" in l1b_star_ds.attrs assert "min_count_threshold" in l1b_star_ds.attrs - assert l1b_star_ds.attrs["start_angle_offset_deg"] == 62.0 - assert l1b_star_ds.attrs["edge_bins_excluded"] == 2 + assert l1b_star_ds.attrs["lo_angle_offset_deg"] == 2.0 + assert l1b_star_ds.attrs["end_bins_excluded"] == 2 assert l1b_star_ds.attrs["min_count_threshold"] == 700 @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") @@ -1979,28 +1973,18 @@ def test_start_and_end_doy_variables( l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b) # Assert - Check that start_doy and end_doy exist as scalars (global values) - assert "start_doy" in l1b_star_ds.data_vars - assert "end_doy" in l1b_star_ds.data_vars - - # Assert - Check they are scalar values (no dimensions) - assert l1b_star_ds["start_doy"].dims == () - assert l1b_star_ds["end_doy"].dims == () + assert "start_doy" in l1b_star_ds.attrs + assert "end_doy" in l1b_star_ds.attrs # Assert - Check values are valid day of year (1.0 to 366.x for leap years) - start_doy = float(l1b_star_ds["start_doy"].values) - end_doy = float(l1b_star_ds["end_doy"].values) + start_doy = float(l1b_star_ds.attrs["start_doy"]) + end_doy = float(l1b_star_ds.attrs["end_doy"]) assert 1.0 <= start_doy <= 367.0 assert 1.0 <= end_doy <= 367.0 # Assert - end_doy should be >= start_doy (data spans 30 seconds) assert end_doy >= start_doy - # Assert - Check attributes - assert l1b_star_ds["start_doy"].attrs["UNITS"] == "day" - assert l1b_star_ds["end_doy"].attrs["UNITS"] == "day" - assert "Fractional day of year" in l1b_star_ds["start_doy"].attrs["CATDESC"] - assert "Fractional day of year" in l1b_star_ds["end_doy"].attrs["CATDESC"] - @patch("imap_processing.lo.l1b.lo_l1b.get_pointing_mid_time") @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") def test_multiple_groups_created( @@ -2056,8 +2040,7 @@ def test_multiple_groups_created( # Assert assert len(l1b_star_ds.coords["epoch"]) == 3 # Check pointing_mid_met is present (scalar value) - assert "pointing_mid_met" in l1b_star_ds.data_vars - assert l1b_star_ds["pointing_mid_met"].dims == () + assert "pointing_mid_met" in l1b_star_ds.attrs # First group epoch should be the first L1A epoch assert l1b_star_ds.coords["epoch"].values[0] == met_to_ttj2000ns([0.0])[0] # Second group epoch should be record 64 @@ -2076,13 +2059,13 @@ def test_star_integration(use_test_repoint_data_csv): ) ) star_path = Path( - "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2025/11/imap_lo_l1a_star_20251110-repoint00044_v001.cdf" + "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2026/01/imap_lo_l1a_star_20260121-repoint00133_v001.cdf" ) spin_path = Path( - "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2025/11/imap_lo_l1a_spin_20251110-repoint00044_v001.cdf" + "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2026/01/imap_lo_l1a_spin_20260121-repoint00133_v001.cdf" ) nhk_path = Path( - "/Users/plummert/Projects/imap/data/prod/imap/lo/l1b/2025/11/imap_lo_l1b_nhk_20251110-repoint00044_v001.cdf" + "/Users/plummert/Projects/imap/data/prod/imap/lo/l1b/2026/01/imap_lo_l1b_nhk_20260121-repoint00133_v001.cdf" ) sci_dependencies = { "imap_lo_l1a_star": load_cdf(star_path), From 35bc5ffdc7105e29f15e53af478adf2a0b20013d Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Mon, 26 Jan 2026 16:45:09 -0700 Subject: [PATCH 08/11] Updates to lo star product Reimplement fractional DOY function --- .../config/imap_lo_l1b_variable_attrs.yaml | 2 + imap_processing/lo/l1b/lo_l1b.py | 37 +++--- imap_processing/spice/time.py | 33 +++--- imap_processing/tests/lo/test_lo_l1b.py | 3 +- imap_processing/tests/spice/test_time.py | 106 ++++++++++++++++++ 5 files changed, 150 insertions(+), 31 deletions(-) diff --git a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml index cb673f9d54..8c457edd37 100644 --- a/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_l1b_variable_attrs.yaml @@ -157,6 +157,7 @@ avg_amplitude: FILLVAL: -1.0000000E+31 UNITS: mV LABLAXIS: Amplitude + DEPEND_0: epoch DEPEND_1: spin_angle count_per_bin: @@ -168,4 +169,5 @@ count_per_bin: VALIDMAX: 100000 UNITS: ' ' LABLAXIS: Sample Count + DEPEND_0: epoch DEPEND_1: spin_angle diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index ab82d6c92e..177637327a 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -2000,7 +2000,7 @@ def calculate_star_sensor_profile_for_group( Number of samples accumulated per bin, shape (720,). """ if len(data) == 0: - return np.full(720, -1.0e31, dtype=np.float64), np.zeros(720, dtype=np.int32) + return np.full(720, np.nan, dtype=np.float64), np.zeros(720, dtype=np.int32) # Determine valid bin ranges for each record use_edge_exclusion = (end_bins_to_exclude > 0) & (counts > end_bins_to_exclude) @@ -2095,36 +2095,39 @@ def calculate_star_sensor_profiles_by_group( np.empty((0, 720), dtype=np.int32), ) - # Keep valid data + # Keep valid data using xarray selection l1a_star = l1a_star.isel(epoch=valid_indices) # Calculate number of groups (include partial groups) n_groups = (n_valid + group_size - 1) // group_size + last_group_size = n_valid % group_size logger.info( f"Processing {n_valid} valid records into {n_groups} groups of {group_size}" ) + if last_group_size != 0: + logger.debug(f"Last group contains {last_group_size} records (partial group)") + + # Assign group labels to the dataset for xarray groupby operations + group_labels = np.repeat(np.arange(n_groups), group_size)[:n_valid] + l1a_star = l1a_star.assign_coords(group=("epoch", group_labels)) + + # Extract first MET for each group using xarray groupby + group_mets = l1a_star["shcoarse"].groupby("group").first().values.astype(np.int64) # Initialize output arrays avg_amplitudes = np.zeros((n_groups, 720), dtype=np.float64) counts_per_bin = np.zeros((n_groups, 720), dtype=np.int32) - group_mets = np.zeros((n_groups,), dtype=np.int64) - - # Process each group - for group_idx in range(n_groups): - start_idx = group_idx * group_size - end_idx = min(start_idx + group_size, n_valid) - - group_data = l1a_star.isel(epoch=slice(start_idx, end_idx)) + # Process each group using xarray groupby + for group_label, group_data in l1a_star.groupby("group"): # Calculate profile for this group avg_amp, count_arr = calculate_star_sensor_profile_for_group( - group_data["data"], group_data["count"], end_bins_to_exclude + group_data["data"].values, group_data["count"].values, end_bins_to_exclude ) - avg_amplitudes[group_idx] = avg_amp - counts_per_bin[group_idx] = count_arr - group_mets[group_idx] = group_data["shcoarse"].values[0] + avg_amplitudes[group_label] = avg_amp + counts_per_bin[group_label] = count_arr return spin_angle, group_mets, avg_amplitudes, counts_per_bin @@ -2200,12 +2203,13 @@ def l1b_star( logger.info(f"Using spin duration from spin data: {spin_duration:.6f} s") # TODO: Read from ancillary config file when available - lo_angle_offset = 2 + lo_angle_offset = 2.0 sc_to_inst_angle_offset = ( 360 * get_spacecraft_to_instrument_spin_phase_offset(SpiceFrame.IMAP_LO) + lo_angle_offset ) end_bins_to_exclude = 2 + min_count_threshold = 700 # Calculate profiles for each 64-spin group ( @@ -2220,6 +2224,7 @@ def l1b_star( group_size=group_size, start_angle_offset=sc_to_inst_angle_offset, end_bins_to_exclude=end_bins_to_exclude, + min_count_threshold=min_count_threshold, ) # Get global epoch times from L1A data for start_doy and end_doy @@ -2293,7 +2298,7 @@ def l1b_star( l1b_star_ds.attrs["spin_duration_sec"] = spin_duration l1b_star_ds.attrs["lo_angle_offset_deg"] = lo_angle_offset l1b_star_ds.attrs["end_bins_excluded"] = end_bins_to_exclude - l1b_star_ds.attrs["min_count_threshold"] = 700 + l1b_star_ds.attrs["min_count_threshold"] = min_count_threshold l1b_star_ds.attrs["group_size"] = group_size logger.info( diff --git a/imap_processing/spice/time.py b/imap_processing/spice/time.py index 1b91972232..e8f06c2d94 100644 --- a/imap_processing/spice/time.py +++ b/imap_processing/spice/time.py @@ -409,31 +409,36 @@ def epoch_to_doy(epoch: np.ndarray) -> npt.NDArray: ) -def epoch_to_fractional_doy(epoch_ttj2000ns: int) -> float: +def epoch_to_fractional_doy(epoch: int | Iterable[int]) -> float | np.ndarray: """ - Convert epoch in TTJ2000ns to floating point day of year. + Convert epoch in TTJ2000ns to floating point day-of-year. + + Uses SPICE's timout function to directly extract day of year and # codespell:ignore + time components, avoiding intermediate datetime parsing. Parameters ---------- - epoch_ttj2000ns : int + epoch : int Epoch in TTJ2000ns format (nanoseconds since J2000). Returns ------- doy : float Floating point day of year (1.0 = Jan 1 00:00:00). + + References + ---------- + https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/timout_c.html """ - # Convert to ephemeris time, then to UTC string - et = ttj2000ns_to_et(epoch_ttj2000ns) - utc_str = et_to_utc(et) # Returns ISO format: "YYYY-MM-DDTHH:MM:SS.sss" + # Convert to ephemeris time (ET/TDB) + et = ttj2000ns_to_et(epoch) - # Parse the datetime (remove trailing 'Z' if present) - dt = datetime.fromisoformat(utc_str.rstrip("Z")) + def single_et_to_fractional_doy(et: float) -> float: # numpydoc ignore=GL08 + # Use SPICE timout to extract DOY and time components # codespell:ignore + # Format: "DOY.####" with ::UTC modifier for UTC-based output + # The ::UTC modifier converts ET to UTC but doesn't appear in output + return float(spiceypy.timout(et, "DOY.#### ::UTC", 30)) # codespell:ignore - # Calculate day of year as floating point - # Day of year starts at 1, so Jan 1 00:00:00 = 1.0 - start_of_year = datetime(dt.year, 1, 1) - delta = dt - start_of_year - doy = 1.0 + delta.total_seconds() / 86400.0 + vectorized_et_to_frac_doy = _vectorize(single_et_to_fractional_doy, otypes=[float]) - return doy + return vectorized_et_to_frac_doy(et) diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index bd33e57ed1..a5e02cafd7 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1631,7 +1631,8 @@ def test_profile_for_group_empty_data(self): # Assert np.testing.assert_array_equal(count_per_bin, np.zeros(720)) - np.testing.assert_array_equal(avg_amplitude, np.full(720, -1.0e31)) + # Empty data returns NaN for all bins (consistent with bins having no samples) + assert np.all(np.isnan(avg_amplitude)) @patch("imap_processing.lo.l1b.lo_l1b.interpolate_repoint_data") def test_profiles_by_group_creates_correct_groups(self, mock_repoint): diff --git a/imap_processing/tests/spice/test_time.py b/imap_processing/tests/spice/test_time.py index 87f5f0e501..2dd62d0bf5 100644 --- a/imap_processing/tests/spice/test_time.py +++ b/imap_processing/tests/spice/test_time.py @@ -8,6 +8,7 @@ from imap_processing.spice.time import ( TICK_DURATION, epoch_to_doy, + epoch_to_fractional_doy, et_to_datetime64, et_to_met, et_to_ttj2000ns, @@ -297,3 +298,108 @@ def test_ttj2000ns_to_met(): ttj2000ns_array = met_to_ttj2000ns(met_array) roundtrip_met_array = ttj2000ns_to_met(ttj2000ns_array) np.testing.assert_array_almost_equal(roundtrip_met_array, met_array) + + +class TestEpochToFractionalDoy: + """Tests for epoch_to_fractional_doy function.""" + + def test_january_first_midnight(self): + """Test that January 1st 00:00:00 returns exactly 1.0.""" + # January 1st at midnight should be DOY 1.0 + utc = "2025-01-01T00:00:00.000" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + assert doy == 1.0 + + def test_january_first_noon(self): + """Test that January 1st 12:00:00 returns 1.5.""" + # January 1st at noon should be DOY 1.5 + utc = "2025-01-01T12:00:00.000" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + np.testing.assert_almost_equal(doy, 1.5, decimal=6) + + def test_known_mid_year_date(self): + """Test a known mid-year date (July 1st = DOY 182 in non-leap year).""" + # July 1st 00:00:00 in 2025 (non-leap year) should be DOY 182.0 + utc = "2025-07-01T00:00:00.000" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + # July 1st is day 182 in a non-leap year + # Jan(31) + Feb(28) + Mar(31) + Apr(30) + May(31) + Jun(30) = 181 days + # So July 1 is day 182 + assert doy == 182.0 + + def test_leap_year_date(self): + """Test leap year handling (Feb 29th exists in 2024).""" + # March 1st 00:00:00 in 2024 (leap year) should be DOY 61.0 + utc = "2024-03-01T00:00:00.000" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + # In leap year: Jan(31) + Feb(29) = 60 days, so March 1 is day 61 + assert doy == 61.0 + + def test_end_of_year(self): + """Test December 31st returns DOY 365 (non-leap) or 366 (leap).""" + # December 31st 00:00:00 in 2025 (non-leap year) should be DOY 365.0 + utc = "2025-12-31T00:00:00.000" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + assert doy == 365.0 + + def test_fractional_time_components(self): + """Test that hours, minutes, and seconds contribute correctly.""" + # January 2nd at 06:30:00 should be DOY 2.0 + 6.5/24 = 2.270833... + utc = "2025-01-02T06:30:00.000" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + # Expected: 2.0 + 6/24 + 30/1440 = 2.0 + 0.25 + 0.02083... = 2.270833... + expected_doy = 2.0 + 6.0 / 24.0 + 30.0 / 1440.0 + np.testing.assert_almost_equal(doy, expected_doy, decimal=4) + + def test_subsecond_precision(self): + """Test that subsecond precision is preserved.""" + # January 1st at 00:00:00.5 should be DOY 1.0 + 0.5/86400 + utc = "2025-01-01T00:00:00.500" + et = spiceypy.str2et(utc) + epoch = int(spiceypy.unitim(et, "ET", "TT") * 1e9) + + doy = epoch_to_fractional_doy(epoch) + + # Expected: 1.0 + 0.5/86400 = 1.00000578703... + expected_doy = 1.0 + 0.5 / 86400.0 + np.testing.assert_almost_equal(doy, expected_doy, decimal=4) + + def test_array_input(self): + """Test that np.array as input works.""" + # Test various dates throughout the year + test_dates = [ + "2025-01-01T00:00:00.000", + "2025-06-15T12:30:45.123", + "2025-12-31T23:59:59.999", + "2024-02-29T12:00:00.000", # Leap year + ] + + ets = spiceypy.str2et(test_dates) + epochs = et_to_ttj2000ns(ets) + doys = epoch_to_fractional_doy(epochs) + assert np.all(doys >= 1.0) + assert np.all(doys < 367.0) From 40c2776d723423ad578ef1284163f00be4a4fc5e Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Tue, 27 Jan 2026 17:29:40 -0700 Subject: [PATCH 09/11] PR feedback rename l1b_star to l1b_prostar --- .../cdf/config/imap_lo_global_cdf_attrs.yaml | 6 ---- imap_processing/lo/l1b/lo_l1b.py | 4 +-- imap_processing/tests/lo/test_lo_l1b.py | 33 ++----------------- 3 files changed, 4 insertions(+), 39 deletions(-) diff --git a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml index 4527e7c5eb..942a317ba9 100644 --- a/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml +++ b/imap_processing/cdf/config/imap_lo_global_cdf_attrs.yaml @@ -83,12 +83,6 @@ imap_lo_l1b_prostar: Logical_source: imap_lo_l1b_prostar Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1B Data -imap_lo_l1b_star: - <<: *instrument_base - Data_type: L1B_star>Level-1B Star Sensor Spin-Averaged Profile - Logical_source: imap_lo_l1b_star - Logical_source_description: IMAP Mission IMAP-Lo Instrument Level-1B Star Sensor Spin-Averaged Profile Data - imap_lo_l1b_nhk: <<: *instrument_base Data_type: L1B_star>Level-1B Nominal Housekeeping diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index 177637327a..b70c34de3c 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -2158,7 +2158,7 @@ def get_sampling_cadence_from_nhk(l1b_nhk: xr.Dataset) -> float: # Get mean value across all epochs (should be relatively constant) sampling_cadence = float(l1b_nhk["ifb_data_interval"].values.mean()) - logger.info(f"Sampling cadence from NHK: {sampling_cadence:.3f} ms") + logger.info(f"Star sensor sampling cadence from NHK: {sampling_cadence:.3f} ms") return sampling_cadence @@ -2189,7 +2189,7 @@ def l1b_star( L1B star sensor dataset with spin_angle, avg_amplitude, count_per_bin, and time range metadata. Each epoch corresponds to a group of records. """ - logical_source = "imap_lo_l1b_star" + logical_source = "imap_lo_l1b_prostar" l1a_star = sci_dependencies["imap_lo_l1a_star"] l1b_nhk = sci_dependencies["imap_lo_l1b_nhk"] spin_data = sci_dependencies["imap_lo_l1a_spin"] diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index a5e02cafd7..64c055db87 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -1,5 +1,4 @@ from collections import namedtuple -from pathlib import Path from unittest.mock import patch import numpy as np @@ -9,7 +8,7 @@ from imap_processing import imap_module_directory from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes -from imap_processing.cdf.utils import load_cdf, write_cdf +from imap_processing.cdf.utils import load_cdf from imap_processing.lo.l1b.lo_l1b import ( calculate_de_rates, calculate_histogram_rates, @@ -1816,7 +1815,7 @@ def test_initializes_with_spin_data( l1b_star_ds = l1b_star(sci_dependencies, attr_mgr_l1b, group_size=64) # Assert - assert l1b_star_ds.attrs["Logical_source"] == "imap_lo_l1b_star" + assert l1b_star_ds.attrs["Logical_source"] == "imap_lo_l1b_prostar" assert "epoch" in l1b_star_ds.coords # 150 records / 64 group_size = 3 groups (64 + 64 + 22) assert len(l1b_star_ds.coords["epoch"]) == 3 @@ -2050,31 +2049,3 @@ def test_multiple_groups_created( assert ( l1b_star_ds.coords["epoch"].values[2] == met_to_ttj2000ns([128 * 15.0])[0] ) - - -def test_star_integration(use_test_repoint_data_csv): - """Temporary integration test for star data.""" - use_test_repoint_data_csv( - Path( - "/Users/plummert/Projects/imap/data/prod/imap/spice/repoint/imap_2026_022_01.repoint" - ) - ) - star_path = Path( - "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2026/01/imap_lo_l1a_star_20260121-repoint00133_v001.cdf" - ) - spin_path = Path( - "/Users/plummert/Projects/imap/data/prod/imap/lo/l1a/2026/01/imap_lo_l1a_spin_20260121-repoint00133_v001.cdf" - ) - nhk_path = Path( - "/Users/plummert/Projects/imap/data/prod/imap/lo/l1b/2026/01/imap_lo_l1b_nhk_20260121-repoint00133_v001.cdf" - ) - sci_dependencies = { - "imap_lo_l1a_star": load_cdf(star_path), - "imap_lo_l1a_spin": load_cdf(spin_path), - "imap_lo_l1b_nhk": load_cdf(nhk_path), - } - anc_dependencies = [] - descriptor = "star" - result = lo_l1b(sci_dependencies, anc_dependencies, descriptor) - assert len(result) == 1 - print(write_cdf(result[0])) From 48d1803f5cff9e6da52e8e5f4a5842113d2dc440 Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Tue, 27 Jan 2026 17:30:06 -0700 Subject: [PATCH 10/11] Apply suggestion from @Copilot Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- imap_processing/spice/time.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/imap_processing/spice/time.py b/imap_processing/spice/time.py index e8f06c2d94..32af85676a 100644 --- a/imap_processing/spice/time.py +++ b/imap_processing/spice/time.py @@ -418,13 +418,16 @@ def epoch_to_fractional_doy(epoch: int | Iterable[int]) -> float | np.ndarray: Parameters ---------- - epoch : int - Epoch in TTJ2000ns format (nanoseconds since J2000). + epoch : int or Iterable[int] + Epoch in TTJ2000ns format (nanoseconds since J2000). Can be a single + integer or an iterable of integers. Returns ------- - doy : float - Floating point day of year (1.0 = Jan 1 00:00:00). + doy : float or numpy.ndarray + Floating point day of year (1.0 = Jan 1 00:00:00). Returns a scalar + when `epoch` is a single integer, or a NumPy array when `epoch` is an + iterable. References ---------- From e558d3737a10ee174fbc9fb6a1dbb937317bce5b Mon Sep 17 00:00:00 2001 From: Tim Plummer Date: Tue, 27 Jan 2026 17:33:35 -0700 Subject: [PATCH 11/11] Precommit fix --- imap_processing/lo/l1b/lo_l1b.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index b70c34de3c..447ae02da3 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -20,9 +20,8 @@ SpiceFrame, cartesian_to_latitudinal, frame_transform, - lo_instrument_pointing, get_spacecraft_to_instrument_spin_phase_offset, - instrument_pointing, + lo_instrument_pointing, ) from imap_processing.spice.repoint import ( get_pointing_mid_time,