CoDICE l2 angular intensity validation (#2331)

* add todos to remove pytest workarounds

Author: lacoak21

imap_processing/cdf/config/imap_codice_global_cdf_attrs.yaml

@@ -241,6 +241,24 @@ imap_codice_l2_lo-sw-species:
   Logical_source: imap_codice_l2_lo-sw-species
   Logical_source_description: IMAP Mission CoDICE Lo Level-2 Sunward Species Intensity Data.
 
+imap_codice_l2_lo-nsw-species:
+  <<: *instrument_base
+  Data_type: L2_lo-nsw-species>Level-2 Lo Non-Sunward Species Intensities Data
+  Logical_source: imap_codice_l2_lo-nsw-species
+  Logical_source_description: IMAP Mission CoDICE Lo Level-2 Non-Sunward Species Intensity Data.
+
+imap_codice_l2_lo-sw-angular:
+  <<: *instrument_base
+  Data_type: L2_lo-sw-species>Level-2 Lo Sunward Angular Intensities Data
+  Logical_source: imap_codice_l2_lo-sw-angular
+  Logical_source_description: IMAP Mission CoDICE Lo Level-2 Sunward Angular Intensity Data.
+
+imap_codice_l2_lo-nsw-angular:
+  <<: *instrument_base
+  Data_type: L2_lo-nsw-species>Level-2 Lo Non-Sunward Angular Intensities Data
+  Logical_source: imap_codice_l2_lo-nsw-angular
+  Logical_source_description: IMAP Mission CoDICE Lo Level-2 Non-Sunward Angular Intensity Data.
+
 imap_codice_l2_hi-omni:
   <<: *instrument_base
   Data_type: L2_hi-omni>Level-2 Hi Omnidirectional Species Intensities Data

imap_processing/codice/codice_l2.py

@@ -31,8 +31,9 @@
     LO_POSITION_TO_ELEVATION_ANGLE,
     LO_SW_ANGULAR_VARIABLE_NAMES,
     LO_SW_PICKUP_ION_SPECIES_VARIABLE_NAMES,
-    LO_SW_SPECIES_VARIABLE_NAMES,
+    LO_SW_SOLAR_WIND_SPECIES_VARIABLE_NAMES,
     NSW_POSITIONS,
+    PIXEL_ORIENTATIONS,
     PUI_POSITIONS,
     SOLAR_WIND_POSITIONS,
     SW_POSITIONS,
@@ -354,11 +355,14 @@ def process_lo_angular_intensity(
         positions,
         average_across_positions=False,
     )
+
     # transform positions to elevation angles
     if positions == SW_POSITIONS:
         pos_to_el = LO_POSITION_TO_ELEVATION_ANGLE["sw"]
+        position_index_to_adjust = 0
     elif positions == NSW_POSITIONS:
         pos_to_el = LO_POSITION_TO_ELEVATION_ANGLE["nsw"]
+        position_index_to_adjust = 9
     else:
         raise ValueError("Unknown positions for elevation angle mapping.")
 
@@ -369,6 +373,8 @@ def process_lo_angular_intensity(
             [pos_to_el[pos] for pos in dataset["inst_az"].data],
         )
     )
+    # add uncertainties to species list
+    species_list = species_list + [f"unc_{var}" for var in species_list]
     # Take the mean across elevation angles and restore the original dimension order
     dataset_converted = (
         dataset[species_list]
@@ -383,8 +389,41 @@ def process_lo_angular_intensity(
     dataset = dataset.assign_coords(
         spin_angle=("spin_sector", dataset["spin_sector"].data * 15.0 + 7.5)
     )
-
     dataset = dataset.drop_vars(species_list).merge(dataset_converted)
+    # Positions 0 and 10 only observe half of the 24 spins for each esa step.
+    # To account for this, we replicate the counts observed in position 0 and 10 for
+    # each esa step to either spin angles 0-11 or 12-23, depending on the pixel
+    # orientation (A/B). See section 11.2.2 of the CoDICE algorithm document
+    a_inds = np.array(
+        [pos for pos, orientation in PIXEL_ORIENTATIONS.items() if orientation == "A"]
+    )
+    b_inds = np.array(
+        [pos for pos, orientation in PIXEL_ORIENTATIONS.items() if orientation == "B"]
+    )
+
+    position_index = position_index_to_adjust
+    for species in species_list:
+        # Determine the correct spin indices based on the position
+        spin_sectors = dataset["spin_sector"].data
+        spin_inds_1 = np.where(spin_sectors >= 12)[0]
+        spin_inds_2 = np.where(spin_sectors < 12)[0]
+        # if position_index is 9, swap the spin indices
+        if position_index == 9:
+            spin_inds_1, spin_inds_2 = spin_inds_2, spin_inds_1
+
+        # Assign the values to the correct positions and spin sectors
+        dataset[species].values[
+            :, a_inds[:, np.newaxis], spin_inds_1, position_index
+        ] = dataset[species].values[
+            :, a_inds[:, np.newaxis], spin_inds_2, position_index
+        ]
+
+        dataset[species].values[
+            :, b_inds[:, np.newaxis], spin_inds_2, position_index
+        ] = dataset[species].values[
+            :, b_inds[:, np.newaxis], spin_inds_1, position_index
+        ]
+
     return dataset
 
 
@@ -777,7 +816,6 @@ def process_codice_l2(
     # This should get science files since ancillary or spice doesn't have data_type
     # as data level.
     file_path = dependencies.get_file_paths(descriptor=descriptor)[0]
-
     # Now form product name from descriptor
     descriptor = ScienceFilePath(file_path).descriptor
     dataset_name = f"imap_codice_l2_{descriptor}"
@@ -790,19 +828,23 @@ def process_codice_l2(
         "imap_codice_l2_lo-nsw-angular",
         "imap_codice_l2_lo-sw-angular",
     ]:
+        cdf_attrs = ImapCdfAttributes()
+        cdf_attrs.add_instrument_global_attrs("codice")
+
         l2_dataset = load_cdf(file_path).copy()
 
         geometric_factor_lookup = get_geometric_factor_lut(dependencies)
         efficiency_lookup = get_efficiency_lut(dependencies)
         geometric_factors = compute_geometric_factors(
             l2_dataset, geometric_factor_lookup
         )
+
         if dataset_name == "imap_codice_l2_lo-sw-species":
             # Filter the efficiency lookup table for solar wind efficiencies
             efficiencies = efficiency_lookup[efficiency_lookup["product"] == "sw"]
             # Calculate the pickup ion sunward solar wind intensities using equation
             # described in section 11.2.3 of algorithm document.
-            process_lo_species_intensity(
+            l2_dataset = process_lo_species_intensity(
                 l2_dataset,
                 LO_SW_PICKUP_ION_SPECIES_VARIABLE_NAMES,
                 geometric_factors,
@@ -811,25 +853,31 @@ def process_codice_l2(
             )
             # Calculate the sunward solar wind species intensities using equation
             # described in section 11.2.3 of algorithm document.
-            process_lo_species_intensity(
+            l2_dataset = process_lo_species_intensity(
                 l2_dataset,
-                LO_SW_SPECIES_VARIABLE_NAMES,
+                LO_SW_SOLAR_WIND_SPECIES_VARIABLE_NAMES,
                 geometric_factors,
                 efficiencies,
                 SOLAR_WIND_POSITIONS,
             )
+            l2_dataset.attrs.update(
+                cdf_attrs.get_global_attributes("imap_codice_l2_lo-sw-species")
+            )
         elif dataset_name == "imap_codice_l2_lo-nsw-species":
             # Filter the efficiency lookup table for non-solar wind efficiencies
             efficiencies = efficiency_lookup[efficiency_lookup["product"] == "nsw"]
             # Calculate the non-sunward species intensities using equation
             # described in section 11.2.3 of algorithm document.
-            process_lo_species_intensity(
+            l2_dataset = process_lo_species_intensity(
                 l2_dataset,
                 LO_NSW_SPECIES_VARIABLE_NAMES,
                 geometric_factors,
                 efficiencies,
                 NSW_POSITIONS,
             )
+            l2_dataset.attrs.update(
+                cdf_attrs.get_global_attributes("imap_codice_l2_lo-nsw-species")
+            )
         elif dataset_name == "imap_codice_l2_lo-sw-angular":
             efficiencies = efficiency_lookup[efficiency_lookup["product"] == "sw"]
             # Calculate the sunward solar wind angular intensities using equation
@@ -841,6 +889,9 @@ def process_codice_l2(
                 efficiencies,
                 SW_POSITIONS,
             )
+            l2_dataset.attrs.update(
+                cdf_attrs.get_global_attributes("imap_codice_l2_lo-sw-angular")
+            )
         if dataset_name == "imap_codice_l2_lo-nsw-angular":
             # Calculate the non sunward angular intensities
             efficiencies = efficiency_lookup[efficiency_lookup["product"] == "nsw"]
@@ -851,7 +902,9 @@ def process_codice_l2(
                 efficiencies,
                 NSW_POSITIONS,
             )
-
+            l2_dataset.attrs.update(
+                cdf_attrs.get_global_attributes("imap_codice_l2_lo-nsw-angular")
+            )
     if dataset_name in [
         "imap_codice_l2_hi-counters-singles",
         "imap_codice_l2_hi-counters-aggregated",

imap_processing/tests/codice/test_codice_l2.py

@@ -21,6 +21,7 @@
     process_lo_species_intensity,
 )
 from imap_processing.codice.constants import (
+    LO_NSW_ANGULAR_VARIABLE_NAMES,
     LO_SW_ANGULAR_VARIABLE_NAMES,
     LO_SW_SOLAR_WIND_SPECIES_VARIABLE_NAMES,
     SW_POSITIONS,
@@ -307,8 +308,12 @@ def test_process_lo_angular_intensity():
             .groupby("group")
             .sum()
         )
+        # Skip checking the first elevations. Those get reassigned and will be
+        # validated below.
         np.testing.assert_allclose(
-            l1b_val_data_processed[var].values, expected_intensity.values, rtol=1e-5
+            l1b_val_data_processed[var].values[:, :, :, 1:],
+            expected_intensity.values[:, :, :, 1:],
+            rtol=1e-5,
         )
     # Check coords
     np.testing.assert_allclose(l1b_val_data_processed["elevation_angle"], [0, 15, 30])
@@ -320,30 +325,100 @@ def test_process_lo_angular_intensity():
 def test_codice_l2_sw_species_intensity(processing_dependencies, mock_get_file_paths):
     sci_input = ScienceInput("imap_codice_l1b_lo-sw-species_20250814_v007.cdf")
     processing_dependencies.add(sci_input)
+
+    l2_val_data = (
+        imap_module_directory
+        / "tests"
+        / "codice"
+        / "data"
+        / "l2_validation"
+        / "imap_codice_l2_lo-sw-species_20250814_v007.cdf"
+    )
+    l2_val_data = load_cdf(l2_val_data)
     ds = process_codice_l2("lo-sw-species", processing_dependencies)
+    for variable in l2_val_data.data_vars:
+        processed_val = ds[variable].values
+        np.testing.assert_allclose(
+            processed_val,
+            l2_val_data[variable].values,
+            rtol=1e-5,
+            err_msg=f"Mismatch in variable '{variable}'",
+        )
     ds.attrs["Data_version"] = "001"
-    write_cdf(ds)
+    write_cdf(ds, istp=False)
 
 
 def test_codice_l2_nsw_species_intensity(processing_dependencies, mock_get_file_paths):
     sci_input = ScienceInput("imap_codice_l1b_lo-nsw-species_20250814_v007.cdf")
     processing_dependencies.add(sci_input)
+    l2_val_data = (
+        imap_module_directory
+        / "tests"
+        / "codice"
+        / "data"
+        / "l2_validation"
+        / "imap_codice_l2_lo-nsw-species_20250814_v007.cdf"
+    )
+    l2_val_data = load_cdf(l2_val_data)
     ds = process_codice_l2("lo-nsw-species", processing_dependencies)
+    for variable in l2_val_data.data_vars:
+        np.testing.assert_allclose(
+            ds[variable].values,
+            l2_val_data[variable].values,
+            rtol=1e-5,
+            err_msg=f"Mismatch in variable '{variable}'",
+        )
     ds.attrs["Data_version"] = "001"
-    write_cdf(ds)
+    write_cdf(ds, istp=False)
 
 
 def test_codice_l2_nsw_angular_intensity(processing_dependencies, mock_get_file_paths):
     sci_input = ScienceInput("imap_codice_l1b_lo-nsw-angular_20250814_v007.cdf")
     processing_dependencies.add(sci_input)
+    l2_val_data = (
+        imap_module_directory
+        / "tests"
+        / "codice"
+        / "data"
+        / "l2_validation"
+        / "imap_codice_l2_lo-nsw-angular_20250814_v007.cdf"
+    )
+    l2_val_data = load_cdf(l2_val_data)
     ds = process_codice_l2("lo-nsw-angular", processing_dependencies)
+    for variable in LO_NSW_ANGULAR_VARIABLE_NAMES:
+        np.testing.assert_allclose(
+            ds[variable].values,
+            l2_val_data[variable].values,
+            rtol=1e-5,
+            err_msg=f"Mismatch in variable '{variable}'",
+        )
     ds.attrs["Data_version"] = "001"
-    write_cdf(ds)
+    # TODO fix attrs
+    write_cdf(ds, istp=False)
 
 
 def test_codice_l2_sw_angular_intensity(processing_dependencies, mock_get_file_paths):
     sci_input = ScienceInput("imap_codice_l1b_lo-sw-angular_20250814_v007.cdf")
     processing_dependencies.add(sci_input)
+    l2_val_data = (
+        imap_module_directory
+        / "tests"
+        / "codice"
+        / "data"
+        / "l2_validation"
+        / "imap_codice_l2_lo-sw-angular_20250814_v007.cdf"
+    )
+    l2_val_data = load_cdf(l2_val_data)
     ds = process_codice_l2("lo-sw-angular", processing_dependencies)
+    for variable in LO_SW_ANGULAR_VARIABLE_NAMES:
+        np.testing.assert_allclose(
+            ds[variable].values,
+            l2_val_data[variable].values,
+            # TODO is 1e-4 ok?
+            rtol=1e-4,
+            err_msg=f"Mismatch in variable '{variable}'",
+        )
+
     ds.attrs["Data_version"] = "001"
-    write_cdf(ds)
+    # TODO fix attrs
+    write_cdf(ds, istp=False)