ENH/MNT: Allow array-like time conversion

.github/workflows/test.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [windows-2019, ubuntu-latest, macos-latest]
+        os: [windows-latest, ubuntu-latest, macos-latest]
         python-version: ['3.9', '3.10', '3.11', '3.12']
     defaults:
       run:

imap_processing/__init__.py

@@ -14,8 +14,6 @@
 # This directory is used by the imap_processing package to find the packet definitions.
 from pathlib import Path
 
-import numpy as np
-
 from imap_processing._version import __version__, __version_tuple__  # noqa: F401
 
 # Eg. imap_module_directory = /usr/local/lib/python3.11/site-packages/imap_processing
@@ -34,6 +32,3 @@
     "swe": ["l0", "l1a", "l1b", "l2"],
     "ultra": ["l0", "l1a", "l1b", "l1c", "l2"],
 }
-
-# Reference start time (launch time or epoch)
-launch_time = np.datetime64("2010-01-01T00:01:06.184", "ns")

imap_processing/cdf/config/imap_hi_variable_attrs.yaml

@@ -70,7 +70,7 @@ default_epoch: &default_epoch
   TIME_BASE: J2000
   TIME_SCALE: Terrestrial Time
   REFERENCE_POSITION: Rotating Earth Geoid
-  dtype: datetime64[ns]
+  dtype: int64
 
 # ------- L1A DE Section -------
 hi_de_ccsds_met:

imap_processing/cdf/utils.py

@@ -16,39 +16,49 @@
 logger = logging.getLogger(__name__)
 
 
-def calc_start_time(
-    shcoarse_time: float,
-    launch_time: Optional[np.datetime64] = imap_processing.launch_time,
-) -> np.datetime64:
-    """
-    Calculate the datetime64 from the CCSDS secondary header information.
+# Reference start time (launch time or epoch)
+# DEFAULT_EPOCH = np.datetime64("2010-01-01T00:01:06.184", "ns")
+IMAP_EPOCH = np.datetime64("2010-01-01T00:00:00", "ns")
+J2000_EPOCH = np.datetime64("2000-01-01T11:58:55.816", "ns")
+
 
-    Since all instrument has SHCOARSE or MET seconds, we need convert it to
-    UTC. Took this from IDEX code.
+def met_to_j2000ns(
+    met: np.typing.ArrayLike,
+    reference_epoch: Optional[np.datetime64] = IMAP_EPOCH,
+) -> np.typing.ArrayLike:
+    """
+    Convert mission elapsed time (MET) to nanoseconds from J2000.
 
     Parameters
     ----------
-    shcoarse_time : float
-        Number of seconds since epoch (nominally the launch time).
-    launch_time : np.datetime64
-        The time of launch to use as the baseline.
+    met : array_like
+        Number of seconds since epoch according to the spacecraft clock.
+    reference_epoch : np.datetime64
+        The time of reference for the mission elapsed time. The standard
+        reference time for IMAP is January 1, 2010 00:00:00 UTC. Per APL's
+        IMAP Timekeeping System Design document.
 
     Returns
     -------
-    np.timedelta64
-        The time of the event.
+    array_like or scalar, int64
+        The mission elapsed time converted to nanoseconds since the J2000 epoch.
 
     Notes
     -----
-    TODO - move this into imap-data-access? How should it be used?
-    This conversion is temporary for now, and will need SPICE in the future.
-    Nick Dutton mentioned that s/c clock start epoch is
-        jan-1-2010-00:01:06.184 ET
-    We will use this for now.
+    This conversion is temporary for now, and will need SPICE in the future to
+    account for spacecraft clock drift.
     """
-    # Get the datetime of Jan 1 2010 as the start date
-    time_delta = np.timedelta64(int(shcoarse_time * 1e9), "ns")
-    return launch_time + time_delta
+    # Mission elapsed time is in seconds, convert to nanoseconds
+    # NOTE: We need to multiply the incoming met by 1e9 first because we could have
+    #       float input and we want to keep ns precision in those floats
+    # NOTE: We need int64 here when running on 32bit systems as plain int will default
+    #       to 32bit and overflow due to the nanosecond multiplication
+    time_array = (np.asarray(met, dtype=float) * 1e9).astype(np.int64)
+    # Calculate the time difference between our reference system and J2000
+    j2000_offset = (
+        (reference_epoch - J2000_EPOCH).astype("timedelta64[ns]").astype(np.int64)
+    )
+    return j2000_offset + time_array
 
 
 def load_cdf(
@@ -64,19 +74,14 @@ def load_cdf(
     remove_xarray_attrs : bool
         Whether to remove the xarray attributes that get injected by the
         cdf_to_xarray function from the output xarray.Dataset. Default is True.
-    **kwargs : {dict} optional
-        Keyword arguments for ``cdf_to_xarray``. This function overrides the
-        ``cdf_to_xarray`` default keyword value `to_datetime=False` with
-        ``to_datetime=True`.
+    **kwargs : dict, optional
+        Keyword arguments for ``cdf_to_xarray``.
 
     Returns
     -------
     dataset : xarray.Dataset
         The ``xarray`` dataset for the CDF file.
     """
-    # TODO: remove this when cdflib is updated to version >1.3.0
-    if "to_datetime" not in kwargs:
-        kwargs["to_datetime"] = True
     dataset = cdf_to_xarray(file_path, kwargs)
 
     # cdf_to_xarray converts single-value attributes to lists
@@ -121,9 +126,9 @@ def write_cdf(dataset: xr.Dataset):
     # Create the filename from the global attributes
     # Logical_source looks like "imap_swe_l2_counts-1min"
     instrument, data_level, descriptor = dataset.attrs["Logical_source"].split("_")[1:]
-    start_time = np.datetime_as_string(dataset["epoch"].values[0], unit="D").replace(
-        "-", ""
-    )
+    # Convert J2000 epoch referenced data to datetime64
+    dt64 = J2000_EPOCH + dataset["epoch"].values[0].astype("timedelta64[ns]")
+    start_time = np.datetime_as_string(dt64, unit="D").replace("-", "")
 
     # Will now accept vXXX or XXX formats, as batch starter sends versions as vXXX.
     r = re.compile(r"v\d{3}")
@@ -159,7 +164,6 @@ def write_cdf(dataset: xr.Dataset):
     xarray_to_cdf(
         dataset,
         str(file_path),
-        datetime64_to_cdftt2000=True,
         terminate_on_warning=True,
     )  # Terminate if not ISTP compliant
 

imap_processing/cli.py

@@ -618,15 +618,15 @@ def do_processing(self, dependencies):
         elif self.data_level == "l1b":
             data_dict = {}
             for dependency in dependencies:
-                dataset = load_cdf(dependency, to_datetime=True)
+                dataset = load_cdf(dependency)
                 data_dict[dataset.attrs["Logical_source"]] = dataset
             dataset = lo_l1b.lo_l1b(data_dict, self.version)
             return [dataset]
 
         elif self.data_level == "l1c":
             data_dict = {}
             for dependency in dependencies:
-                dataset = load_cdf(dependency, to_datetime=True)
+                dataset = load_cdf(dependency)
                 data_dict[dataset.attrs["Logical_source"]] = dataset
             dataset = lo_l1c.lo_l1c(data_dict, self.version)
             return [dataset]

imap_processing/codice/codice_l1a.py

@@ -22,7 +22,7 @@
 
 from imap_processing import imap_module_directory
 from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes
-from imap_processing.cdf.utils import calc_start_time
+from imap_processing.cdf.utils import IMAP_EPOCH, met_to_j2000ns
 from imap_processing.codice import constants
 from imap_processing.codice.codice_l0 import decom_packets
 from imap_processing.codice.utils import CODICEAPID, create_hskp_dataset
@@ -78,18 +78,16 @@ def __init__(self, table_id: int, plan_id: int, plan_step: int, view_id: int):
         self.plan_step = plan_step
         self.view_id = view_id
 
-    def create_science_dataset(
-        self, start_time: np.datetime64, data_version: str
-    ) -> xr.Dataset:
+    def create_science_dataset(self, met: np.int64, data_version: str) -> xr.Dataset:
         """
         Create an ``xarray`` dataset for the unpacked science data.
 
         The dataset can then be written to a CDF file.
 
         Parameters
         ----------
-        start_time : numpy.datetime64
-            The start time of the packet, used to determine epoch data variable.
+        met : numpy.int64
+            The mission elapsed time of the packet, used to determine epoch data.
         data_version : str
             Version of the data product being created.
 
@@ -106,10 +104,7 @@ def create_science_dataset(
 
         # Define coordinates
         epoch = xr.DataArray(
-            [
-                start_time,
-                start_time + np.timedelta64(1, "s"),
-            ],  # TODO: Fix after SIT-3 (see note below)
+            met_to_j2000ns(met),  # TODO: Fix after SIT-3 (see note below)
             name="epoch",
             dims=["epoch"],
             attrs=cdf_attrs.get_variable_attributes("epoch_attrs"),
@@ -380,11 +375,8 @@ def process_codice_l1a(file_path: Path | str, data_version: str) -> xr.Dataset:
                 packets = sort_by_time(grouped_data[apid], "SHCOARSE")
 
                 # Determine the start time of the packet
-                start_time = calc_start_time(
-                    packets[0].data["ACQ_START_SECONDS"].raw_value,
-                    launch_time=np.datetime64("2010-01-01T00:01:06.184", "ns"),
-                )
-
+                met = packets[0].data["ACQ_START_SECONDS"].raw_value
+                met = [met, met + 1]  # TODO: Remove after SIT-3
                 # Extract the data
                 science_values = packets[0].data["DATA"].raw_value
 
@@ -397,7 +389,7 @@ def process_codice_l1a(file_path: Path | str, data_version: str) -> xr.Dataset:
                 pipeline.get_acquisition_times()
                 pipeline.get_data_products(apid)
                 pipeline.unpack_science_data(science_values)
-                dataset = pipeline.create_science_dataset(start_time, data_version)
+                dataset = pipeline.create_science_dataset(met, data_version)
 
     # TODO: Temporary workaround in order to create hi data products in absence
     #       of simulated data. This is essentially the same process as is for
@@ -417,15 +409,15 @@ def process_codice_l1a(file_path: Path | str, data_version: str) -> xr.Dataset:
             apid = CODICEAPID.COD_HI_SECT_SPECIES_COUNTS
             table_id, plan_id, plan_step, view_id = (1, 0, 0, 6)
 
-        start_time = np.datetime64(
-            "2024-04-29T00:00:00", "ns"
-        )  # Using this to match the other data products
+        met0 = (np.datetime64("2024-04-29T00:00") - IMAP_EPOCH).astype("timedelta64[s]")
+        met0 = met0.astype(np.int64)
+        met = [met0, met0 + 1]  # Using this to match the other data products
         science_values = ""  # Currently don't have simulated data for this
 
         pipeline = CoDICEL1aPipeline(table_id, plan_id, plan_step, view_id)
         pipeline.get_data_products(apid)
         pipeline.unpack_science_data(science_values)
-        dataset = pipeline.create_science_dataset(start_time, data_version)
+        dataset = pipeline.create_science_dataset(met, data_version)
 
     # Write dataset to CDF
     logger.info(f"\nFinal data product:\n{dataset}\n")

imap_processing/codice/utils.py

@@ -13,7 +13,7 @@
 
 from imap_processing.cdf.global_attrs import ConstantCoordinates
 from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes
-from imap_processing.cdf.utils import calc_start_time
+from imap_processing.cdf.utils import met_to_j2000ns
 
 
 class CODICEAPID(IntEnum):
@@ -129,12 +129,10 @@ def create_hskp_dataset(packets, data_version: str) -> xr.Dataset:
 
     # TODO: Is there a way to get the attrs from the YAML-based method?
     epoch = xr.DataArray(
-        [
-            calc_start_time(
-                item, launch_time=np.datetime64("2010-01-01T00:01:06.184", "ns")
-            )
-            for item in metadata_arrays["SHCOARSE"]
-        ],
+        met_to_j2000ns(
+            metadata_arrays["SHCOARSE"],
+            reference_epoch=np.datetime64("2010-01-01T00:01:06.184", "ns"),
+        ),
         name="epoch",
         dims=["epoch"],
         attrs=ConstantCoordinates.EPOCH,

imap_processing/glows/l1a/glows_l1a.py

@@ -8,7 +8,7 @@
 import xarray as xr
 
 from imap_processing.cdf.global_attrs import ConstantCoordinates
-from imap_processing.cdf.utils import calc_start_time
+from imap_processing.cdf.utils import J2000_EPOCH, met_to_j2000ns
 from imap_processing.glows import __version__, glows_cdf_attrs
 from imap_processing.glows.l0.decom_glows import decom_packets
 from imap_processing.glows.l0.glows_l0_data import DirectEventL0
@@ -46,7 +46,7 @@ def glows_l1a(packet_filepath: Path, data_version: str) -> list[xr.Dataset]:
         hist_l1a = HistogramL1A(hist)
         # Split by IMAP start time
         # TODO: Should this be MET?
-        hist_day = calc_start_time(hist.SEC).astype("datetime64[D]")
+        hist_day = (J2000_EPOCH + met_to_j2000ns(hist.SEC)).astype("datetime64[D]")
         hists_by_day[hist_day].append(hist_l1a)
 
     # Generate CDF files for each day
@@ -85,7 +85,7 @@ def process_de_l0(
     de_by_day = dict()
 
     for de in de_l0:
-        de_day = calc_start_time(de.MET).astype("datetime64[D]")
+        de_day = (J2000_EPOCH + met_to_j2000ns(de.MET)).astype("datetime64[D]")
         if de_day not in de_by_day:
             de_by_day[de_day] = [DirectEventL1A(de)]
         elif de.SEQ != 0:
@@ -163,7 +163,7 @@ def generate_de_dataset(
 
     for index, de in enumerate(de_l1a_list):
         # Set the timestamp to the first timestamp of the direct event list
-        epoch_time = calc_start_time(de.l0.MET).astype("datetime64[ns]")
+        epoch_time = met_to_j2000ns(de.l0.MET).astype("datetime64[ns]")
 
         # determine if the length of the direct_events numpy array is long enough,
         # and extend the direct_events length dimension if necessary.
@@ -332,7 +332,7 @@ def generate_histogram_dataset(
 
     for index, hist in enumerate(hist_l1a_list):
         # TODO: Should this be MET?
-        epoch_time = calc_start_time(hist.imap_start_time.to_seconds())
+        epoch_time = met_to_j2000ns(hist.imap_start_time.to_seconds())
         hist_data[index] = hist.histograms
 
         support_data["flags_set_onboard"].append(hist.flags["flags_set_onboard"])

imap_processing/hi/l1a/histogram.py

@@ -7,7 +7,7 @@
 from space_packet_parser.parser import Packet
 
 from imap_processing.cdf.global_attrs import ConstantCoordinates
-from imap_processing.cdf.utils import calc_start_time
+from imap_processing.cdf.utils import met_to_j2000ns
 from imap_processing.hi import hi_cdf_attrs
 
 # TODO: Verify that these names are OK for counter variables in the CDF
@@ -57,9 +57,7 @@ def create_dataset(packets: list[Packet]) -> xr.Dataset:
 
     # unpack the packets data into the Dataset
     for i_epoch, packet in enumerate(packets):
-        dataset.epoch.data[i_epoch] = calc_start_time(
-            packet.data["CCSDS_MET"].raw_value
-        )
+        dataset.epoch.data[i_epoch] = met_to_j2000ns(packet.data["CCSDS_MET"].raw_value)
         dataset.ccsds_met[i_epoch] = packet.data["CCSDS_MET"].raw_value
         dataset.esa_step[i_epoch] = packet.data["ESA_STEP"].raw_value
 

imap_processing/hi/l1a/science_direct_event.py

@@ -4,9 +4,10 @@
 import xarray as xr
 from space_packet_parser.parser import Packet
 
-from imap_processing import imap_module_directory, launch_time
+from imap_processing import imap_module_directory
 from imap_processing.cdf.cdf_attribute_manager import CdfAttributeManager
 from imap_processing.cdf.global_attrs import ConstantCoordinates
+from imap_processing.cdf.utils import met_to_j2000ns
 
 # TODO: read LOOKED_UP_DURATION_OF_TICK from
 # instrument status summary later. This value
@@ -20,25 +21,6 @@
 MICROSECOND_TO_NS = 1e3
 
 
-def get_direct_event_time(time_in_ns: int) -> np.datetime64:
-    """
-    Create MET(Mission Elapsed Time) time using input times.
-
-    Parameters
-    ----------
-    time_in_ns : int
-        Time in nanoseconds.
-
-    Returns
-    -------
-    met_datetime : numpy.datetime64
-        Human-readable MET time.
-    """
-    met_datetime = launch_time + np.timedelta64(int(time_in_ns), "ns")
-
-    return met_datetime
-
-
 def parse_direct_event(event_data: str) -> dict:
     """
     Parse event data.
@@ -267,7 +249,7 @@ def create_dataset(de_data_list: list, packet_met_time: list) -> xr.Dataset:
             + event["de_tag"] * LOOKED_UP_DURATION_OF_TICK * MICROSECOND_TO_NS
         )
         data_dict["event_met"].append(de_met_in_ns)
-        data_dict["epoch"].append(get_direct_event_time(de_met_in_ns))
+        data_dict["epoch"].append(met_to_j2000ns(de_met_in_ns / 1e9))
         data_dict["esa_stepping_num"].append(current_esa_step)
         # start_bitmask_data is 1, 2, 3 for detector A, B, C
         # respectively. This is used to identify which detector