Ultra geo culling (IMAP-Science-Operations-Center#1999)

Author: laspsandoval

imap_processing/tests/ultra/unit/test_ultra_l1c_culling.py

@@ -0,0 +1,97 @@
+"""Tests Culling for ULTRA L1c."""
+
+import astropy_healpix.healpy as hp
+import numpy as np
+import pytest
+import spiceypy
+
+from imap_processing.spice.geometry import SpiceBody
+from imap_processing.ultra.l1c.ultra_l1c_culling import compute_culling_mask
+
+
+@pytest.mark.external_kernel
+@pytest.mark.usefixtures("_unset_metakernel_path")
+def test_compute_culling_mask(furnish_kernels, spice_test_data_path):
+    """Tests compute_culling_mask function."""
+
+    planet_radii_km = {
+        "EARTH": 6378.137,
+    }
+
+    kernels = [
+        "imap_science_100.tf",
+        "sim_1yr_imap_pointing_frame.bc",
+        "imap_spk_demo.bsp",
+    ]
+
+    keepout_radius_km = 30 * planet_radii_km["EARTH"]
+
+    # Corresponds to 2025-11-28T00:00:00
+    et_start = 817561854.185627
+    et_end = 817644684.1856259
+    step_seconds = 1800  # 30 minutes
+    et_steps = np.arange(et_start, et_end, step_seconds)
+
+    spiceypy.kclear()
+
+    with furnish_kernels(kernels):
+        mask, _ = compute_culling_mask(et_steps, keepout_radius_km)
+
+    assert mask.shape[0] == len(et_steps)
+    assert mask.shape[1] == hp.nside2npix(128)
+
+    # Check that some pixels are masked out
+    assert not np.all(mask)
+    assert np.any(mask)
+
+
+@pytest.mark.external_kernel
+@pytest.mark.usefixtures("_unset_metakernel_path")
+def test_compare_sincpt_with_culling_mask_deterministic(furnish_kernels):
+    """Compare culling mask output for the closest-to-Earth pixel with sincpt."""
+
+    with furnish_kernels(
+        [
+            "imap_science_100.tf",
+            "imap_sclk_0000.tsc",
+            "sim_1yr_imap_pointing_frame.bc",
+            "imap_spk_demo.bsp",
+            "earth_1962_240827_2124_combined.bpc",
+            "pck00011.tpc",
+            "naif0012.tls",
+            "de440s.bsp",
+        ]
+    ):
+        et = np.array([817561854.185627])
+        keepout_radius_km = 6378.1  # Earth radius
+        nside = 128
+
+        # Compute culling mask and IMAP-to-Earth unit vector
+        mask, unit_vectors = compute_culling_mask(
+            et, keepout_radius_km, observer=SpiceBody.EARTH, nside=nside
+        )
+
+        # Computes the 3D unit vectors pointing to the centers of all HEALPix pixels
+        pixel_vecs_dps = np.column_stack(
+            hp.pix2vec(nside, np.arange(hp.nside2npix(nside)), nest=False)
+        )
+        # Find the HEALPix pixel direction closest to the direction from IMAP to Earth
+        closest_idx = np.argmax(np.dot(pixel_vecs_dps, unit_vectors[0]))
+
+        # Transform closest pixel vector to J2000
+        rot_dps_to_j2000 = spiceypy.pxform("IMAP_DPS", "J2000", et[0])
+        pixel_vec_j2000 = np.dot(rot_dps_to_j2000, pixel_vecs_dps[closest_idx])
+
+        with spiceypy.no_found_check():
+            result = spiceypy.sincpt(
+                method="ELLIPSOID",
+                target="EARTH",
+                et=et[0],
+                fixref="IAU_EARTH",
+                abcorr="NONE",
+                obsrvr="IMAP",
+                dref="J2000",
+                dvec=pixel_vec_j2000,
+            )
+
+        assert result[-1]  # Check if the ray intersects Earth

imap_processing/ultra/l1c/ultra_l1c_culling.py

@@ -0,0 +1,85 @@
+"""Culling for ULTRA L1c."""
+
+import astropy_healpix.healpy as hp
+import numpy as np
+from numpy.typing import NDArray
+
+from imap_processing.spice.geometry import (
+    SpiceBody,
+    SpiceFrame,
+    imap_state,
+)
+
+
+def compute_culling_mask(
+    et: NDArray,
+    keepout_radius_km: float,
+    observer: SpiceBody = SpiceBody.EARTH,
+    nside: int = 128,
+    nested: bool = False,
+) -> tuple[NDArray, NDArray]:
+    """
+    Compute a mask for HEALPix pixels within a keep-out radius of the target body.
+
+    Parameters
+    ----------
+    et : NDArray
+        Ephemeris times in TDB seconds past J2000.
+    keepout_radius_km : float
+        Radius (in km) within which HEALPix pixels will be excluded.
+    observer : SpiceBody, optional
+        Body from which IMAP is observed.
+    nside : int, optional
+        HEALPix NSIDE resolution. Default is 128.
+    nested : bool, optional
+        Whether to use NESTED indexing.
+
+    Returns
+    -------
+    mask : tuple[NDArray, NDArray]
+        Boolean array of shape (len(et), npix).
+    unit_target_vecs : NDArray
+        Unit vectors from IMAP to the target body
+        (e.g., Earth), shape (len(et), 3).
+    """
+    # Compute number of HEALPix pixels
+    npix = hp.nside2npix(nside)
+
+    # Compute IMAP to Earth position in the pointing frame.
+    state = imap_state(et, ref_frame=SpiceFrame.IMAP_DPS, observer=observer)
+    # Flip to get vector from IMAP to Earth
+    # position.shape = (len(et), 3)
+    position = -state[:, :3]
+
+    # Distance from IMAP to target (e.g. Earth) (km):
+    # distance.shape = (len(et),)
+    distance = np.linalg.norm(position, axis=1)  # shape (len(et),)
+
+    # Calculate the keepout angle (radians).
+    # keepout_angle.shape = (len(et),)
+    keepout_angle = np.arcsin(keepout_radius_km / distance)  # radians
+
+    # Calculate the direction from IMAP to Earth. (shape: [N, 3])
+    # unit_target_vecs.shape = (len(et), 3)
+    unit_target_vecs = position / distance[:, np.newaxis]
+
+    # Get pixel unit vectors pointing from the center of the
+    # HEALPix sphere to the center of each pixel on the sky.
+    pixel_vecs = np.column_stack(
+        hp.pix2vec(nside, np.arange(npix), nest=nested)
+    )  # shape: (npix, 3)
+
+    # Returns cos(theta) where theta is the separation angle between:
+    # (1) vector from IMAP to Earth
+    # (2) vector from IMAP to HEALPix pixel center
+    # If theta is within the keepout angle, then the pixel is culled.
+    cos_sep = np.dot(unit_target_vecs, pixel_vecs.T)  # shape (N, npix)
+    cos_sep = np.clip(cos_sep, -1.0, 1.0)
+    # Get theta here.
+    sep_angle = np.arccos(cos_sep)
+
+    # Exclude pixels within the keepout angle.
+    # mask.shape = (len(et), npix)
+    mask = sep_angle > keepout_angle[:, np.newaxis]
+
+    return mask, unit_target_vecs