Use pset midpoint in map projection at l2

imap_processing/ena_maps/ena_maps.py

@@ -22,7 +22,7 @@
 # so we define an enum to handle the coordinate names.
 from imap_processing.ena_maps.utils.coordinates import CoordNames
 from imap_processing.spice import geometry
-from imap_processing.spice.time import ttj2000ns_to_et
+from imap_processing.spice.time import met_to_ttj2000ns, ttj2000ns_to_et
 
 logger = logging.getLogger(__name__)
 
@@ -136,14 +136,15 @@ def match_coords_to_indices(
     if isinstance(input_object, PointingSet) and isinstance(output_object, PointingSet):
         raise ValueError("Cannot match indices between two PointingSet objects.")
 
-    # If event_et is not specified, use epoch of the PointingSet, if present.
+    # If event_et is not specified, use the midpoint of the PointingSet, if
+    # present.
     # The epoch will be in units of terrestrial time (TT) J2000 nanoseconds,
     # which must be converted to ephemeris time (ET) for SPICE.
     if event_et is None:
         if isinstance(input_object, PointingSet):
-            event_et = ttj2000ns_to_et(input_object.epoch)
+            event_et = input_object.midpoint_j2000_et
         elif isinstance(output_object, PointingSet):
-            event_et = ttj2000ns_to_et(output_object.epoch)
+            event_et = output_object.midpoint_j2000_et
         else:
             raise ValueError(
                 "Event time must be specified if both objects are SkyMaps."
@@ -301,6 +302,19 @@ def epoch(self) -> int:
         """
         return self.data["epoch"].values[0]
 
+    @property
+    def midpoint_j2000_et(self) -> float:
+        """
+        The midpoint of the pointing in ET.
+
+        Returns
+        -------
+        midpoint: int
+            The midpoint value [J2000 ET] of the pointing set.
+        """
+        epoch_delta = self.data["epoch_delta"].values[0]
+        return float(ttj2000ns_to_et(self.epoch + epoch_delta / 2))
+
     @property
     def unwrapped_dims_dict(self) -> dict[str, tuple[str, ...]]:
         """
@@ -696,6 +710,21 @@ def __init__(self, dataset: xr.Dataset):
         # Update az_el_points using the base class method
         self.update_az_el_points()
 
+    @property
+    def midpoint_j2000_et(self) -> float:
+        """
+        The midpoint of the pointing in ET.
+
+        Returns
+        -------
+        midpoint: int
+            The midpoint value [J2000 ET] of the pointing set.
+        """
+        epoch_delta = met_to_ttj2000ns(
+            self.data["pointing_end_met"].data - self.data["pointing_start_met"].data
+        )
+        return float(ttj2000ns_to_et(self.epoch + epoch_delta / 2))
+
 
 # Define the Map classes
 class AbstractSkyMap(ABC):

imap_processing/tests/ena_maps/test_ena_maps.py

@@ -20,6 +20,7 @@
 from imap_processing.ena_maps.utils import spatial_utils
 from imap_processing.ena_maps.utils.coordinates import CoordNames
 from imap_processing.spice import geometry
+from imap_processing.spice.time import met_to_ttj2000ns, ttj2000ns_to_et
 
 
 @pytest.fixture(autouse=True, scope="module")
@@ -76,7 +77,10 @@ def test_instantiate(self):
                 ultra_pset.num_points,
                 hp.nside2npix(self.nside),
             )
-
+            # check that the midpoint_j2000_et property is equal to the expected value
+            assert ultra_pset.midpoint_j2000_et == ttj2000ns_to_et(
+                ultra_pset.epoch + self.l1c_pset_products[0].epoch_delta[0] / 2
+            )
             # Check the repr exists
             assert "UltraPointingSet" in repr(ultra_pset)
 
@@ -148,12 +152,15 @@ class TestHiPointingSet:
     def test_init(self, hi_pset_cdf_path):
         """Test coverage for __init__ method."""
         pset_ds = load_cdf(hi_pset_cdf_path)
+        delta = 10
+        pset_ds["epoch_delta"] = ("epoch", np.array([delta]))  # Add dummy epoch delta
         hi_pset = ena_maps.HiPointingSet(pset_ds)
         assert isinstance(hi_pset, ena_maps.HiPointingSet)
         assert hi_pset.spice_reference_frame == geometry.SpiceFrame.IMAP_HAE
         assert hi_pset.num_points == 3600
         np.testing.assert_array_equal(hi_pset.az_el_points.shape, (3600, 2))
-
+        # check that the midpoint_j2000_et property is equal to the expected value
+        assert hi_pset.midpoint_j2000_et == ttj2000ns_to_et(hi_pset.epoch + delta / 2)
         for var_name in ["exposure_factor", "bg_rate", "bg_rate_sys_err"]:
             assert var_name in hi_pset.data
 
@@ -164,7 +171,9 @@ def test_from_cdf(self, hi_pset_cdf_path):
 
     def test_plays_nice_with_rectangular_sky_map(self, hi_pset_cdf_path):
         """Test that HiPointingSet works with RectangularSkyMap"""
-        hi_pset = ena_maps.HiPointingSet(hi_pset_cdf_path)
+        hi_ds = load_cdf(hi_pset_cdf_path)
+        hi_ds["epoch_delta"] = ("epoch", np.array([0]))  # Add dummy epoch delta
+        hi_pset = ena_maps.HiPointingSet(hi_ds)
         rect_map = ena_maps.RectangularSkyMap(
             spacing_deg=2, spice_frame=geometry.SpiceFrame.IMAP_HAE
         )
@@ -213,6 +222,16 @@ def lo_pset_ds():
         dims=["epoch"],
         name="epoch",
     )
+    dataset.coords["pointing_start_met"] = xr.DataArray(
+        [1],
+        dims=["epoch"],
+        name="epoch",
+    )
+    dataset.coords["pointing_end_met"] = xr.DataArray(
+        [10],
+        dims=["epoch"],
+        name="epoch",
+    )
     dataset.coords["spin_angle"] = xr.DataArray(
         [i for i in range(3600)],
         dims=["spin_angle"],
@@ -254,6 +273,14 @@ def test_init(self, lo_pset_ds):
         assert lo_pset.num_points == 144000
         np.testing.assert_array_equal(lo_pset.az_el_points.shape, (144000, 2))
 
+        # check that the midpoint_j2000_et property is equal to the expected value
+        assert lo_pset.midpoint_j2000_et == ttj2000ns_to_et(
+            lo_pset.epoch
+            + met_to_ttj2000ns(
+                lo_pset_ds.pointing_end_met - lo_pset_ds.pointing_start_met
+            )
+            / 2
+        )
         for var_name in ["exposure_time", "h_counts"]:
             assert var_name in lo_pset.data
 
@@ -295,6 +322,7 @@ def create_hi_pset_with_multidim_coords(
             HiPointingSet with multi-dimensional az_el_points.
         """
         pset_ds = load_cdf(hi_pset_cdf_path)
+        pset_ds["epoch_delta"] = ("epoch", np.array([0]))  # Add dummy epoch delta
         hi_pset = ena_maps.HiPointingSet(pset_ds)
         hi_pset.data["hae_longitude"] = xr.DataArray(
             np.random.uniform(0, 360, shape),

imap_processing/tests/ultra/mock_data.py

@@ -167,6 +167,7 @@ def get_binomial_counts(distance_scaling, lat_bin, central_lat_bin):
                 ],
                 sensitivity,
             ),
+            "epoch_delta": ([CoordNames.TIME.value], np.array([10], dtype=np.float64)),
         },
         coords={
             CoordNames.TIME.value: [
@@ -405,6 +406,7 @@ def mock_l1c_pset_product_healpix(
                 [CoordNames.TIME.value, CoordNames.HEALPIX_INDEX.value],
                 np.full((1, npix), ImapPSETUltraFlags.NONE.value, dtype=np.uint16),
             ),
+            "epoch_delta": ([CoordNames.TIME.value], np.array([10], dtype=np.float64)),
         },
         coords={
             CoordNames.TIME.value: [