Mag L2 Calibration and offsets (IMAP-Science-Operations-Center#1616)

* L2 first pass

* Setting up empty mag L2

* First pass at calibration

* Updates

* updates

* Finished first draft

* Fixing filenames

* Updating docs

* Move config file to python file, add additional info

* Fix docstring

Author: maxinelasp

imap_processing/cli.py

@@ -842,6 +842,7 @@ def do_processing(
 
         if self.data_level == "l2":
             # TODO: Overwrite dependencies with versions from offsets file
+            # TODO: Ensure that parent_files attribute works with that
             input_data = load_cdf(science_files[0])
             # TODO: use ancillary from input
             calibration_dataset = load_cdf(
@@ -850,10 +851,18 @@ def do_processing(
                 / "mag"
                 / "validation"
                 / "calibration"
-                / "imap_mag_l1b-calibration_20240229_v001.cdf"
+                / "imap_mag_l2-calibration-matrices_20251017_v004.cdf"
+            )
+
+            offset_dataset = load_cdf(
+                Path(__file__).parent
+                / "tests"
+                / "mag"
+                / "validation"
+                / "calibration"
+                / "imap_mag_l2-offsets-norm_20251017_20251017_v001.cdf"
             )
             # TODO: Test data missing
-            offset_dataset = xr.Dataset()
             datasets = mag_l2(calibration_dataset, offset_dataset, input_data)
 
         return datasets

imap_processing/mag/imap_mag_sdc-configuration_v001.yaml

@@ -0,0 +1,8 @@
+"""File describing configuration changes that are not expected to change often."""
+
+# interpolation method used in L1C processing
+L1C_INTERPOLATION_METHOD = "linear_filtered"
+
+# Setting for L2 output
+# if this setting is updated, the dependency settings for L2 will also need to change.
+ALWAYS_OUTPUT_MAGO = True

imap_processing/mag/imap_mag_sdc_configuration_v001.py

@@ -334,13 +334,14 @@ def calibrate_vector(
     updated_vector : numpy.ndarray
         Calibrated vector.
     """
-    updated_vector: np.ndarray = input_vector.copy()
+    updated_vector: np.ndarray = input_vector.copy().astype(np.float64)
     if input_vector[3] % 1 != 0:
         raise ValueError("Range must be an integer.")
 
     range = int(input_vector[3])
-    x_y_z = input_vector[:3]
+    x_y_z = updated_vector[:3]
     updated_vector[:3] = np.dot(calibration_matrix.values[:, :, range], x_y_z)
+
     return updated_vector
 
 

imap_processing/mag/l1b/mag_l1b.py

@@ -1,14 +1,13 @@
 """MAG L1C processing module."""
 
 import logging
-from pathlib import Path
 from typing import Optional
 
 import numpy as np
 import xarray as xr
-import yaml
 
 from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes
+from imap_processing.mag import imap_mag_sdc_configuration_v001 as configuration
 from imap_processing.mag.constants import ModeFlags, VecSec
 from imap_processing.mag.l1c.interpolation_methods import InterpolationFunction
 
@@ -59,12 +58,7 @@ def mag_l1c(
         first_input_dataset, second_input_dataset
     )
 
-    with open(
-        Path(__file__).parent.parent / "imap_mag_sdc-configuration_v001.yaml"
-    ) as f:
-        configuration = yaml.safe_load(f)
-
-    interp_function = InterpolationFunction[configuration["L1C_interpolation_method"]]
+    interp_function = InterpolationFunction[configuration.L1C_INTERPOLATION_METHOD]
     if normal_mode_dataset and burst_mode_dataset:
         full_interpolated_timeline = process_mag_l1c(
             normal_mode_dataset, burst_mode_dataset, interp_function

imap_processing/mag/l1c/mag_l1c.py

@@ -4,26 +4,58 @@
 import xarray as xr
 
 from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes
+from imap_processing.mag import imap_mag_sdc_configuration_v001 as configuration
 from imap_processing.mag.constants import DataMode
+from imap_processing.mag.l1b.mag_l1b import calibrate_vector
 from imap_processing.mag.l2.mag_l2_data import MagL2
 
 
 def mag_l2(
-    calibration_dataset: xr.Dataset,
-    offset_dataset: xr.Dataset,
+    calibration_datasets: list[xr.Dataset],
+    offsets_dataset: xr.Dataset,
     input_data: xr.Dataset,
 ) -> list[xr.Dataset]:
     """
     Complete MAG L2 processing.
 
+    Processing uses 4 data input sources:
+    1. Calibration dataset
+        These calibration files are ancillary files and may require multiple files to
+        cover the entire timespan. They are not expected to change often. They are used
+        to provide rotation matrices to correct the frame of the vectors. The same
+        file(s) are used for both burst and norm calculations.
+    2. Offsets dataset
+        This is one, hand-created file which must correspond exactly to an L1B
+        (for burst) or L1C (for norm) data file. For each vector, this file includes
+        offsets, timedelta, and quality flags. The offsets are added to the vectors,
+        the timedelta is used to correct the epoch time, and the quality flags are
+        directly passed into the output file.
+    3. Input data
+        This is the L1B or L1C data file. It is used to provide the vectors and epoch
+        time. It should always be MAGo in the nominal case, but it is possible that we
+        will switch permanently to using MAGi (in the case of sensor failure, for
+        example.) The offsets dataset and the input
+        data are tightly related, so the input data filename is actually retrieved from
+        the offset dataset to ensure they always match.
+    4. sdc-configuration
+        This is a local configuration file for changes we never expect to make in
+        flight. This is in the IMAP local repo because changes to these settings will
+        require other code updates to validate the changes. In L2, the only setting used
+        is "always_output_mago", which indicates whether we should always output MAGo.
+        Note that if this ever is set to False, we will need to update the dependency
+        system to set MAGi files as an upstream dependency.
+
     Input data can be burst or normal mode, but MUST match the file in offset_dataset.
     TODO: retrieve the file from offset_dataset in cli.py.
+    Calibration dataset is the same for all runs.
+
+    MAGi data is not used unless we indicate it.
 
     Parameters
     ----------
-    calibration_dataset : xr.Dataset
-        Calibration ancillary file input.
-    offset_dataset : xr.Dataset
+    calibration_datasets : list[xr.Dataset]
+        Calibration ancillary file inputs.
+    offsets_dataset : xr.Dataset
         Offset ancillary file input.
     input_data : xr.Dataset
         Input data from MAG L1C or L1B.
@@ -36,39 +68,73 @@ def mag_l2(
     """
     # TODO we may need to combine multiple calibration datasets into one timeline.
 
+    always_output_mago = configuration.ALWAYS_OUTPUT_MAGO
+
+    # TODO Check that the input file matches the offsets file
+    if not np.array_equal(input_data["epoch"].data, offsets_dataset["epoch"].data):
+        raise ValueError("Input file and offsets file must have the same timestamps.")
+
+    calibration_matrix = retrieve_matrix_from_l2_calibration(
+        calibration_datasets, always_output_mago
+    )
+
+    vectors = np.apply_along_axis(
+        func1d=calibrate_vector,
+        axis=1,
+        arr=input_data["vectors"].data,
+        calibration_matrix=calibration_matrix,
+    )
+
     basic_test_data = MagL2(
-        input_data["vectors"].data[:, :3],  # level 2 vectors don't include range
+        vectors[:, :3],  # level 2 vectors don't include range
         input_data["epoch"].data,
         input_data["vectors"].data[:, 3],
         {},
         np.zeros(len(input_data["epoch"].data)),
         np.zeros(len(input_data["epoch"].data)),
         DataMode.NORM,
+        offsets=offsets_dataset["offsets"].data,
+        timedelta=offsets_dataset["timedeltas"].data,
     )
     attributes = ImapCdfAttributes()
     attributes.add_instrument_global_attrs("mag")
     # temporarily point to l1c
     attributes.add_instrument_variable_attrs("mag", "l1c")
-
     return [basic_test_data.generate_dataset(attributes)]
 
 
-def apply_calibration_matrix(
-    calibration_dataset: xr.Dataset, vectors: np.ndarray
-) -> np.ndarray:
+def retrieve_matrix_from_l2_calibration(
+    calibration_datasets: list[xr.Dataset], use_mago: bool = True
+) -> xr.DataArray:
     """
-    Apply the calibration file to the vectors to rotate them in space.
+    Get the calibration matrix for the file.
 
     Parameters
     ----------
-    calibration_dataset : xr.Dataset
-        Ancillary file input for calibration.
-    vectors : np.ndarray
-        (n, 4) array of vectors to rotate and timeshift.
+    calibration_datasets : list[xr.Dataset]
+        Ancillary file inputs for calibration.
+    use_mago : bool
+        Use the MAGo calibration matrix. Default is True.
 
     Returns
     -------
     np.ndarray
-        Rotated and timeshifted vectors.
+        Calibration matrix in the shape (3, 3, 4) to rotate vectors.
     """
-    raise NotImplementedError
+    # TODO: allow for multiple inputs
+    if isinstance(calibration_datasets, list):
+        calibration_dataset = calibration_datasets[0]
+        if len(calibration_datasets) > 1:
+            raise NotImplementedError
+    else:
+        calibration_dataset = calibration_datasets
+
+    if use_mago:
+        calibration_data = calibration_dataset["URFTOORFO"]
+    else:
+        calibration_data = calibration_dataset["URFTOORFI"]
+
+    # TODO will need to combine multiple files here
+    # TODO: Check validity of the calibration file?
+
+    return calibration_data

imap_processing/mag/l2/mag_l2.py

@@ -1,6 +1,6 @@
 """Data structures for MAG L2 and L1D processing."""
 
-from dataclasses import dataclass, field
+from dataclasses import InitVar, dataclass, field
 from enum import Enum
 
 import numpy as np
@@ -26,6 +26,8 @@ class MagL2:
 
     Since L2 and L1D should have the same structure, this can be used for either level.
 
+    Some of the methods are also static, so they can be used in i-ALiRT processing.
+
     Attributes
     ----------
     vectors: np.ndarray
@@ -57,9 +59,27 @@ class MagL2:
     data_mode: DataMode
     magnitude: np.ndarray = field(init=False)
     is_l1d: bool = False
+    offsets: InitVar[np.ndarray] = None
+    timedelta: InitVar[np.ndarray] = None
+
+    def __post_init__(self, offsets: np.ndarray, timedelta: np.ndarray) -> None:
+        """
+        Calculate the magnitude of the vectors after initialization.
+
+        Parameters
+        ----------
+        offsets : np.ndarray
+            Offsets to apply to the vectors. Should be of shape (n, 3) where n is the
+            number of vectors.
+        timedelta : np.ndarray
+            Time deltas to shift the timestamps by. Should be of length n.
+            Given in seconds.
+        """
+        if offsets is not None:
+            self.vectors = self.apply_offsets(self.vectors, offsets)
+        if timedelta is not None:
+            self.epoch = self.shift_timestamps(self.epoch, timedelta)
 
-    def __post_init__(self) -> None:
-        """Calculate the magnitude of the vectors after initialization."""
         self.magnitude = self.calculate_magnitude(self.vectors)
 
     @staticmethod
@@ -84,6 +104,69 @@ def calculate_magnitude(
         """
         return np.zeros(vectors.shape[0])  # type: ignore
 
+    @staticmethod
+    def apply_offsets(vectors: np.ndarray, offsets: np.ndarray) -> np.ndarray:
+        """
+        Apply the offsets to the vectors by adding them together.
+
+        These offsets are used to shift the vectors in the x, y, and z directions.
+        They can either be provided through a custom offsets datafile, or calculated
+        using a gradiometry algorithm.
+
+        Parameters
+        ----------
+        vectors : np.ndarray
+            Array of vectors to apply the offsets to. Should be of shape (n, 3) where n
+            is the number of vectors.
+        offsets : np.ndarray
+            Array of offsets to apply to the vectors. Should be of shape (n, 3) where n
+            is the number of vectors.
+
+        Returns
+        -------
+        np.ndarray
+            Array of vectors with offsets applied. Should be of shape (n, 3).
+        """
+        if vectors.shape[0] != offsets.shape[0]:
+            raise ValueError("Vectors and offsets must have the same length.")
+
+        offset_vectors: np.ndarray = vectors[:, :3] + offsets
+
+        # TODO: CDF files don't have NaNs. Emailed MAG to ask what this will look like.
+        # Any values where offsets is nan must also be nan
+        offset_vectors[np.isnan(offsets).any(axis=1)] = np.nan
+
+        return offset_vectors
+
+    @staticmethod
+    def shift_timestamps(epoch: np.ndarray, timedelta: np.ndarray) -> np.ndarray:
+        """
+        Shift the timestamps by the given timedelta.
+
+        If timedelta is positive, the epochs are shifted forward in time.
+
+        Parameters
+        ----------
+        epoch : np.ndarray
+            Array of timestamps to shift. Should be of length n.
+        timedelta : np.ndarray
+            Array of time deltas to shift the timestamps by. Should be the same length
+            as epoch. Given in seconds.
+
+        Returns
+        -------
+        np.ndarray
+            Shifted timestamps.
+        """
+        if epoch.shape[0] != timedelta.shape[0]:
+            raise ValueError(
+                "Input Epoch and offsets timedeltas must be the same length."
+            )
+
+        timedelta_ns = timedelta * 1e9
+        shifted_timestamps = epoch + timedelta_ns
+        return shifted_timestamps
+
     def truncate_to_24h(self, timestamp: str) -> None:
         """
         Truncate all data to a 24 hour period.

imap_processing/mag/l2/mag_l2_data.py

@@ -72,6 +72,27 @@ def mag_test_l1b_calibration_data():
     return calibration_data
 
 
+@pytest.fixture
+def mag_test_l2_data():
+    imap_dir = Path(__file__).parent
+    cal_file = (
+        imap_dir
+        / "validation"
+        / "calibration"
+        / "imap_mag_l2-calibration-matrices_20251017_v004.cdf"
+    )
+    calibration_data = load_cdf(cal_file)
+
+    offsets_data = load_cdf(
+        imap_dir
+        / "validation"
+        / "calibration"
+        / "imap_mag_l2-offsets-norm_20251017_20251017_v001.cdf"
+    )
+
+    return calibration_data, offsets_data
+
+
 def mag_generate_l1b_from_csv(df, logical_source):
     length = len(df.index)
     dataset = mag_l1a_dataset_generator(length)