diff --git a/imap_processing/ena_maps/ena_maps.py b/imap_processing/ena_maps/ena_maps.py index ed7b67523..77ef15702 100644 --- a/imap_processing/ena_maps/ena_maps.py +++ b/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): diff --git a/imap_processing/tests/ena_maps/test_ena_maps.py b/imap_processing/tests/ena_maps/test_ena_maps.py index 617aa1035..7a6049c30 100644 --- a/imap_processing/tests/ena_maps/test_ena_maps.py +++ b/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), diff --git a/imap_processing/tests/ultra/mock_data.py b/imap_processing/tests/ultra/mock_data.py index 4a4bae451..5a410e182 100644 --- a/imap_processing/tests/ultra/mock_data.py +++ b/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: [