Mag l1c updates (IMAP-Science-Operations-Center#1542)

* Initial setup

* updates

* gap finding and timeline creation

* First pass of interpolation stuff

* finished code

* CDF file creation

* CDF file stuff, adding documentation

* removing changes to cli

* add test function

* Addressing PR

* Fix bug

* updating interpolation methods

* Updating interpolation methods

* renaming

* Switching to enums for flags

* Updates to pass rate through to interpolation steps

* updates

* Fixing tests

* fix cic filter

* fix cic filter

* Updates to L1c

* Fixing tests

* PR comment updates

Author: maxinelasp

imap_processing/cli.py

@@ -819,18 +819,18 @@ def do_processing(
             datasets = [mag_l1b(input_data, self.version)]
 
         if self.data_level == "l1c":
-            # L1C depends on matching norm/burst files: eg burst-magi and norm-magi or
-            # burst-mago and norm-mago
-            if len(dependency_list) != 2:
+            input_data = [load_cdf(dep) for dep in dependencies]
+            # Input datasets can be in any order, and are validated within mag_l1c
+            if len(input_data) == 1:
+                datasets = [mag_l1c(input_data[0], self.version)]
+            elif len(input_data) == 2:
+                datasets = [mag_l1c(input_data[0], self.version, input_data[1])]
+            else:
                 raise ValueError(
                     f"Invalid dependencies found for MAG L1C:"
-                    f"{dependency_list}. Expected two dependencies."
+                    f"{dependencies}. Expected one or two dependencies."
                 )
 
-            input_data = [load_cdf(dep.imap_file_paths[0]) for dep in dependency_list]
-            # Input datasets can be in any order
-            datasets = [mag_l1c(input_data[0], input_data[1], self.version)]
-
         if self.data_level == "l2":
             # TODO: Overwrite dependencies with versions from offsets file
             input_data = load_cdf(dependencies[0])

imap_processing/mag/constants.py

@@ -59,6 +59,23 @@ class PrimarySensor(Enum):
     MAGI = 1
 
 
+class VecSec(Enum):
+    """Enum for all valid vector rates (Vectors per second)."""
+
+    ONE_VEC_PER_S = 1
+    TWO_VECS_PER_S = 2
+    FOUR_VECS_PER_S = 4
+    EIGHT_VECS_PER_S = 8
+    SIXTEEN_VECS_PER_S = 16
+    THIRTY_TWO_VECS_PER_S = 32
+    SIXTY_FOUR_VECS_PER_S = 64
+    ONE_TWENTY_EIGHT_VECS_PER_S = 128
+
+
+# Possible sensor rates
+POSSIBLE_RATES = [e.value for e in VecSec]
+
+
 class ModeFlags(Enum):
     """Enum for MAG mode flags: burst and normal (BURST + NORM)."""
 

imap_processing/mag/imap_mag_sdc-configuration_v001.yaml

@@ -2,5 +2,5 @@
 # but may change.
 
 # interpolation method used in L1C processing
-L1C_interpolation_method: "linear"
+L1C_interpolation_method: "linear_filtered"
 

imap_processing/mag/l1c/interpolation_methods.py

@@ -1,16 +1,25 @@
 # mypy: ignore-errors
 """Module containing interpolation methods for MAG L1C."""
 
+import logging
 from enum import Enum
+from typing import Optional
 
 import numpy as np
 from scipy.interpolate import make_interp_spline
+from scipy.signal import lfilter
+
+from imap_processing.mag.constants import POSSIBLE_RATES, VecSec
+
+logger = logging.getLogger(__name__)
 
 
 def linear(
     input_vectors: np.ndarray,
     input_timestamps: np.ndarray,
     output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec] = None,
+    output_rate: Optional[VecSec] = None,
 ) -> np.ndarray:
     """
     Linear interpolation of input vectors to output timestamps.
@@ -25,6 +34,10 @@ def linear(
         vectors.
     output_timestamps : numpy.ndarray
         Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Not required for this interpolation method.
+    output_rate : VecSec, optional
+        Not required for this interpolation method.
 
     Returns
     -------
@@ -40,6 +53,8 @@ def quadratic(
     input_vectors: np.ndarray,
     input_timestamps: np.ndarray,
     output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec] = None,
+    output_rate: Optional[VecSec] = None,
 ) -> np.ndarray:
     """
     Quadratic interpolation of input vectors to output timestamps.
@@ -54,6 +69,10 @@ def quadratic(
         vectors.
     output_timestamps : numpy.ndarray
         Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Not required for this interpolation method.
+    output_rate : VecSec, optional
+        Not required for this interpolation method.
 
     Returns
     -------
@@ -69,6 +88,8 @@ def cubic(
     input_vectors: np.ndarray,
     input_timestamps: np.ndarray,
     output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec] = None,
+    output_rate: Optional[VecSec] = None,
 ) -> np.ndarray:
     """
     Cubic interpolation of input vectors to output timestamps.
@@ -83,6 +104,10 @@ def cubic(
         vectors.
     output_timestamps : numpy.ndarray
         Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Not required for this interpolation method.
+    output_rate : VecSec, optional
+        Not required for this interpolation method.
 
     Returns
     -------
@@ -94,10 +119,105 @@ def cubic(
     return spline(output_timestamps)
 
 
+def estimate_rate(timestamps: np.ndarray) -> VecSec:
+    """
+    Given a set of timestamps, estimate the rate of the timestamps.
+
+    This rate will be one of the defined rates in the VecSec enum. The calculation
+    assumes there are no significant gaps in the timestamps.
+
+    Parameters
+    ----------
+    timestamps : numpy.ndarray
+        1D array of timestamps to estimate the rate of.
+
+    Returns
+    -------
+    VecSec
+        Estimated rate of the timestamps.
+    """
+    samples_per_second = timestamps.shape[0] / (timestamps[-1] - timestamps[0]) * 1e9
+    per_second = VecSec(
+        POSSIBLE_RATES[(np.abs(POSSIBLE_RATES - samples_per_second)).argmin()]
+    )
+
+    return per_second
+
+
+def cic_filter(
+    input_vectors: np.ndarray,
+    input_timestamps: np.ndarray,
+    output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec],
+    output_rate: Optional[VecSec],
+):
+    """
+    Apply CIC filter to data before interpolating.
+
+    The filtering uses a Cascaded integrator-comb (CIC) filter which is used in FSW to
+    filter down the raw data to telemetered data.
+
+    This assumes that the input_vectors and input_timestamps are downsampled to
+    the output_timestamps rate. Neither input_timestamps nor output_timestamps should
+    have significant gaps.
+
+    After the CIC filter is applied, the timestamps at the beginning and end of the
+    output are invalid. Therefore, we must pass in extra values and remove them after
+    applying the filter. This needs about double the input values to work.
+
+    Parameters
+    ----------
+    input_vectors : numpy.ndarray
+        Input vectors of shape (n, 3) where n is equal to the number of input
+        timestamps. Contains x, y, z components of the vector.
+    input_timestamps : numpy.ndarray
+        Input timestamps of shape (n,) which correspond to the timestamps of the input
+        vectors.
+    output_timestamps : numpy.ndarray
+        Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Expected rate of input timestamps.
+    output_rate : VecSec, optional
+        Expected rate of output timestamps.
+
+    Returns
+    -------
+    input_filtered, vectors_filtered : tuple[numpy.ndarray]
+        Filtered input timestamps and filtered input vectors.
+    """
+    # output rate should always be higher
+    input_rate = estimate_rate(input_timestamps) if input_rate is None else input_rate
+    output_rate = (
+        estimate_rate(output_timestamps) if output_rate is None else output_rate
+    )
+
+    if input_rate.value <= output_rate.value:
+        raise ValueError(
+            f"Burst mode input rate {input_rate} should never be less than "
+            f"the normal mode output rate {output_rate}. "
+            f"Both rates are required"
+        )
+
+    decimation_factor = int(input_rate.value / output_rate.value)
+    cic1 = np.ones(decimation_factor)
+    cic1 = cic1 / decimation_factor
+    cic2 = np.convolve(cic1, cic1)
+    delay = (len(cic2) - 1) // 2
+    input_filtered = input_timestamps
+    if delay != 0:
+        input_filtered = input_timestamps[:-delay]
+
+    vectors_filtered = lfilter(cic2, 1, input_vectors, axis=0)[delay:]
+
+    return input_filtered, vectors_filtered
+
+
 def linear_filtered(
     input_vectors: np.ndarray,
     input_timestamps: np.ndarray,
     output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec] = None,
+    output_rate: Optional[VecSec] = None,
 ) -> np.ndarray:
     """
     Linear filtered interpolation of input vectors to output timestamps.
@@ -112,20 +232,31 @@ def linear_filtered(
         vectors.
     output_timestamps : numpy.ndarray
         Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Expected rate of input timestamps to be passed into the CIC filter. If not
+        provided, this will be estimated.
+    output_rate : VecSec, optional
+        Expected rate of output timestamps to be passed into the CIC filter. If not
+        provided, this will be estimated.
 
     Returns
     -------
     numpy.ndarray
         Interpolated vectors of shape (m, 3) where m is equal to the number of output
         timestamps. Contains x, y, z components of the vector.
     """
-    pass
+    input_filtered, vectors_filtered = cic_filter(
+        input_vectors, input_timestamps, output_timestamps, input_rate, output_rate
+    )
+    return linear(vectors_filtered, input_filtered, output_timestamps)
 
 
 def quadratic_filtered(
     input_vectors: np.ndarray,
     input_timestamps: np.ndarray,
     output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec] = None,
+    output_rate: Optional[VecSec] = None,
 ) -> np.ndarray:
     """
     Quadratic filtered interpolation of input vectors to output timestamps.
@@ -140,20 +271,31 @@ def quadratic_filtered(
         vectors.
     output_timestamps : numpy.ndarray
         Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Expected rate of input timestamps to be passed into the CIC filter. If not
+        provided, this will be estimated.
+    output_rate : VecSec, optional
+        Expected rate of output timestamps to be passed into the CIC filter. If not
+        provided, this will be estimated.
 
     Returns
     -------
     numpy.ndarray
         Interpolated vectors of shape (m, 3) where m is equal to the number of output
         timestamps. Contains x, y, z components of the vector.
     """
-    pass
+    input_filtered, vectors_filtered = cic_filter(
+        input_vectors, input_timestamps, output_timestamps, input_rate, output_rate
+    )
+    return quadratic(vectors_filtered, input_filtered, output_timestamps)
 
 
 def cubic_filtered(
     input_vectors: np.ndarray,
     input_timestamps: np.ndarray,
     output_timestamps: np.ndarray,
+    input_rate: Optional[VecSec] = None,
+    output_rate: Optional[VecSec] = None,
 ) -> np.ndarray:
     """
     Cubic filtered interpolation of input vectors to output timestamps.
@@ -168,14 +310,23 @@ def cubic_filtered(
         vectors.
     output_timestamps : numpy.ndarray
         Output timestamps of shape (m,) to generate interpolated vectors for.
+    input_rate : VecSec, optional
+        Expected rate of input timestamps to be passed into the CIC filter. If not
+        provided, this will be estimated.
+    output_rate : VecSec, optional
+        Expected rate of output timestamps to be passed into the CIC filter. If not
+        provided, this will be estimated.
 
     Returns
     -------
     numpy.ndarray
         Interpolated vectors of shape (m, 3) where m is equal to the number of output
         timestamps. Contains x, y, z components of the vector.
     """
-    pass
+    input_filtered, vectors_filtered = cic_filter(
+        input_vectors, input_timestamps, output_timestamps, input_rate, output_rate
+    )
+    return cubic(vectors_filtered, input_filtered, output_timestamps)
 
 
 class InterpolationFunction(Enum):