diff --git a/imap_processing/cli.py b/imap_processing/cli.py index 9f67cf60f..595decf8e 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -67,6 +67,7 @@ from imap_processing.idex.idex_l1b import idex_l1b from imap_processing.idex.idex_l2a import idex_l2a from imap_processing.idex.idex_l2b import idex_l2b +from imap_processing.lo import lo_pivot_kernel from imap_processing.lo.constants import LoConstants from imap_processing.lo.l1a import lo_l1a from imap_processing.lo.l1b import lo_l1b @@ -474,6 +475,33 @@ def _resolve_version(self, descriptor: str) -> Version: logger.warning(msg) return self.version_map.get(descriptor, self._fallback_version) + def _resolve_kernel_minor_version(self, descriptor: str) -> int: + """ + Return the minor version to use in the filename of a generated kernel. + + Parameters + ---------- + descriptor : str + The descriptor of the kernel job, e.g. "pointing-attitude". + + Returns + ------- + int + The minor version for the kernel filename. + """ + resolved_version = self._resolve_version(descriptor) + if resolved_version is None: + raise ValueError( + f"No version provided for {descriptor} processing. " + f"Provide a version for the '{descriptor}' descriptor in " + "the dependency JSON's version block, or a fallback --version." + ) + return ( + resolved_version.minor + if isinstance(resolved_version, Version) + else int(resolved_version.lstrip("v")) + ) + def upload_products(self, products: list[Path]) -> None: """ Upload data products to the IMAP SDC. @@ -1431,9 +1459,48 @@ def pre_processing(self) -> ProcessingInputCollection: return filtered_dependencies + def _generate_pivot_kernel( + self, dependencies: ProcessingInputCollection + ) -> list[Path]: + """ + Generate the Lo pivot platform CK for the pointing given by repointing. + + Parameters + ---------- + dependencies : ProcessingInputCollection + Object containing dependencies to process. + + Returns + ------- + list[Path] + The generated Lo pivot kernel. + """ + if self.repointing is None: + raise ValueError( + "repointing must be provided for pivot-ckernel processing." + ) + nhk_files = dependencies.get_file_paths( + source="lo", data_type="l1b", descriptor="nhk" + ) + if len(nhk_files) != 1: + raise ValueError( + f"Unexpected dependencies found for IMAP-Lo pivot-ckernel: " + f"{nhk_files}. Expected exactly one L1B NHK file." + ) + # The repoint table provides the pointing start and end times. + if not dependencies.get_file_paths(data_type=RepointInput.data_type): + raise ValueError( + "A repoint table dependency is required for IMAP-Lo pivot-ckernel " + "processing." + ) + minor_version = self._resolve_kernel_minor_version(self.descriptor) + return lo_pivot_kernel.generate_lo_pivot_kernel( + nhk_files[0], self.repointing, minor_version + ) + def do_processing( self, dependencies: ProcessingInputCollection - ) -> list[xr.Dataset]: + ) -> list[xr.Dataset | Path]: """ Perform IMAP-Lo specific processing. @@ -1444,12 +1511,15 @@ def do_processing( Returns ------- - dataset : xr.Dataset - Xr.Dataset of output files. + datasets : list[xarray.Dataset | Path] + The list of processed products. """ print(f"Processing IMAP-Lo {self.data_level}") - datasets: list[xr.Dataset] = [] - if self.data_level == "l1a": + datasets: list[xr.Dataset | Path] = [] + if self.data_level == "l1b" and self.descriptor == "pivot-ckernel": + datasets.extend(self._generate_pivot_kernel(dependencies)) + + elif self.data_level == "l1a": # L1A packet / products are 1 to 1. Should only have # one dependency file science_files = dependencies.get_file_paths(source="lo", data_type="l0") @@ -1860,18 +1930,7 @@ def do_processing( data_type=SPICESource.SPICE.value ) ah_paths = [path for path in spice_inputs if ".ah" in path.suffixes] - resolved_version = self._resolve_version(self.descriptor) - if resolved_version is None: - raise ValueError( - "No version provided for pointing-attitude processing. " - "Provide a version for the 'pointing-attitude' descriptor in " - "the dependency JSON's version block, or a fallback --version." - ) - minor_version = ( - resolved_version.minor - if isinstance(resolved_version, Version) - else int(resolved_version.lstrip("v")) - ) + minor_version = self._resolve_kernel_minor_version(self.descriptor) pointing_kernel_paths = pointing_frame.generate_pointing_attitude_kernel( ah_paths, self.start_date, minor_version ) diff --git a/imap_processing/lo/l1b/lo_l1b.py b/imap_processing/lo/l1b/lo_l1b.py index d73ff02c9..c48e681a7 100644 --- a/imap_processing/lo/l1b/lo_l1b.py +++ b/imap_processing/lo/l1b/lo_l1b.py @@ -1862,6 +1862,44 @@ def get_pivot_angle_from_nhk(ds_nhk: xr.Dataset) -> float: return ds_nhk["pcc_cumulative_cnt_pri"].isel(epoch=nitems // 2).item() +def get_median_pivot_angle(ds_nhk: xr.Dataset) -> float: + """ + Get the median pivot angle from the NHK dataset. + + The median of ``pcc_coarse_pot_pri`` is taken over the samples between + ``PIVOT_HK_HOUR_RANGE`` hours after the first NHK sample, which avoids the + pivot platform motion at the start of a pointing. + + Parameters + ---------- + ds_nhk : xr.Dataset + The NHK dataset containing pivot angle information. + + Returns + ------- + pivot_angle : float + The median pivot angle [degrees], or NaN if the dataset has no records + or there are no valid samples within the time range. + """ + if ds_nhk.sizes.get("epoch", 0) == 0: + return np.nan + + hk_epoch_ets = ttj2000ns_to_et(ds_nhk["epoch"]) + start_et_hk = ( + hk_epoch_ets[0] + timedelta(hours=c.PIVOT_HK_HOUR_RANGE[0]).total_seconds() + ) + end_et_hk = ( + hk_epoch_ets[0] + timedelta(hours=c.PIVOT_HK_HOUR_RANGE[1]).total_seconds() + ) + + coarse_pot_pri = ds_nhk["pcc_coarse_pot_pri"].values + return float( + np.nanmedian( + coarse_pot_pri[(hk_epoch_ets >= start_et_hk) & (hk_epoch_ets <= end_et_hk)] + ) + ) + + def _get_esa_level_indices(epochs: np.ndarray, anc_dependencies: list) -> np.ndarray: """ Get the ESA level indices (reswept indices) for the given epochs. @@ -2485,18 +2523,7 @@ def l1b_bgrates_and_goodtimes( # noqa: PLR0912 pivot: float = 90.0 cdf_hk = sci_dependencies.get("imap_lo_l1b_nhk") if cdf_hk is not None and "pcc_coarse_pot_pri" in cdf_hk: - hk_epoch_ets = ttj2000ns_to_et(cdf_hk["epoch"]) - start_et_hk = ( - hk_epoch_ets[0] + timedelta(hours=c.PIVOT_HK_HOUR_RANGE[0]).total_seconds() - ) - end_et_hk = ( - hk_epoch_ets[0] + timedelta(hours=c.PIVOT_HK_HOUR_RANGE[1]).total_seconds() - ) - - coarse_pot_pri = cdf_hk["pcc_coarse_pot_pri"].values - pivot = np.nanmedian( # type: ignore - coarse_pot_pri[(hk_epoch_ets >= start_et_hk) & (hk_epoch_ets <= end_et_hk)] - ) + pivot = get_median_pivot_angle(cdf_hk) if np.isnan(pivot): pivot = 90.0 diff --git a/imap_processing/lo/lo_pivot_kernel.py b/imap_processing/lo/lo_pivot_kernel.py new file mode 100644 index 000000000..0ab3873b8 --- /dev/null +++ b/imap_processing/lo/lo_pivot_kernel.py @@ -0,0 +1,227 @@ +"""Generate the IMAP-Lo pivot platform attitude kernel (CK).""" + +import logging +import os +import tempfile +from datetime import datetime, timezone +from pathlib import Path + +import imap_data_access +import numpy as np +import spiceypy + +from imap_processing.cdf.utils import load_cdf +from imap_processing.lo.l1b.lo_l1b import get_median_pivot_angle +from imap_processing.spice.geometry import SpiceFrame +from imap_processing.spice.pointing_frame import ( + POINTING_SEGMENT_DTYPE, + segment_ck_filename, + write_constant_attitude_ck, +) +from imap_processing.spice.repoint import get_pointing_times_from_id +from imap_processing.spice.time import ( + et_to_utc, + met_to_sclkticks, + met_to_utc, + sct_to_et, +) + +logger = logging.getLogger(__name__) + +LO_PIVOT_KERNEL_PREFIX = "imap_lopivot" +LO_PIVOT_KERNEL_EXTENSION = "bc" + + +def generate_lo_pivot_kernel( + l1b_nhk_path: Path, repointing: str, minor_version: int +) -> list[Path]: + """ + Generate the IMAP-Lo pivot platform CK for a single pointing. + + The kernel contains one constant-attitude segment covering the pointing, + from the end of its repoint maneuver to the start of the next repoint + maneuver. There is no coverage during repoint maneuvers, which is when + the pivot platform is moved. + + Parameters + ---------- + l1b_nhk_path : Path + Lo L1B NHK file for the pointing. + repointing : str + The pointing to generate the kernel for, in the format 'repoint#####'. + minor_version : int + Minor version, from the batch command, to use for the output + kernel filename. + + Returns + ------- + kernel_paths : list[Path] + Location of the new Lo pivot kernel. + + Raises + ------ + ValueError + If `l1b_nhk_path` is not an NHK file for `repointing`. + FileExistsError + If the output kernel already exists. + + Notes + ----- + Kernels required to be furnished: + + - Latest NAIF leapseconds kernel (naif0012.tls) + - The latest IMAP sclk (imap_sclk_NNNN.tsc) + - The latest IMAP frame kernel (imap_###.tf), which defines IMAP_LO_BASE + + The repoint table must also be set (`imap_processing.spice.repoint`), as it + gives the pointing start and end times. + """ + repoint_id = imap_data_access.ScienceFilePath(l1b_nhk_path.name).repointing + if repoint_id is None or f"repoint{repoint_id:05d}" != repointing: + raise ValueError(f"{l1b_nhk_path.name} is not an NHK file for {repointing}.") + + segment, pivot_angle = calculate_pivot_segment(l1b_nhk_path, repoint_id) + + kernel_filename = segment_ck_filename( + f"{LO_PIVOT_KERNEL_PREFIX}-{repointing}", + segment, + minor_version, + LO_PIVOT_KERNEL_EXTENSION, + ) + kernel_path = imap_data_access.SPICEFilePath(kernel_filename).construct_path() + # open_spice_ck_file would append a duplicate segment to an existing file. + if kernel_path.exists(): + raise FileExistsError(f"Lo pivot kernel already exists: {kernel_path}") + kernel_path.parent.mkdir(parents=True, exist_ok=True) + + # Write the kernel in a temporary directory and only publish it once it is + # complete, so a failed write never leaves a partial kernel at the output + # path. os.link does not replace an existing file, preserving the + # no-overwrite behavior even if the kernel appeared during the write. + with tempfile.TemporaryDirectory(dir=kernel_path.parent) as tmp_dir: + tmp_kernel_path = Path(tmp_dir) / kernel_path.name + write_lo_pivot_ck(tmp_kernel_path, segment, pivot_angle, l1b_nhk_path.name) + os.link(tmp_kernel_path, kernel_path) + return [kernel_path] + + +def calculate_pivot_segment( + l1b_nhk_path: Path, repoint_id: int +) -> tuple[np.ndarray, float]: + """ + Calculate the data for the single segment of a Lo pivot kernel. + + Parameters + ---------- + l1b_nhk_path : Path + Lo L1B NHK file for the pointing. + repoint_id : int + Repoint ID of the pointing. + + Returns + ------- + segment : numpy.ndarray + Structured array of POINTING_SEGMENT_DTYPE with one element. The + quaternion rotates vectors from the IMAP_LO_BASE frame into the + IMAP_LO frame. + pivot_angle : float + Pivot angle [degrees] of the pointing. + + Raises + ------ + ValueError + If the NHK file has no valid pivot angle samples. + """ + pointing_start_met, pointing_end_met = get_pointing_times_from_id(repoint_id) + # Use the same pivot angle as the goodtimes product. + pivot_angle = get_median_pivot_angle(load_cdf(l1b_nhk_path)) + if np.isnan(pivot_angle): + raise ValueError(f"No valid pivot angle samples in {l1b_nhk_path.name}.") + logger.info( + f"repoint{repoint_id:05d} ({met_to_utc(pointing_start_met)}, " + f"{met_to_utc(pointing_end_met)}): pivot angle {pivot_angle:.4f} deg" + ) + + segment = np.zeros(1, dtype=POINTING_SEGMENT_DTYPE) + segment[0]["pointing_id"] = repoint_id + segment[0]["start_sclk_ticks"] = met_to_sclkticks(pointing_start_met) + segment[0]["end_sclk_ticks"] = met_to_sclkticks(pointing_end_met) + segment[0]["quaternion"] = pivot_angle_to_quaternion(pivot_angle) + return segment, pivot_angle + + +def pivot_angle_to_quaternion(pivot_angle: float) -> np.ndarray: + """ + Get the SPICE quaternion rotating IMAP_LO_BASE vectors into IMAP_LO. + + The IMAP_LO frame is the IMAP_LO_BASE frame rotated about its +X axis by + the pivot angle. The rotation is consistent with + `imap_processing.spice.geometry.get_lo_pivot_boresight`: the IMAP_LO + boresight, -Y, expressed in IMAP_LO_BASE is the pivot boresight. + + Parameters + ---------- + pivot_angle : float + The pivot angle [degrees]. + + Returns + ------- + quaternion : numpy.ndarray + SPICE-style quaternion, shape (4,). + """ + # spiceypy.rotate returns the matrix that rotates the coordinate frame by + # the angle about the axis, i.e. it maps IMAP_LO_BASE vectors into IMAP_LO. + rotation_matrix = spiceypy.rotate(np.deg2rad(pivot_angle), 1) + return np.asarray(spiceypy.m2q(rotation_matrix)) + + +def write_lo_pivot_ck( + kernel_path: Path, + segment_data: np.ndarray, + pivot_angle: float, + parent_file: str, +) -> None: + """ + Write the Lo pivot CK, recording the pivot angle in the comments. + + Parameters + ---------- + kernel_path : pathlib.Path + Location to write the CK kernel. + segment_data : numpy.ndarray + Structured array of POINTING_SEGMENT_DTYPE with one element. + pivot_angle : float + Pivot angle [degrees] of the pointing. + parent_file : str + Filename of the NHK file the pivot angle was derived from. + """ + segment = segment_data[0] + start_utc = et_to_utc(sct_to_et(segment["start_sclk_ticks"])) + end_utc = et_to_utc(sct_to_et(segment["end_sclk_ticks"])) + comments = [ + "CK FOR IMAP_LO FRAME (IMAP-LO PIVOT PLATFORM)", + "==================================================================", + "", + f"Original file name: {kernel_path.name}", + f"Creation date: {datetime.now(timezone.utc).strftime('%Y-%m-%d')}", + f"Parent files: {[parent_file]}", + "", + "The IMAP_LO frame is the IMAP_LO_BASE frame rotated about +X by the", + "pivot angle, constant over the pointing.", + "", + f"Repoint ID: repoint{segment['pointing_id']:05d}", + f"Pointing start (UTC): {start_utc}", + f"Pointing end (UTC): {end_utc}", + f"Pivot angle (deg): {pivot_angle:.4f}", + "", + ] + + logger.debug(f"Writing Lo pivot kernel: {kernel_path}") + write_constant_attitude_ck( + kernel_path, + segment_data, + SpiceFrame.IMAP_LO, + SpiceFrame.IMAP_LO_BASE, + comments, + ) + logger.debug(f"Finished writing Lo pivot kernel: {kernel_path}") diff --git a/imap_processing/spice/pointing_frame.py b/imap_processing/spice/pointing_frame.py index 2702bb120..5d9ba09d7 100644 --- a/imap_processing/spice/pointing_frame.py +++ b/imap_processing/spice/pointing_frame.py @@ -66,20 +66,9 @@ def generate_pointing_attitude_kernel( if len(pointing_segments) == 0: raise ValueError("No Pointings covered by input dependencies.") - # get the start and end yyyy_doy strings - start_datetime = spiceypy.et2datetime( - sct_to_et(pointing_segments[0]["start_sclk_ticks"]) - ) - end_datetime = spiceypy.et2datetime( - sct_to_et(pointing_segments[-1]["end_sclk_ticks"]) - ) sorted_ck_paths = list(sorted(imap_attitude_cks, key=lambda x: x.name)) - version_str = str(Version(None, minor_version)).lstrip("v") - pointing_kernel_path = ( - sorted_ck_paths[-1].parent / f"imap_dps_" - f"{start_datetime.strftime('%Y_%j')}_" - f"{end_datetime.strftime('%Y_%j')}_" - f"{version_str}.ah.bc" + pointing_kernel_path = sorted_ck_paths[-1].parent / segment_ck_filename( + "imap_dps", pointing_segments, minor_version, "ah.bc" ) write_pointing_frame_ck( pointing_kernel_path, pointing_segments, [p.name for p in imap_attitude_cks] @@ -87,6 +76,45 @@ def generate_pointing_attitude_kernel( return [pointing_kernel_path] +def segment_ck_filename( + prefix: str, segment_data: np.ndarray, minor_version: int, extension: str +) -> str: + """ + Build a CK filename from the time coverage of its segments. + + The filename has the form ``___.`` + where the dates are the start of the first segment and the end of the + last segment. + + Parameters + ---------- + prefix : str + Leading portion of the filename, e.g. "imap_dps". + segment_data : np.ndarray + Numpy structured array of POINTING_SEGMENT_DTYPE, sorted by time. + minor_version : int + Minor version to use in the filename. + extension : str + Filename extension without the leading ".", e.g. "ah.bc". + + Returns + ------- + filename : str + The CK filename. + """ + start_datetime = spiceypy.et2datetime( + sct_to_et(segment_data[0]["start_sclk_ticks"]) + ) + end_datetime = spiceypy.et2datetime(sct_to_et(segment_data[-1]["end_sclk_ticks"])) + version_str = str(Version(None, minor_version)).lstrip("v") + return ( + f"{prefix}_" + f"{start_datetime.strftime('%Y_%j')}_" + f"{end_datetime.strftime('%Y_%j')}_" + f"{version_str}.{extension}" + ) + + @contextmanager def open_spice_ck_file(pointing_frame_path: Path) -> Generator[int, None, None]: """ @@ -109,7 +137,10 @@ def open_spice_ck_file(pointing_frame_path: Path) -> Generator[int, None, None]: try: yield handle finally: - spiceypy.ckcls(handle) + # dafcls rather than ckcls: ckcls also raises SPICE(NOSEGMENTSFOUND) + # when no segment was written, which would leave the file open and + # replace the original error. + spiceypy.dafcls(handle) def write_pointing_frame_ck( @@ -142,8 +173,54 @@ def write_pointing_frame_ck( ] logger.debug(f"Writing pointing attitude kernel: {pointing_kernel_path}") + write_constant_attitude_ck( + pointing_kernel_path, + segment_data, + SpiceFrame.IMAP_DPS, + SpiceFrame.ECLIPJ2000, + comments, + ) + logger.debug(f"Finished writing pointing attitude kernel: {pointing_kernel_path}") + + +def write_constant_attitude_ck( + kernel_path: Path, + segment_data: np.ndarray, + frame: SpiceFrame, + reference_frame: SpiceFrame, + comments: list[str], +) -> None: + """ + Write a CK with a single constant attitude record per segment. + + Parameters + ---------- + kernel_path : pathlib.Path + Location to write the CK kernel. + segment_data : np.ndarray + Numpy structured array with the following dtypes: + ("start_sclk_ticks", np.float64), + ("end_sclk_ticks", np.float64), + ("quaternion", np.float64, (4,)), + ("pointing_id", np.uint32), + frame : SpiceFrame + Frame whose orientation the CK defines. Also used as the segment + identifier. + reference_frame : SpiceFrame + Frame the orientation is relative to. The quaternion in each segment + rotates vectors from this frame into `frame`. + comments : list[str] + Lines to write to the comment area of the CK. + + Raises + ------ + ValueError + If `segment_data` is empty. + """ + if len(segment_data) == 0: + raise ValueError(f"No segments to write to {kernel_path.name}.") - with open_spice_ck_file(pointing_kernel_path) as handle: + with open_spice_ck_file(kernel_path) as handle: # Write the comments to the file spiceypy.dafac(handle, comments) @@ -157,12 +234,12 @@ def write_pointing_frame_ck( segment["start_sclk_ticks"], # End time of the segment. segment["end_sclk_ticks"], - # Pointing frame ID. - SpiceFrame.IMAP_DPS.value, + # Frame ID. + frame.value, # Reference frame. - SpiceFrame.ECLIPJ2000.name, # Reference frame + reference_frame.name, # Identifier. - SpiceFrame.IMAP_DPS.name, + frame.name, # Number of pointing intervals. 1, # Start times of individual pointing records within segment. @@ -171,18 +248,16 @@ def write_pointing_frame_ck( # End times of individual pointing records within segment. # Since there is only a single record this is equal to sclk_endtim. np.array([segment["end_sclk_ticks"]]), # Single stop time - # Average quaternion. + # Constant quaternion for the segment. segment["quaternion"], - # Angular velocity vectors. The IMAP_DPS frame is quasi-inertial - # for each pointing so each segment has zeros here. + # Angular velocity vectors. The attitude is constant over each + # segment so each segment has zeros here. np.array([0.0, 0.0, 0.0]), # The number of seconds per encoded spacecraft clock # tick for each interval. np.array([TICK_DURATION]), ) - logger.debug(f"Finished writing pointing attitude kernel: {pointing_kernel_path}") - def calculate_pointing_attitude_segments( ck_paths: list[Path], diff --git a/imap_processing/tests/lo/test_lo_l1b.py b/imap_processing/tests/lo/test_lo_l1b.py index c4405e323..331b75dbd 100644 --- a/imap_processing/tests/lo/test_lo_l1b.py +++ b/imap_processing/tests/lo/test_lo_l1b.py @@ -24,6 +24,7 @@ create_datasets, filter_valid_star_records, get_avg_spin_durations_per_cycle, + get_median_pivot_angle, get_pivot_angle_from_nhk, get_sampling_cadence_from_nhk, get_spin_start_times, @@ -2213,6 +2214,43 @@ def test_get_pivot_angle_from_nhk(): assert pivot_angle == expected_pivot_angle +def _pivot_nhk(minutes: np.ndarray, pivot: np.ndarray) -> xr.Dataset: + """Make an NHK dataset sampled at the given minutes after an arbitrary t0.""" + t0_ttj2000ns = 8.2e17 + return xr.Dataset( + {"pcc_coarse_pot_pri": ("epoch", np.asarray(pivot, dtype=np.float64))}, + coords={"epoch": (t0_ttj2000ns + minutes * 60e9).astype(np.int64)}, + ) + + +def test_get_median_pivot_angle(furnish_kernels): + """Median over 0.5 h to 22.5 h after the first NHK sample.""" + minutes = np.arange(0, 24 * 60, 1.0) + pivot = np.full(minutes.shape, 75.0) + pivot[minutes < 30] = 90.0 # Pivot still moving at the start + pivot[minutes > 22.5 * 60] = 105.0 # Next repoint + pivot[100] = np.nan + with furnish_kernels(["naif0012.tls"]): + assert get_median_pivot_angle(_pivot_nhk(minutes, pivot)) == 75.0 + + +def test_get_median_pivot_angle_no_samples(furnish_kernels): + """NaN when no valid samples are in the time range.""" + minutes = np.arange(0, 20, 1.0) # Ends before the 0.5 h start + with furnish_kernels(["naif0012.tls"]): + pivot = get_median_pivot_angle(_pivot_nhk(minutes, np.full(20, 75.0))) + assert np.isnan(pivot) + + +def test_get_median_pivot_angle_empty(): + """NaN, not an IndexError, for an NHK dataset with no records.""" + empty = xr.Dataset( + {"pcc_coarse_pot_pri": ("epoch", np.array([], dtype=np.float64))}, + coords={"epoch": np.array([], dtype=np.int64)}, + ) + assert np.isnan(get_median_pivot_angle(empty)) + + def test_l1b_bgrates_and_goodtimes_basic(anc_dependencies, attr_mgr_l1b): """Test basic functionality of l1b_bgrates_and_goodtimes.""" # Arrange - Create a simple L1B histogram rates dataset diff --git a/imap_processing/tests/lo/test_lo_pivot_kernel.py b/imap_processing/tests/lo/test_lo_pivot_kernel.py new file mode 100644 index 000000000..9c5541f78 --- /dev/null +++ b/imap_processing/tests/lo/test_lo_pivot_kernel.py @@ -0,0 +1,266 @@ +"""Tests for the IMAP-Lo pivot platform kernel generation.""" + +from pathlib import Path + +import numpy as np +import pytest +import spiceypy +import xarray as xr + +from imap_processing.lo import lo_pivot_kernel +from imap_processing.lo.lo_pivot_kernel import ( + calculate_pivot_segment, + generate_lo_pivot_kernel, + pivot_angle_to_quaternion, +) +from imap_processing.spice.geometry import ( + SpiceFrame, + get_lo_pivot_boresight, + instrument_pointing, + lo_instrument_pointing, +) +from imap_processing.spice.repoint import get_pointing_times_from_id +from imap_processing.spice.time import ( + et_to_met, + met_to_sclkticks, + met_to_ttj2000ns, + sct_to_et, + str_to_et, +) + + +@pytest.fixture +def furnish_lo_pivot_kernels(furnish_kernels): + """Furnish the kernels needed to write and read the Lo pivot kernel.""" + with furnish_kernels(["naif0012.tls", "imap_sclk_0000.tsc", "imap_130.tf"]): + yield + + +def make_nhk(met: np.ndarray, pivot: np.ndarray) -> xr.Dataset: + """Make a minimal Lo L1B NHK dataset.""" + return xr.Dataset( + {"pcc_coarse_pot_pri": ("epoch", np.asarray(pivot, dtype=np.float64))}, + coords={"epoch": met_to_ttj2000ns(np.asarray(met, dtype=np.float64))}, + ) + + +@pytest.fixture +def pointings(furnish_lo_pivot_kernels, use_fake_repoint_data_for_time): + """ + Fake repoint table with two short pointings on 2025-11-10 (100 and 101) and + a pointing running from 2025-11-10 into 2025-11-11 (102). + + Returns the start and end MET of each pointing, keyed by repoint id. + """ + t0 = float(et_to_met(str_to_et("2025-11-10T01:00:00"))) + repoint_starts = t0 + np.array([0, 6, 12, 36]) * 3600.0 + use_fake_repoint_data_for_time(repoint_starts, repoint_id_start=100) + return { + repoint_id: get_pointing_times_from_id(repoint_id) + for repoint_id in [100, 101, 102] + } + + +@pytest.fixture +def nhk_files(pointings, monkeypatch, tmp_path): + """ + Lo L1B NHK files for each pointing, served by a mocked load_cdf. + + Each NHK file samples the pivot every 10 s from 10 minutes before the + pointing (the pivot move) to the end of the pointing. + """ + pivots = {100: 90.0, 101: 75.0, 102: 105.0} + datasets = {} + paths = {} + for repoint_id, (start_met, end_met) in pointings.items(): + met = np.arange(start_met - 600, end_met, 10.0) + pivot = np.where(met < start_met, 60.0, pivots[repoint_id]) + name = f"imap_lo_l1b_nhk_20251110-repoint{repoint_id:05d}_v001.cdf" + datasets[name] = make_nhk(met, pivot) + paths[repoint_id] = tmp_path / name + monkeypatch.setattr( + lo_pivot_kernel, "load_cdf", lambda path: datasets[Path(path).name] + ) + return {"paths": paths, "pivots": pivots} + + +@pytest.mark.parametrize("pivot_angle", [0.0, 60.0, 90.0, 105.0, 160.0]) +def test_pivot_angle_to_quaternion(pivot_angle): + """The IMAP_LO boresight rotated into IMAP_LO_BASE is the pivot boresight.""" + base_to_lo = spiceypy.q2m(pivot_angle_to_quaternion(pivot_angle)) + boresight_in_base = np.asarray(base_to_lo).T @ np.array([0, -1, 0]) + np.testing.assert_allclose( + boresight_in_base, get_lo_pivot_boresight(pivot_angle), atol=1e-12 + ) + + +def test_calculate_pivot_segment(nhk_files, pointings): + """One segment with the pointing coverage and the pointing's pivot angle.""" + segment, pivot_angle = calculate_pivot_segment(nhk_files["paths"][101], 101) + assert segment.shape == (1,) + assert segment[0]["pointing_id"] == 101 + assert pivot_angle == 75.0 + np.testing.assert_allclose( + [segment[0]["start_sclk_ticks"], segment[0]["end_sclk_ticks"]], + met_to_sclkticks(np.array(pointings[101])), + rtol=0, + atol=1, + ) + + +def read_comments(kernel_path: Path) -> str: + """Read the comment area of a CK.""" + handle = spiceypy.dafopr(str(kernel_path)) + try: + _, comments, _ = spiceypy.dafec(handle, 30, 200) + finally: + spiceypy.dafcls(handle) + return "\n".join(comments) + + +def test_generate_lo_pivot_kernel(nhk_files, pointings, tmp_path): + """The kernel covers the pointing with the IMAP_LO attitude.""" + kernel_paths = generate_lo_pivot_kernel(nhk_files["paths"][101], "repoint00101", 3) + assert kernel_paths == [ + tmp_path / "imap/spice/ck/imap_lopivot-repoint00101_2025_314_2025_314_003.bc" + ] + kernel_path = kernel_paths[0] + + # Coverage is exactly the pointing, excluding the repoint maneuvers. + cover = spiceypy.ckcov( + str(kernel_path), SpiceFrame.IMAP_LO.value, False, "INTERVAL", 0, "SCLK" + ) + assert spiceypy.wncard(cover) == 1 + start_met, end_met = pointings[101] + np.testing.assert_allclose( + spiceypy.wnfetd(cover, 0), + met_to_sclkticks(np.array([start_met, end_met])), + rtol=0, + atol=1, + ) + + spiceypy.furnsh(str(kernel_path)) + try: + for et in sct_to_et(met_to_sclkticks(np.linspace(start_met, end_met, 3))): + boresight = spiceypy.mxv( + spiceypy.pxform("IMAP_LO", "IMAP_LO_BASE", et), [0, -1, 0] + ) + np.testing.assert_allclose( + boresight, get_lo_pivot_boresight(75.0), atol=1e-12 + ) + # With the kernel, the IMAP_LO frame gives the same pointing as the + # pivot-angle based calculation. + et = sct_to_et(met_to_sclkticks(np.mean(pointings[101]))) + np.testing.assert_allclose( + instrument_pointing( + et, SpiceFrame.IMAP_LO, SpiceFrame.IMAP_LO_BASE, cartesian=True + ), + lo_instrument_pointing(et, 75.0, SpiceFrame.IMAP_LO_BASE, cartesian=True), + atol=1e-12, + ) + finally: + spiceypy.unload(str(kernel_path)) + + comments = read_comments(kernel_path) + assert kernel_path.name in comments + assert "imap_lo_l1b_nhk_20251110-repoint00101_v001.cdf" in comments + assert "Repoint ID: repoint00101" in comments + assert "Pivot angle (deg): 75.0000" in comments + + +def test_generate_lo_pivot_kernel_same_day_pointings(nhk_files): + """Short pointings on the same day produce distinct kernels.""" + kernel_100 = generate_lo_pivot_kernel(nhk_files["paths"][100], "repoint00100", 1) + kernel_101 = generate_lo_pivot_kernel(nhk_files["paths"][101], "repoint00101", 1) + assert kernel_100[0].name == "imap_lopivot-repoint00100_2025_314_2025_314_001.bc" + assert kernel_101[0].name == "imap_lopivot-repoint00101_2025_314_2025_314_001.bc" + assert "Pivot angle (deg): 90.0000" in read_comments(kernel_100[0]) + + +def test_generate_lo_pivot_kernel_multi_day_pointing(nhk_files): + """The end date in the filename is the end of the pointing.""" + kernel_path = generate_lo_pivot_kernel(nhk_files["paths"][102], "repoint00102", 1) + assert kernel_path[0].name == "imap_lopivot-repoint00102_2025_314_2025_315_001.bc" + + +def test_generate_lo_pivot_kernel_exists(nhk_files): + """An existing kernel is not overwritten or appended to.""" + nhk_path = nhk_files["paths"][100] + kernel_path = generate_lo_pivot_kernel(nhk_path, "repoint00100", 1)[0] + original = kernel_path.read_bytes() + with pytest.raises(FileExistsError, match=kernel_path.name): + generate_lo_pivot_kernel(nhk_path, "repoint00100", 1) + assert kernel_path.read_bytes() == original + + +def test_generate_lo_pivot_kernel_wrong_repointing(nhk_files): + """The NHK file must be for the requested pointing.""" + with pytest.raises(ValueError, match="is not an NHK file for repoint00101"): + generate_lo_pivot_kernel(nhk_files["paths"][100], "repoint00101", 1) + + +def test_calculate_pivot_segment_no_samples(nhk_files, monkeypatch): + """No valid pivot samples is an error, not a 90 deg default.""" + nhk = make_nhk(np.arange(0, 3600 * 2, 10.0), np.full(720, np.nan)) + monkeypatch.setattr(lo_pivot_kernel, "load_cdf", lambda path: nhk) + with pytest.raises(ValueError, match="No valid pivot angle samples"): + calculate_pivot_segment(nhk_files["paths"][100], 100) + + +def test_calculate_pivot_segment_empty_nhk(nhk_files, monkeypatch): + """An NHK file with no records is the same error as no valid samples.""" + nhk = make_nhk(np.array([]), np.array([])) + monkeypatch.setattr(lo_pivot_kernel, "load_cdf", lambda path: nhk) + with pytest.raises(ValueError, match="No valid pivot angle samples"): + calculate_pivot_segment(nhk_files["paths"][100], 100) + + +def test_generate_lo_pivot_kernel_write_failure(nhk_files, monkeypatch, tmp_path): + """A failed write leaves no partial kernel, so a retry succeeds.""" + + def failing_ckw02(*args, **kwargs): + raise spiceypy.utils.exceptions.SpiceyError("simulated write failure") + + ckopn = spiceypy.ckopn + handles = [] + + def recording_ckopn(*args): + handles.append(ckopn(*args)) + return handles[-1] + + with monkeypatch.context() as m: + m.setattr(spiceypy, "ckw02", failing_ckw02) + m.setattr(spiceypy, "ckopn", recording_ckopn) + with pytest.raises( + spiceypy.utils.exceptions.SpiceyError, match="simulated write failure" + ): + generate_lo_pivot_kernel(nhk_files["paths"][100], "repoint00100", 1) + + # The CK was closed, which Windows needs to delete the temporary file. + with pytest.raises(spiceypy.utils.exceptions.SpiceyError): + spiceypy.dafhsf(handles[0]) + + ck_dir = tmp_path / "imap/spice/ck" + assert list(ck_dir.iterdir()) == [] + + kernel_path = generate_lo_pivot_kernel(nhk_files["paths"][100], "repoint00100", 1) + assert [p.name for p in ck_dir.iterdir()] == [kernel_path[0].name] + + +def test_generate_lo_pivot_kernel_appears_during_write(nhk_files, monkeypatch): + """A kernel created by another process mid-write is not overwritten.""" + write = lo_pivot_kernel.write_lo_pivot_ck + final_path = {} + + def write_then_race(kernel_path, *args): + write(kernel_path, *args) + final_path["path"] = kernel_path.parent.parent / kernel_path.name + final_path["path"].write_bytes(b"other process") + + monkeypatch.setattr(lo_pivot_kernel, "write_lo_pivot_ck", write_then_race) + with pytest.raises(FileExistsError): + generate_lo_pivot_kernel(nhk_files["paths"][100], "repoint00100", 1) + assert final_path["path"].read_bytes() == b"other process" + assert [p.name for p in final_path["path"].parent.iterdir()] == [ + final_path["path"].name + ] diff --git a/imap_processing/tests/spice/test_pointing_frame.py b/imap_processing/tests/spice/test_pointing_frame.py index 76ecd5e3c..8611e2324 100644 --- a/imap_processing/tests/spice/test_pointing_frame.py +++ b/imap_processing/tests/spice/test_pointing_frame.py @@ -20,6 +20,8 @@ _mean_spin_axis, calculate_pointing_attitude_segments, generate_pointing_attitude_kernel, + open_spice_ck_file, + write_constant_attitude_ck, write_pointing_frame_ck, ) from imap_processing.spice.time import TICK_DURATION, met_to_sclkticks, sct_to_et @@ -185,6 +187,34 @@ def test_write_pointing_frame_ck( assert parent_file in lines[5] +def test_open_spice_ck_file_error_closes_file(tmp_path): + """An error before any segment is written closes the CK and is re-raised. + + ckcls would raise SPICE(NOSEGMENTSFOUND) here, leaving the file open and + replacing the original error. + """ + ck_path = tmp_path / "empty.bc" + with pytest.raises(RuntimeError, match="original error"): + with open_spice_ck_file(ck_path) as handle: + raise RuntimeError("original error") + with pytest.raises(spiceypy.utils.exceptions.SpiceyError): + spiceypy.dafhsf(handle) + + +def test_write_constant_attitude_ck_no_segments(tmp_path): + """No segments is an error, and no file is created.""" + ck_path = tmp_path / "empty.bc" + with pytest.raises(ValueError, match="No segments to write"): + write_constant_attitude_ck( + ck_path, + np.zeros(0, dtype=POINTING_SEGMENT_DTYPE), + SpiceFrame.IMAP_DPS, + SpiceFrame.ECLIPJ2000, + ["comment"], + ) + assert not ck_path.exists() + + @pytest.mark.external_test_data def test_mean_spin_axis(furnish_flight_ah_kernels): """Tests _mean_spin_axis function.""" diff --git a/imap_processing/tests/test_cli.py b/imap_processing/tests/test_cli.py index 024889ec9..5af0abb81 100644 --- a/imap_processing/tests/test_cli.py +++ b/imap_processing/tests/test_cli.py @@ -39,6 +39,9 @@ main, ) from imap_processing.spice import config as spice_config +from imap_processing.spice.geometry import get_lo_pivot_boresight +from imap_processing.spice.time import met_to_sclkticks, met_to_ttj2000ns, sct_to_et +from imap_processing.tests.conftest import generate_repoint_data @pytest.fixture(autouse=True) @@ -841,6 +844,189 @@ def test_spacecraft_pointing_kernel_no_version( assert mock_spacecraft_pointing.call_count == 0 +LO_PIVOT_DEPENDENCY_FILES = ( + '[{"type": "science","files": [' + '"imap_lo_l1b_nhk_20251110-repoint00100_v001.cdf"]}, ' + '{"type": "spice","files": ["naif0012.tls", "imap_sclk_0005.tsc", ' + '"imap_130.tf"]}, ' + '{"type": "repoint","files": ["imap_2025_315_01.repoint"]}]' +) + + +@mock.patch( + "imap_processing.cli.lo_pivot_kernel.generate_lo_pivot_kernel", autospec=True +) +def test_lo_pivot_kernel(mock_lo_pivot, mock_instrument_dependencies): + """Test coverage for the cli.Lo pivot-ckernel job""" + dependency_str = json.dumps( + { + "dependency": json.loads(LO_PIVOT_DEPENDENCY_FILES), + "version": {"pivot-ckernel": {"major_version": None, "minor_version": 4}}, + } + ) + input_collection = ProcessingInputCollection() + input_collection.deserialize(LO_PIVOT_DEPENDENCY_FILES) + kernel_path = Path("imap_lopivot-repoint00100_2025_314_2025_314_004.bc") + mock_lo_pivot.return_value = [kernel_path] + + instrument = Lo( + "l1b", "pivot-ckernel", dependency_str, None, "repoint00100", "v001", False + ) + products = instrument.do_processing(input_collection) + + assert products == [kernel_path] + assert mock_lo_pivot.call_count == 1 + nhk_path, repointing, minor_version = mock_lo_pivot.call_args[0] + assert nhk_path.name == "imap_lo_l1b_nhk_20251110-repoint00100_v001.cdf" + assert repointing == "repoint00100" + assert minor_version == 4 + + +@mock.patch( + "imap_processing.cli.lo_pivot_kernel.generate_lo_pivot_kernel", autospec=True +) +def test_lo_pivot_kernel_no_repointing(mock_lo_pivot, mock_instrument_dependencies): + """The cli.Lo pivot-ckernel job requires a repointing""" + input_collection = ProcessingInputCollection() + input_collection.deserialize(LO_PIVOT_DEPENDENCY_FILES) + instrument = Lo( + "l1b", + "pivot-ckernel", + LO_PIVOT_DEPENDENCY_FILES, + "20251110", + None, + "v001", + False, + ) + with pytest.raises(ValueError, match="repointing must be provided"): + instrument.do_processing(input_collection) + assert mock_lo_pivot.call_count == 0 + + +@mock.patch( + "imap_processing.cli.lo_pivot_kernel.generate_lo_pivot_kernel", autospec=True +) +def test_lo_pivot_kernel_multiple_nhk(mock_lo_pivot, mock_instrument_dependencies): + """The cli.Lo pivot-ckernel job requires exactly one NHK file""" + input_collection = ProcessingInputCollection() + input_collection.deserialize( + '[{"type": "science","files": [' + '"imap_lo_l1b_nhk_20251110-repoint00100_v001.cdf", ' + '"imap_lo_l1b_nhk_20251110-repoint00101_v001.cdf"]}]' + ) + instrument = Lo( + "l1b", + "pivot-ckernel", + LO_PIVOT_DEPENDENCY_FILES, + None, + "repoint00100", + "v001", + False, + ) + with pytest.raises(ValueError, match="Expected exactly one L1B NHK file"): + instrument.do_processing(input_collection) + assert mock_lo_pivot.call_count == 0 + + +@mock.patch( + "imap_processing.cli.lo_pivot_kernel.generate_lo_pivot_kernel", autospec=True +) +def test_lo_pivot_kernel_no_repoint_table(mock_lo_pivot, mock_instrument_dependencies): + """The cli.Lo pivot-ckernel job requires a repoint table""" + input_collection = ProcessingInputCollection() + input_collection.deserialize( + '[{"type": "science","files": [' + '"imap_lo_l1b_nhk_20251110-repoint00100_v001.cdf"]}]' + ) + instrument = Lo( + "l1b", + "pivot-ckernel", + LO_PIVOT_DEPENDENCY_FILES, + None, + "repoint00100", + "v001", + False, + ) + with pytest.raises(ValueError, match="A repoint table dependency is required"): + instrument.do_processing(input_collection) + assert mock_lo_pivot.call_count == 0 + + +def test_lo_pivot_kernel_process(monkeypatch, tmp_path, spice_test_data_path): + """Run the full cli.Lo pivot-ckernel job, including pre-processing. + + The dependencies already exist locally, so pre-processing does not download + anything but does furnish the kernels and set the repoint table. + """ + # Restore the module-level tables that pre-processing sets. + monkeypatch.setattr(spice_config, "_repoint_table_path", None) + monkeypatch.setattr(spice_config, "_spin_table_paths", []) + + spice_dir = tmp_path / "imap" / "spice" + for subdir, kernel in [ + ("lsk", "naif0012.tls"), + ("sclk", "imap_sclk_0000.tsc"), + ("fk", "imap_130.tf"), + ]: + (spice_dir / subdir).mkdir(parents=True) + shutil.copy(spice_test_data_path / kernel, spice_dir / subdir / kernel) + + # Repoints 100 to 102, a few hours apart, so repoint 101 is a full pointing. + repoint_starts = 500_433_000.0 + np.array([0, 6, 12]) * 3600.0 + (spice_dir / "repoint").mkdir() + generate_repoint_data(repoint_starts, repoint_id_start=100).to_csv( + spice_dir / "repoint" / "imap_2025_315_01.repoint", index=False + ) + + nhk_name = "imap_lo_l1b_nhk_20251110-repoint00101_v001.cdf" + nhk_path = imap_data_access.ScienceFilePath(nhk_name).construct_path() + nhk_path.parent.mkdir(parents=True) + nhk_path.touch() + + def fake_nhk(path): + # Built on read, once pre-processing has furnished the kernels. + met = np.arange(repoint_starts[1], repoint_starts[2], 60.0) + return xr.Dataset( + {"pcc_coarse_pot_pri": ("epoch", np.full(met.shape, 75.0))}, + coords={"epoch": met_to_ttj2000ns(met)}, + ) + + monkeypatch.setattr("imap_processing.lo.lo_pivot_kernel.load_cdf", fake_nhk) + + dependency_str = json.dumps( + { + "dependency": [ + {"type": "science", "files": [nhk_name]}, + { + "type": "spice", + "files": ["naif0012.tls", "imap_sclk_0000.tsc", "imap_130.tf"], + }, + {"type": "repoint", "files": ["imap_2025_315_01.repoint"]}, + ], + "version": {"pivot-ckernel": {"major_version": None, "minor_version": 2}}, + } + ) + Lo( + "l1b", "pivot-ckernel", dependency_str, None, "repoint00101", "v001", False + ).process() + + kernels = list((spice_dir / "ck").glob("imap_lopivot-repoint00101_*_002.bc")) + assert len(kernels) == 1 + with spiceypy.KernelPool( + [ + str(spice_dir / "lsk" / "naif0012.tls"), + str(spice_dir / "sclk" / "imap_sclk_0000.tsc"), + str(spice_dir / "fk" / "imap_130.tf"), + str(kernels[0]), + ] + ): + et = sct_to_et(met_to_sclkticks(repoint_starts[1] + 3 * 3600)) + boresight = spiceypy.mxv( + spiceypy.pxform("IMAP_LO", "IMAP_LO_BASE", et), [0, -1, 0] + ) + np.testing.assert_allclose(boresight, get_lo_pivot_boresight(75.0), atol=1e-12) + + @mock.patch("imap_processing.cli.ultra_l1a.ultra_l1a") def test_ultra_l1a(mock_ultra_l1a, mock_instrument_dependencies): """Test coverage for cli.Ultra class with l1a data level""" diff --git a/poetry.lock b/poetry.lock index 8a383fd6e..493953c9c 100644 --- a/poetry.lock +++ b/poetry.lock @@ -782,14 +782,14 @@ files = [ [[package]] name = "imap-data-access" -version = "0.42.0" +version = "1.3.0" description = "IMAP SDC Data Access" optional = false python-versions = "*" groups = ["main"] files = [ - {file = "imap_data_access-0.42.0-py3-none-any.whl", hash = "sha256:5eb6081dd9e6bcac2be60d5bb288a66d8d52922124d52e6d1ac7020beff2bc4b"}, - {file = "imap_data_access-0.42.0.tar.gz", hash = "sha256:64dbd378dbe57e35eee157b513ee4929b279056983c19bd1481eeb90dac07fbf"}, + {file = "imap_data_access-1.3.0-py3-none-any.whl", hash = "sha256:c8733496b8216f99312295240da3a0daa18344e114781bb1a58caab6813198ab"}, + {file = "imap_data_access-1.3.0.tar.gz", hash = "sha256:8ef034878c3f3145b6c0930a366106edd12ee20a9d5dab209d6160e2400a67c3"}, ] [package.dependencies] @@ -3162,4 +3162,4 @@ tools = ["openpyxl", "pandas"] [metadata] lock-version = "2.1" python-versions = ">=3.10,<3.15" -content-hash = "c2ec2d113479399719d05401d45196ed69e78b5fd54c647edd9ee29bdbcb382d" +content-hash = "f34f41059bd79accf06652aa077c863ccf9a92543282abad9fc5252c9ee5b7fb" diff --git a/pyproject.toml b/pyproject.toml index 9d3145f30..dcf0da60c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,7 @@ classifiers = [ dependencies = [ "astropy-healpix>=1.0", "cdflib>=1.3.11", - "imap-data-access>=0.42.0", + "imap-data-access>=1.3.0", "space_packet_parser>=6.0.0", "spiceypy>=6.0.0", "xarray>=2024.10.0,<2026",