SWE: add uncertainty in L1B and L2

imap_processing/cdf/config/imap_swe_l1b_variable_attrs.yaml

@@ -110,6 +110,26 @@ science_data:
     Metadata field acq_duration is 17 uint. Max value of uint17 is 131071.
     Dividing max counts by acq_duration gave validmax
 
+counts_stat_uncert:
+  CATDESC: Counts statistical uncertainty
+  FIELDNAM: Counts Stats Uncertainty
+  DEPEND_0: epoch
+  DEPEND_1: esa_step
+  DEPEND_2: spin_sector
+  DEPEND_3: cem_id
+  LABL_PTR_1: esa_step_label
+  LABL_PTR_2: spin_sector_label
+  LABL_PTR_3: cem_id_label
+  DISPLAY_TYPE: spectrogram
+  FILLVAL: -1.0000000E+31
+  FORMAT: E19.5
+  UNITS: counts
+  VALIDMIN: 0.0
+  VALIDMAX: 1.7976931348623157e+308 # TODO: find actual value
+  VAR_TYPE: data
+  SCALETYP: linear
+  DICT_KEY: SPASE>Particle>ParticleType:Electron,ParticleQuantity:Counts,Qualifier:Uncertainty
+
 acquisition_time:
   CATDESC: Acquisition time organized by voltage step and spin sector
   DEPEND_0: epoch

imap_processing/cdf/config/imap_swe_l2_variable_attrs.yaml

@@ -221,6 +221,25 @@ flux:
   VAR_TYPE: data
   DICT_KEY: SPASE>Particle>ParticleType:Electron,ParticleQuantity:NumberFlux,Qualifier:Array
 
+counts_stat_uncert:
+  CATDESC: Counts statistical uncertainty
+  FIELDNAM: Counts Stats Uncertainty
+  DEPEND_0: epoch
+  DEPEND_1: esa_step
+  DEPEND_2: spin_sector
+  DEPEND_3: cem_id
+  LABL_PTR_1: esa_step_label
+  LABL_PTR_2: spin_sector_label
+  LABL_PTR_3: cem_id_label
+  DISPLAY_TYPE: spectrogram
+  FILLVAL: -1.0000000E+31
+  FORMAT: E19.5
+  UNITS: counts
+  VALIDMIN: 0.0
+  VALIDMAX: 1.7976931348623157e+308 # TODO: find actual value
+  VAR_TYPE: data
+  SCALETYP: linear
+  DICT_KEY: SPASE>Particle>ParticleType:Electron,ParticleQuantity:Counts,Qualifier:Uncertainty
 
 acquisition_time:
   CATDESC: Acquisition time organized by ESA step and spin sector

imap_processing/swe/l1b/swe_l1b.py

@@ -761,6 +761,10 @@ def swe_l1b_science(dependencies: ProcessingInputCollection) -> xr.Dataset:
 
     count_rate = convert_counts_to_rate(inflight_applied_count, acq_duration)
 
+    # Statistical uncertainty is sqrt(decompressed counts)
+    # TODO: Update this if SWE like to include deatime correciton. 
+    counts_stat_uncert = np.sqrt(populated_data["science_data"])
+
     # Store ESA energies of full cycle for L2 purposes.
     esa_energies = get_esa_energy_pattern(esa_lut_files[0])
     # Repeat energies to be in the same shape as the science data
@@ -894,6 +898,11 @@ def swe_l1b_science(dependencies: ProcessingInputCollection) -> xr.Dataset:
         dims=["epoch", "esa_step", "spin_sector", "cem_id"],
         attrs=cdf_attrs.get_variable_attributes("science_data"),
     )
+    science_dataset["counts_stat_uncert"] = xr.DataArray(
+        counts_stat_uncert,
+        dims=["epoch", "esa_step", "spin_sector", "cem_id"],
+        attrs=cdf_attrs.get_variable_attributes("counts_stat_uncert"),
+    )
     science_dataset["acquisition_time"] = xr.DataArray(
         acq_time,
         dims=["epoch", "esa_step", "spin_sector"],

imap_processing/swe/l2/swe_l2.py

@@ -132,7 +132,7 @@ def calculate_flux(
 
 
 def put_data_into_angle_bins(
-    data: np.ndarray, angle_bin_indices: npt.NDArray[np.int_]
+    data: np.ndarray, angle_bin_indices: npt.NDArray[np.int_], is_unc: bool = False
 ) -> npt.NDArray:
     """
     Put data in its angle bins.
@@ -142,10 +142,11 @@ def put_data_into_angle_bins(
     full cycle, it assigns data to the corresponding angle bin
     based on the provided indices.
 
-    Since multiple data points may fall into the same angle bin,
-    the function accumulates values and computes the average across
-    all 7 CEMs, ensuring that each bin contains a representative
-    mean value while maintaining the 7 CEM structure.
+    Since multiple data points can fall into the same angle bin,
+    this function assigns each data point to its bin and sums
+    the values within that bin. If the data is uncertainty data,
+    it computes the combined uncertainty for the bin; otherwise,
+    it calculates the averages.
 
     Parameters
     ----------
@@ -155,6 +156,10 @@ def put_data_into_angle_bins(
     angle_bin_indices : numpy.ndarray
         Indices of angle bins to put data in. Shape:
         (full_cycle_data, N_ESA_STEPS, N_ANGLE_BINS).
+    is_unc : bool
+        Whether data is uncertainty data or not. If yes, sqrt(sum(data)).
+        Otherwise, find mean of data.
+        Default to False.
 
     Returns
     -------
@@ -177,9 +182,16 @@ def put_data_into_angle_bins(
     time_indices = np.arange(data.shape[0])[:, None, None]
     energy_indices = np.arange(swe_constants.N_ESA_STEPS)[None, :, None]
 
-    # Use np.add.at() to accumulate values into bins
+    # Use np.add.at() to put values into bins and add values in the bins into one.
     np.add.at(binned_data, (time_indices, energy_indices, angle_bin_indices), data)
 
+    if is_unc:
+        # Calculate new uncertainty of each bin data(s).
+        # Per SWE instruction:
+        #   At L1B, 'data' is result from sqrt(counts). Now in L2, average
+        #   uncertainty data using this formula sqrt(sum(binned_data)).
+        return np.sqrt(binned_data)
+
     # Count occurrences in each bin to compute the mean.
     # Ensure float dtype for division
     bin_counts = np.zeros_like(binned_data, dtype=float)
@@ -418,5 +430,15 @@ def swe_l2(l1b_dataset: xr.Dataset) -> xr.Dataset:
         dims=["epoch", "energy", "inst_az", "inst_el"],
         attrs=cdf_attributes.get_variable_attributes("phase_space_density"),
     )
+    # Put uncertainty data into its spin angle bins and calculate new uncertainty
+    counts_stat_uncert_binned = put_data_into_angle_bins(
+        l1b_dataset["counts_stat_uncert"].data, spin_angle_bins_indices, is_unc=True
+    )
+    dataset["counts_stat_uncert"] = xr.DataArray(
+        counts_stat_uncert_binned,
+        name="counts_stat_uncert",
+        dims=["epoch", "energy", "inst_az", "inst_el"],
+        attrs=cdf_attributes.get_variable_attributes("counts_stat_uncert"),
+    )
 
     return dataset