From 63416d0a2cd01931c74e840a76b41b4ecad1f34c Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Wed, 11 Mar 2026 07:40:26 -0600
Subject: [PATCH 01/10] counts at higher nside

---
 imap_processing/ultra/constants.py            |  5 +++++
 imap_processing/ultra/l1c/helio_pset.py       |  8 ++++++--
 imap_processing/ultra/l1c/spacecraft_pset.py  | 11 +++++++----
 imap_processing/ultra/utils/ultra_l1_utils.py |  8 +++++++-
 4 files changed, 25 insertions(+), 7 deletions(-)

diff --git a/imap_processing/ultra/constants.py b/imap_processing/ultra/constants.py
index d66b456af..81a2f97e8 100644
--- a/imap_processing/ultra/constants.py
+++ b/imap_processing/ultra/constants.py
@@ -92,6 +92,11 @@ class UltraConstants:
         300.0,
         1e5,
     ]
+    # Counts at l1c are sampled at a finer resolution to help maintain the
+    # pointing accuracy for each event.  Since count maps are a binned integral
+    # quantity, they necessarily require a non-spun approach per Pointing, unlike
+    # exposure time and sensitivities
+    L1C_COUNTS_NSIDE = 128
 
     PSET_ENERGY_BIN_EDGES: ClassVar[list] = [
         3.0,
diff --git a/imap_processing/ultra/l1c/helio_pset.py b/imap_processing/ultra/l1c/helio_pset.py
index 58e9e02a3..3aab888e2 100644
--- a/imap_processing/ultra/l1c/helio_pset.py
+++ b/imap_processing/ultra/l1c/helio_pset.py
@@ -2,6 +2,7 @@
 
 import logging
 
+import astropy_healpix.healpy as hp
 import numpy as np
 import xarray as xr
 
@@ -160,12 +161,14 @@ def calculate_helio_pset(
         vhat_dps_helio,
         species_dataset["energy_heliosphere"].values,
         intervals,
-        nside=nside,
+        nside=UltraConstants.L1C_COUNTS_NSIDE,
     )
     helio_pset_quality_flags = np.full(
         n_pix, ImapPSETUltraFlags.NONE.value, dtype=np.uint16
     )
-    healpix = np.arange(n_pix)
+    counts_healpix = np.arange(n_pix)
+    # Determine nside for non "counts" variables from the lookup table
+    healpix = hp.nside2npix(nside)
 
     logger.info("Calculating spacecraft exposure times with deadtime correction.")
     exposure_time, deadtime_ratios = get_spacecraft_exposure_times(
@@ -236,6 +239,7 @@ def calculate_helio_pset(
     pset_dict["background_rates"] = background_rates[np.newaxis, ...]
     pset_dict["exposure_factor"] = exposure_time[np.newaxis, ...]
     pset_dict["pixel_index"] = healpix
+    pset_dict["counts_pixel_index"] = counts_healpix
     pset_dict["energy_bin_delta"] = np.diff(intervals, axis=1).squeeze()[
         np.newaxis, ...
     ]
diff --git a/imap_processing/ultra/l1c/spacecraft_pset.py b/imap_processing/ultra/l1c/spacecraft_pset.py
index 59b1b5706..5b9df3af0 100644
--- a/imap_processing/ultra/l1c/spacecraft_pset.py
+++ b/imap_processing/ultra/l1c/spacecraft_pset.py
@@ -139,15 +139,17 @@ def calculate_spacecraft_pset(
             reject_scattering,
         )
     )
-    # Determine nside from the lookup table
-    nside = hp.npix2nside(for_indices_by_spin_phase.sizes["pixel"])
     counts, latitude, longitude, n_pix = get_spacecraft_histogram(
         vhat_dps_spacecraft,
         species_dataset["energy_spacecraft"].values,
         intervals,
-        nside=nside,
+        nside=UltraConstants.L1C_COUNTS_NSIDE,
     )
-    healpix = np.arange(n_pix)
+    counts_healpix = np.arange(n_pix)
+    # Determine nside for non "counts" variables from the lookup table
+    nside = hp.npix2nside(for_indices_by_spin_phase.sizes["pixel"])
+    healpix = hp.nside2npix(nside)
+
     # Get the start and stop times of the pointing period
     repoint_id = species_dataset.attrs.get("Repointing", None)
     if repoint_id is None:
@@ -225,6 +227,7 @@ def calculate_spacecraft_pset(
     pset_dict["background_rates"] = background_rates[np.newaxis, ...]
     pset_dict["exposure_factor"] = exposure_pointing[np.newaxis, ...]
     pset_dict["pixel_index"] = healpix
+    pset_dict["counts_pixel_index"] = counts_healpix
     pset_dict["energy_bin_delta"] = np.diff(intervals, axis=1).squeeze()[
         np.newaxis, ...
     ]
diff --git a/imap_processing/ultra/utils/ultra_l1_utils.py b/imap_processing/ultra/utils/ultra_l1_utils.py
index 7efc4d8a5..12f06faea 100644
--- a/imap_processing/ultra/utils/ultra_l1_utils.py
+++ b/imap_processing/ultra/utils/ultra_l1_utils.py
@@ -46,6 +46,7 @@ def create_dataset(  # noqa: PLR0912
         coords = {
             "epoch": data_dict["epoch"],
             "pixel_index": data_dict["pixel_index"],
+            "counts_pixel_index": data_dict["counts_pixel_index"],
             "energy_bin_geometric_mean": data_dict["energy_bin_geometric_mean"],
             "spin_phase_step": data_dict["spin_phase_step"],
         }
@@ -155,7 +156,6 @@ def create_dataset(  # noqa: PLR0912
                 attrs=cdf_manager.get_variable_attributes(key, check_schema=False),
             )
         elif key in {
-            "counts",
             "background_rates",
             "exposure_factor",
             "helio_exposure_factor",
@@ -165,6 +165,12 @@ def create_dataset(  # noqa: PLR0912
                 dims=["epoch", "energy_bin_geometric_mean", "pixel_index"],
                 attrs=cdf_manager.get_variable_attributes(key, check_schema=False),
             )
+        elif key in {"counts"}:
+            dataset[key] = xr.DataArray(
+                data,
+                dims=["epoch", "energy_bin_geometric_mean", "counts_pixel_index"],
+                attrs=cdf_manager.get_variable_attributes(key, check_schema=False),
+            )
         elif key in {
             "geometric_function",
             "scatter_theta",

From 68ce19bee015e6456b88f2563a52d8338367ce49 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Wed, 11 Mar 2026 10:17:23 -0600
Subject: [PATCH 02/10] variable attrsr

---
 imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml b/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
index acfbe6fc1..374d4dbf9 100644
--- a/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
+++ b/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
@@ -137,7 +137,7 @@ counts:
   <<: *default
   CATDESC: Counts for a spin.
   DEPEND_1: energy_bin_geometric_mean
-  DEPEND_2: pixel_index
+  DEPEND_2: counts_pixel_index
   FIELDNAM: counts
   LABLAXIS: counts
   # TODO: come back to format

From b2e10e289e5e29b6d10711ab24e87db417b13624 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Mon, 23 Mar 2026 12:24:11 -0600
Subject: [PATCH 03/10] l2 tests

---
 imap_processing/ena_maps/ena_maps.py          | 65 +++++++++++++++++++
 imap_processing/ena_maps/utils/coordinates.py |  1 +
 .../tests/ena_maps/test_ena_maps.py           | 33 ++++++++--
 imap_processing/tests/ultra/mock_data.py      | 33 ++++++----
 4 files changed, 114 insertions(+), 18 deletions(-)

diff --git a/imap_processing/ena_maps/ena_maps.py b/imap_processing/ena_maps/ena_maps.py
index ed7b67523..23c2249a9 100644
--- a/imap_processing/ena_maps/ena_maps.py
+++ b/imap_processing/ena_maps/ena_maps.py
@@ -558,6 +558,8 @@ def __init__(
             np.stack((azimuth_pixel_center, elevation_pixel_center), axis=-1),
             dims=[CoordNames.GENERIC_PIXEL.value, CoordNames.AZ_EL_VECTOR.value],
         )
+        # downsample counts variable to match the nside of the pointing set
+        self.downsample_counts()
 
     @property
     def num_points(self) -> int:
@@ -600,6 +602,69 @@ def __repr__(self) -> str:
             f"num_points={self.num_points})"
         )
 
+    def downsample_counts(self) -> None:
+        """
+        Downsample the counts variable to match the pset nside.
+
+        Counts at l1c are sampled at a finer resolution to help maintain the
+        pointing accuracy for each event. Since count maps are a binned integral
+        quantity, they necessarily require a non-spun approach per Pointing, unlike
+        exposure time and sensitivities. We need to downsample the counts from the
+        nside of the input pset counts variable (e.g. 128) to the nside of the pset.
+        """
+        pset_data = self.data
+        counts_n_pix = pset_data.sizes["counts_pixel_index"]
+        if counts_n_pix != hp.nside2npix(self.nside):
+            # Raise an error if the nside the counts were sampled at is lower than the
+            # nside of the output map. We never want counts to be upsampled.
+            if counts_n_pix < hp.nside2npix(self.nside):
+                raise ValueError(
+                    f"Counts in the input PSET are sampled at nside "
+                    f"{hp.npix2nside(counts_n_pix)}, and the pset is {self.nside}. "
+                    f"This would require upsampling the counts, which we do not want."
+                )
+            counts_nside = hp.npix2nside(counts_n_pix)
+            pset_n_pix = hp.nside2npix(self.nside)
+            n_energy_bins = pset_data.sizes["energy_bin_geometric_mean"]
+            downsampled_counts = np.zeros((n_energy_bins, pset_n_pix))
+            order_diff = int(np.log2(counts_nside // self.nside))
+            for i in range(n_energy_bins):
+                counts = pset_data["counts"].values[0, i]
+                # Convert to nested ordering if necessary. In nested ordering, the
+                # pixels that need to be binned together to go from the counts nside to
+                # the pset nside are contiguous in the array.
+                counts_n = (
+                    counts
+                    if self.nested
+                    else counts[hp.ring2nest(counts_nside, np.arange(counts_n_pix))]
+                )
+                # reshape the counts by the amount pixels to bin together which is
+                # 4**order_diff because each step in order multiplies the pixel count
+                # by 4
+                binned_counts_n = counts_n.reshape((-1, 4**order_diff)).sum(axis=1)
+                # convert back to ring ordering if necessary and store in the
+                # downsampled counts array
+                downsampled_counts[i] = (
+                    binned_counts_n
+                    if self.nested
+                    else binned_counts_n[
+                        hp.nest2ring(self.nside, np.arange(pset_n_pix))
+                    ]
+                )
+
+            self.data["counts"] = xr.DataArray(
+                downsampled_counts[np.newaxis, :, :],
+                dims=(
+                    *self.data["counts"].dims[:-1],
+                    CoordNames.HEALPIX_INDEX.value,
+                ),
+            )
+        else:
+            # Update the counts variable with the correct dims
+            self.data["counts"] = self.data["counts"].rename(
+                {CoordNames.COUNTS_HEALPIX_INDEX.value: CoordNames.HEALPIX_INDEX.value}
+            )
+
 
 class LoHiBasePointingSet(PointingSet):
     """
diff --git a/imap_processing/ena_maps/utils/coordinates.py b/imap_processing/ena_maps/utils/coordinates.py
index 079e27a4c..4ad217f2f 100644
--- a/imap_processing/ena_maps/utils/coordinates.py
+++ b/imap_processing/ena_maps/utils/coordinates.py
@@ -12,6 +12,7 @@ class CoordNames(Enum):
     ENERGY_ULTRA_L1C = "energy_bin_geometric_mean"
     ENERGY_L2 = "energy"
     HEALPIX_INDEX = "pixel_index"
+    COUNTS_HEALPIX_INDEX = "counts_pixel_index"
 
     # The names of the az/el angular coordinates may differ between L1C and L2 data
     AZIMUTH_L1C = "longitude"
diff --git a/imap_processing/tests/ena_maps/test_ena_maps.py b/imap_processing/tests/ena_maps/test_ena_maps.py
index 617aa1035..42c40553a 100644
--- a/imap_processing/tests/ena_maps/test_ena_maps.py
+++ b/imap_processing/tests/ena_maps/test_ena_maps.py
@@ -8,6 +8,7 @@
 from copy import deepcopy
 from pathlib import Path
 from unittest import mock
+from unittest.mock import patch
 
 import astropy_healpix.healpy as hp
 import numpy as np
@@ -101,10 +102,14 @@ def test_init_cdf(
 
         cdf_filepath = write_cdf(ultra_pset, istp=False)
 
-        ultra_pset_from_dataset = ena_maps.UltraPointingSet(ultra_pset)
-
-        ultra_pset_from_str = ena_maps.UltraPointingSet(cdf_filepath)
-        ultra_pset_from_path = ena_maps.UltraPointingSet(Path(cdf_filepath))
+        # Mock the downsample_counts method to avoid dimension bugs
+        # since this is a dummy cdf, the dimensions are not present.
+        with patch.object(
+            ena_maps.UltraPointingSet, "downsample_counts", lambda self: None
+        ):
+            ultra_pset_from_str = ena_maps.UltraPointingSet(cdf_filepath)
+            ultra_pset_from_path = ena_maps.UltraPointingSet(Path(cdf_filepath))
+            ultra_pset_from_dataset = ena_maps.UltraPointingSet(ultra_pset)
 
         np.testing.assert_allclose(
             ultra_pset_from_dataset.data["counts"].values,
@@ -118,7 +123,6 @@ def test_init_cdf(
             rtol=1e-6,
         )
 
-    @pytest.mark.usefixtures("_setup_ultra_l1c_pset_products")
     @pytest.mark.usefixtures("_setup_ultra_l1c_pset_products")
     def test_different_spacing_raises_error(self):
         """Test that different spaced az/el from the L1C dataset raises ValueError"""
@@ -135,6 +139,25 @@ def test_different_spacing_raises_error(self):
                 spice_reference_frame=geometry.SpiceFrame.IMAP_DPS,
             )
 
+    def test_downsample_counts(self):
+        ultra_pset = self.l1c_pset_products[0]
+
+        # First check that counts are at a finer resolution than the spatial grid
+        # Verify counts start at a finer resolution than the spatial grid
+        counts_nside_before = hp.npix2nside(ultra_pset["counts"].shape[-1])
+        assert counts_nside_before != ultra_pset["exposure_factor"].shape[-1]
+        pset = ena_maps.UltraPointingSet(ultra_pset)
+
+        # Verify counts are now at the same resolution as the spatial grid
+        counts_nside_after = hp.npix2nside(pset.data["counts"].shape[-1])
+        assert counts_nside_after == pset.nside
+
+        # Verify counts are conserved
+        # (sum before == sum after, since power=-2 keeps sum invariant)
+        np.testing.assert_allclose(
+            pset.data["counts"].values.sum(), ultra_pset["counts"].values.sum()
+        )
+
 
 @pytest.fixture(scope="module")
 def hi_pset_cdf_path(imap_tests_path):
diff --git a/imap_processing/tests/ultra/mock_data.py b/imap_processing/tests/ultra/mock_data.py
index 4a4bae451..35cab46e5 100644
--- a/imap_processing/tests/ultra/mock_data.py
+++ b/imap_processing/tests/ultra/mock_data.py
@@ -193,7 +193,8 @@ def get_binomial_counts(distance_scaling, lat_bin, central_lat_bin):
 # TODO: Add ability to mock with/without energy dim to exposure_factor
 # The Helio frame L1C will have the energy dimension, but the spacecraft frame will not.
 def mock_l1c_pset_product_healpix(
-    nside: int = DEFAULT_HEALPIX_NSIDE_L1C,
+    nside: int = 32,
+    counts_nside: int = 128,
     stripe_center_lat: int = 0,
     width_scale: float = 10.0,
     counts_scaling_params: tuple[int, float] = (100, 0.01),
@@ -264,19 +265,19 @@ def mock_l1c_pset_product_healpix(
     energy_bin_delta = np.diff(energy_intervals, axis=1).squeeze()
     num_energy_bins = len(energy_bin_midpoints)
     npix = hp.nside2npix(nside)
-    counts = np.zeros(npix)
-    exposure_time = np.zeros(npix)
-
+    counts_npix = hp.nside2npix(counts_nside)
+    # counts are binned at a higher resolution than the other variables. See L1c
+    # code for more details.
     # Get latitude for each healpix pixel
-    pix_indices = np.arange(npix)
-    lon_pix, lat_pix = hp.pix2ang(nside, pix_indices, lonlat=True)
-
-    counts = np.zeros(shape=(num_energy_bins, npix))
+    counts_pix_indices = np.arange(counts_npix)
+    counts_lon_pix, counts_lat_pix = hp.pix2ang(
+        counts_nside, counts_pix_indices, lonlat=True
+    )
 
     # Calculate probability based on distance from target latitude
-    lat_diff = np.abs(lat_pix - stripe_center_lat)
+    counts_lat_diff = np.abs(counts_lat_pix - stripe_center_lat)
     prob_scaling_factor = counts_scaling_params[1] * np.exp(
-        -(lat_diff**2) / (2 * width_scale**2)
+        -(counts_lat_diff**2) / (2 * width_scale**2)
     )
     # Generate counts using binomial distribution
     rng = np.random.default_rng(seed=42)
@@ -286,6 +287,12 @@ def mock_l1c_pset_product_healpix(
             for _ in range(num_energy_bins)
         ]
     )
+    # Get latitude for each healpix pixel
+    pix_indices = np.arange(npix)
+    lon_pix, lat_pix = hp.pix2ang(nside, pix_indices, lonlat=True)
+
+    # Calculate probability based on distance from target latitude
+    lat_diff = np.abs(lat_pix - stripe_center_lat)
 
     # Generate exposure times using gaussian distribution, but wider
     prob_scaling_factor_exptime = counts_scaling_params[1] * np.exp(
@@ -317,7 +324,7 @@ def mock_l1c_pset_product_healpix(
     counts = counts.astype(int)
     # add an epoch dimension
     counts = np.expand_dims(counts, axis=0)
-    ones_ds = np.ones_like(counts)[0]  # pointing independent
+    ones_ds = np.ones((num_energy_bins, npix))  # pointing independent
     sensitivity = ones_ds
     geometric_function = ones_ds
     efficiency = ones_ds
@@ -338,7 +345,7 @@ def mock_l1c_pset_product_healpix(
                 [
                     CoordNames.TIME.value,
                     CoordNames.ENERGY_ULTRA_L1C.value,
-                    CoordNames.HEALPIX_INDEX.value,
+                    CoordNames.COUNTS_HEALPIX_INDEX.value,
                 ],
                 counts,
             ),
@@ -348,7 +355,7 @@ def mock_l1c_pset_product_healpix(
                     CoordNames.ENERGY_ULTRA_L1C.value,
                     CoordNames.HEALPIX_INDEX.value,
                 ],
-                np.full_like(counts, 0.05, dtype=float),
+                np.full((1, num_energy_bins, npix), 0.05, dtype=float),
             ),
             "exposure_factor": (
                 exposure_dims,  # special case: optionally energy dependent exposure

From d9a420f2054031e544b2d7f3d0472c84bd9bdd97 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Mon, 23 Mar 2026 12:34:12 -0600
Subject: [PATCH 04/10] comments

---
 .../cdf/config/imap_ultra_l1c_variable_attrs.yaml         | 3 +++
 imap_processing/tests/ena_maps/test_ena_maps.py           | 8 +++-----
 imap_processing/ultra/constants.py                        | 5 +----
 3 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml b/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
index 374d4dbf9..287f691ba 100644
--- a/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
+++ b/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
@@ -136,6 +136,9 @@ scatter_threshold:
 counts:
   <<: *default
   CATDESC: Counts for a spin.
+  VAR_NOTES: Counts healpix maps are sampled at finer resolution to help maintain the pointing accuracy for each event.
+    Since count maps are a binned integral quantity, they necessarily require a non-spun approach per Pointing, unlike
+    exposure time and sensitivities.
   DEPEND_1: energy_bin_geometric_mean
   DEPEND_2: counts_pixel_index
   FIELDNAM: counts
diff --git a/imap_processing/tests/ena_maps/test_ena_maps.py b/imap_processing/tests/ena_maps/test_ena_maps.py
index 42c40553a..6b138450d 100644
--- a/imap_processing/tests/ena_maps/test_ena_maps.py
+++ b/imap_processing/tests/ena_maps/test_ena_maps.py
@@ -142,18 +142,16 @@ def test_different_spacing_raises_error(self):
     def test_downsample_counts(self):
         ultra_pset = self.l1c_pset_products[0]
 
-        # First check that counts are at a finer resolution than the spatial grid
-        # Verify counts start at a finer resolution than the spatial grid
+        # First check that counts are at a finer resolution than exposure factor
         counts_nside_before = hp.npix2nside(ultra_pset["counts"].shape[-1])
         assert counts_nside_before != ultra_pset["exposure_factor"].shape[-1]
         pset = ena_maps.UltraPointingSet(ultra_pset)
 
-        # Verify counts are now at the same resolution as the spatial grid
+        # Verify counts are now at the same resolution as pset
         counts_nside_after = hp.npix2nside(pset.data["counts"].shape[-1])
         assert counts_nside_after == pset.nside
 
-        # Verify counts are conserved
-        # (sum before == sum after, since power=-2 keeps sum invariant)
+        # Verify counts are the same after downsampling.
         np.testing.assert_allclose(
             pset.data["counts"].values.sum(), ultra_pset["counts"].values.sum()
         )
diff --git a/imap_processing/ultra/constants.py b/imap_processing/ultra/constants.py
index 81a2f97e8..858dba983 100644
--- a/imap_processing/ultra/constants.py
+++ b/imap_processing/ultra/constants.py
@@ -92,10 +92,7 @@ class UltraConstants:
         300.0,
         1e5,
     ]
-    # Counts at l1c are sampled at a finer resolution to help maintain the
-    # pointing accuracy for each event.  Since count maps are a binned integral
-    # quantity, they necessarily require a non-spun approach per Pointing, unlike
-    # exposure time and sensitivities
+    # Counts at l1c are sampled at a finer resolution.
     L1C_COUNTS_NSIDE = 128
 
     PSET_ENERGY_BIN_EDGES: ClassVar[list] = [

From 12e335356806361f94fa7aeece780ce468d97533 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Mon, 23 Mar 2026 12:43:29 -0600
Subject: [PATCH 05/10] fix counts npix

---
 imap_processing/ultra/l1c/helio_pset.py       | 7 ++++---
 imap_processing/ultra/l1c/spacecraft_pset.py  | 9 +++++----
 imap_processing/ultra/utils/ultra_l1_utils.py | 1 +
 3 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/imap_processing/ultra/l1c/helio_pset.py b/imap_processing/ultra/l1c/helio_pset.py
index 3aab888e2..0d5f20194 100644
--- a/imap_processing/ultra/l1c/helio_pset.py
+++ b/imap_processing/ultra/l1c/helio_pset.py
@@ -157,18 +157,19 @@ def calculate_helio_pset(
         )
     )
 
-    counts, latitude, longitude, n_pix = get_spacecraft_histogram(
+    counts, latitude, longitude, counts_n_pix = get_spacecraft_histogram(
         vhat_dps_helio,
         species_dataset["energy_heliosphere"].values,
         intervals,
         nside=UltraConstants.L1C_COUNTS_NSIDE,
     )
+    n_pix = hp.nside2npix(nside)
     helio_pset_quality_flags = np.full(
         n_pix, ImapPSETUltraFlags.NONE.value, dtype=np.uint16
     )
-    counts_healpix = np.arange(n_pix)
+    counts_healpix = np.arange(counts_n_pix)
     # Determine nside for non "counts" variables from the lookup table
-    healpix = hp.nside2npix(nside)
+    healpix = np.arange(n_pix)
 
     logger.info("Calculating spacecraft exposure times with deadtime correction.")
     exposure_time, deadtime_ratios = get_spacecraft_exposure_times(
diff --git a/imap_processing/ultra/l1c/spacecraft_pset.py b/imap_processing/ultra/l1c/spacecraft_pset.py
index 5b9df3af0..4cb781d88 100644
--- a/imap_processing/ultra/l1c/spacecraft_pset.py
+++ b/imap_processing/ultra/l1c/spacecraft_pset.py
@@ -139,16 +139,17 @@ def calculate_spacecraft_pset(
             reject_scattering,
         )
     )
-    counts, latitude, longitude, n_pix = get_spacecraft_histogram(
+    counts, latitude, longitude, counts_n_pix = get_spacecraft_histogram(
         vhat_dps_spacecraft,
         species_dataset["energy_spacecraft"].values,
         intervals,
         nside=UltraConstants.L1C_COUNTS_NSIDE,
     )
-    counts_healpix = np.arange(n_pix)
+    counts_healpix = np.arange(counts_n_pix)
     # Determine nside for non "counts" variables from the lookup table
-    nside = hp.npix2nside(for_indices_by_spin_phase.sizes["pixel"])
-    healpix = hp.nside2npix(nside)
+    n_pix = for_indices_by_spin_phase.sizes["pixel"]
+    nside = hp.npix2nside(n_pix)
+    healpix = np.arange(n_pix)
 
     # Get the start and stop times of the pointing period
     repoint_id = species_dataset.attrs.get("Repointing", None)
diff --git a/imap_processing/ultra/utils/ultra_l1_utils.py b/imap_processing/ultra/utils/ultra_l1_utils.py
index 12f06faea..76ec67451 100644
--- a/imap_processing/ultra/utils/ultra_l1_utils.py
+++ b/imap_processing/ultra/utils/ultra_l1_utils.py
@@ -108,6 +108,7 @@ def create_dataset(  # noqa: PLR0912
             "spin_number",
             "energy_bin_geometric_mean",
             "pixel_index",
+            "counts_pixel_index",
             "spin_phase_step",
         ]:
             # update attrs

From d6bc40a7f15200618df2c6bf0e5c57ca6c3d25fd Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Mon, 23 Mar 2026 16:23:11 -0600
Subject: [PATCH 06/10] pr comments

---
 .../config/imap_ultra_l1c_variable_attrs.yaml | 12 +++++
 imap_processing/ena_maps/ena_maps.py          | 50 ++++++++++---------
 .../ultra/unit/test_ultra_l1c_pset_bins.py    | 12 +----
 imap_processing/ultra/l1c/helio_pset.py       |  6 ++-
 imap_processing/ultra/l1c/spacecraft_pset.py  |  6 ++-
 .../ultra/l1c/ultra_l1c_pset_bins.py          | 12 +----
 6 files changed, 53 insertions(+), 45 deletions(-)

diff --git a/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml b/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
index 287f691ba..2db856954 100644
--- a/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
+++ b/imap_processing/cdf/config/imap_ultra_l1c_variable_attrs.yaml
@@ -247,6 +247,18 @@ pixel_index:
   VALIDMAX: *max_int
   VAR_TYPE: support_data
 
+counts_pixel_index:
+  FILLVAL: *min_int
+  CATDESC: Counts variable HEALPix pixel index.
+  DISPLAY_TYPE: no_plot
+  FIELDNAM: counts_pixel_index
+  FORMAT: I12
+  LABLAXIS: Counts Pixel Index
+  UNITS: " "
+  VALIDMIN: *min_int
+  VALIDMAX: *max_int
+  VAR_TYPE: support_data
+
 spin_phase_step:
   FILLVAL: *min_int
   CATDESC: Spin phase step index (1ms resolution).
diff --git a/imap_processing/ena_maps/ena_maps.py b/imap_processing/ena_maps/ena_maps.py
index 23c2249a9..395b61918 100644
--- a/imap_processing/ena_maps/ena_maps.py
+++ b/imap_processing/ena_maps/ena_maps.py
@@ -626,34 +626,38 @@ def downsample_counts(self) -> None:
             counts_nside = hp.npix2nside(counts_n_pix)
             pset_n_pix = hp.nside2npix(self.nside)
             n_energy_bins = pset_data.sizes["energy_bin_geometric_mean"]
-            downsampled_counts = np.zeros((n_energy_bins, pset_n_pix))
             order_diff = int(np.log2(counts_nside // self.nside))
-            for i in range(n_energy_bins):
-                counts = pset_data["counts"].values[0, i]
-                # Convert to nested ordering if necessary. In nested ordering, the
-                # pixels that need to be binned together to go from the counts nside to
-                # the pset nside are contiguous in the array.
-                counts_n = (
-                    counts
-                    if self.nested
-                    else counts[hp.ring2nest(counts_nside, np.arange(counts_n_pix))]
-                )
-                # reshape the counts by the amount pixels to bin together which is
-                # 4**order_diff because each step in order multiplies the pixel count
-                # by 4
-                binned_counts_n = counts_n.reshape((-1, 4**order_diff)).sum(axis=1)
+            counts = pset_data["counts"].values[
+                0
+            ]  # shape: (n_energy_bins, counts_n_pix)
+            # Get counts in nested ordering. In nested ordering, the
+            # pixels that need to be binned together to go from the counts nside to
+            # the pset nside are contiguous in the array.
+            if not self.nested:
+                counts_n = counts[
+                    :, hp.ring2nest(counts_nside, np.arange(counts_n_pix))
+                ]
+            else:
+                counts_n = counts
+
+            # reshape the counts by the amount pixels to bin together which is
+            # 4**order_diff because each step in order multiplies the pixel count
+            # by 4
+            # Shape: (n_energy_bins, pset_n_pix, 4**order_diff) ->
+            # (n_energy_bins, pset_n_pix)
+            binned_counts_n = counts_n.reshape(
+                (n_energy_bins, pset_n_pix, 4**order_diff)
+            ).sum(axis=-1)
+
+            if not self.nested:
                 # convert back to ring ordering if necessary and store in the
                 # downsampled counts array
-                downsampled_counts[i] = (
-                    binned_counts_n
-                    if self.nested
-                    else binned_counts_n[
-                        hp.nest2ring(self.nside, np.arange(pset_n_pix))
-                    ]
-                )
+                binned_counts_n = binned_counts_n[
+                    :, hp.nest2ring(self.nside, np.arange(pset_n_pix))
+                ]
 
             self.data["counts"] = xr.DataArray(
-                downsampled_counts[np.newaxis, :, :],
+                binned_counts_n[np.newaxis, :, :],
                 dims=(
                     *self.data["counts"].dims[:-1],
                     CoordNames.HEALPIX_INDEX.value,
diff --git a/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py b/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py
index 655af6306..e42cc33dd 100644
--- a/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py
+++ b/imap_processing/tests/ultra/unit/test_ultra_l1c_pset_bins.py
@@ -118,13 +118,9 @@ def test_get_spacecraft_histogram(test_data):
     energy_bin_edges, _, _ = build_energy_bins()
     subset_energy_bin_edges = energy_bin_edges[:3]
 
-    hist, latitude, longitude, n_pix = get_spacecraft_histogram(
-        v, energy, subset_energy_bin_edges, nside=1
-    )
+    hist, n_pix = get_spacecraft_histogram(v, energy, subset_energy_bin_edges, nside=1)
     assert hist.shape == (len(subset_energy_bin_edges), hp.nside2npix(1))
     assert n_pix == hp.nside2npix(1)
-    assert latitude.shape == (n_pix,)
-    assert longitude.shape == (n_pix,)
 
     # Spot check that 2 counts are in the second energy bin
     assert np.sum(hist[2, :]) == 2
@@ -135,14 +131,10 @@ def test_get_spacecraft_histogram(test_data):
         (2.5, 4.137),
         (3.385, 5.057),
     ]
-    hist, latitude, longitude, n_pix = get_spacecraft_histogram(
-        v, energy, overlapping_bins, nside=1
-    )
+    hist, n_pix = get_spacecraft_histogram(v, energy, overlapping_bins, nside=1)
     # Spot check that 3 counts are in the third energy bin
     assert np.sum(hist[2, :]) == 3
     assert n_pix == hp.nside2npix(1)
-    assert latitude.shape == (n_pix,)
-    assert longitude.shape == (n_pix,)
 
 
 def mock_imap_state(time, ref_frame):
diff --git a/imap_processing/ultra/l1c/helio_pset.py b/imap_processing/ultra/l1c/helio_pset.py
index 0d5f20194..96f6c2f56 100644
--- a/imap_processing/ultra/l1c/helio_pset.py
+++ b/imap_processing/ultra/l1c/helio_pset.py
@@ -157,7 +157,7 @@ def calculate_helio_pset(
         )
     )
 
-    counts, latitude, longitude, counts_n_pix = get_spacecraft_histogram(
+    counts, counts_n_pix = get_spacecraft_histogram(
         vhat_dps_helio,
         species_dataset["energy_heliosphere"].values,
         intervals,
@@ -171,6 +171,10 @@ def calculate_helio_pset(
     # Determine nside for non "counts" variables from the lookup table
     healpix = np.arange(n_pix)
 
+    # Calculate the corresponding longitude (az) latitude (el)
+    # center coordinates
+    longitude, latitude = hp.pix2ang(nside, healpix, lonlat=True)
+
     logger.info("Calculating spacecraft exposure times with deadtime correction.")
     exposure_time, deadtime_ratios = get_spacecraft_exposure_times(
         rates_dataset,
diff --git a/imap_processing/ultra/l1c/spacecraft_pset.py b/imap_processing/ultra/l1c/spacecraft_pset.py
index 4cb781d88..fcd96a577 100644
--- a/imap_processing/ultra/l1c/spacecraft_pset.py
+++ b/imap_processing/ultra/l1c/spacecraft_pset.py
@@ -139,7 +139,7 @@ def calculate_spacecraft_pset(
             reject_scattering,
         )
     )
-    counts, latitude, longitude, counts_n_pix = get_spacecraft_histogram(
+    counts, counts_n_pix = get_spacecraft_histogram(
         vhat_dps_spacecraft,
         species_dataset["energy_spacecraft"].values,
         intervals,
@@ -151,6 +151,10 @@ def calculate_spacecraft_pset(
     nside = hp.npix2nside(n_pix)
     healpix = np.arange(n_pix)
 
+    # Calculate the corresponding longitude (az) latitude (el)
+    # center coordinates
+    longitude, latitude = hp.pix2ang(nside, healpix, lonlat=True)
+
     # Get the start and stop times of the pointing period
     repoint_id = species_dataset.attrs.get("Repointing", None)
     if repoint_id is None:
diff --git a/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py b/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py
index 507ca75c5..69d925f2a 100644
--- a/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py
+++ b/imap_processing/ultra/l1c/ultra_l1c_pset_bins.py
@@ -79,7 +79,7 @@ def get_spacecraft_histogram(
     energy_bin_edges: list[tuple[float, float]],
     nside: int = 128,
     nested: bool = False,
-) -> tuple[NDArray, NDArray, NDArray, NDArray]:
+) -> tuple[NDArray, NDArray]:
     """
     Compute a 2D histogram of the particle data using HEALPix binning.
 
@@ -101,10 +101,6 @@ def get_spacecraft_histogram(
     -------
     hist : np.ndarray
         A 2D histogram array with shape (n_pix, n_energy_bins).
-    latitude : np.ndarray
-        Array of latitude values.
-    longitude : np.ndarray
-        Array of longitude values.
     n_pix : int
         Number of healpix pixels.
 
@@ -126,10 +122,6 @@ def get_spacecraft_histogram(
     # Compute number of HEALPix pixels that cover the sphere
     n_pix = hp.nside2npix(nside)
 
-    # Calculate the corresponding longitude (az) latitude (el)
-    # center coordinates
-    longitude, latitude = hp.pix2ang(nside, np.arange(n_pix), lonlat=True)
-
     # Get HEALPix pixel indices for each event
     # HEALPix expects latitude in [-90, 90] so we don't need to change elevation
     hpix_idx = hp.ang2pix(nside, az, el, nest=nested, lonlat=True)
@@ -143,7 +135,7 @@ def get_spacecraft_histogram(
         # Only count the events that fall within the energy bin
         hist[i, :] += np.bincount(hpix_idx[mask], minlength=n_pix).astype(np.float64)
 
-    return hist, latitude, longitude, n_pix
+    return hist, n_pix
 
 
 def get_spacecraft_count_rate_uncertainty(hist: NDArray, exposure: NDArray) -> NDArray:

From 4e54f6e380d33f237d9610f3c66878fda1e40d1d Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Mon, 23 Mar 2026 16:29:12 -0600
Subject: [PATCH 07/10] pr comments

---
 imap_processing/ena_maps/ena_maps.py            | 6 +++---
 imap_processing/tests/ena_maps/test_ena_maps.py | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/imap_processing/ena_maps/ena_maps.py b/imap_processing/ena_maps/ena_maps.py
index 395b61918..80900dca0 100644
--- a/imap_processing/ena_maps/ena_maps.py
+++ b/imap_processing/ena_maps/ena_maps.py
@@ -614,17 +614,17 @@ def downsample_counts(self) -> None:
         """
         pset_data = self.data
         counts_n_pix = pset_data.sizes["counts_pixel_index"]
-        if counts_n_pix != hp.nside2npix(self.nside):
+        pset_n_pix = hp.nside2npix(self.nside)
+        if counts_n_pix != pset_n_pix:
             # Raise an error if the nside the counts were sampled at is lower than the
             # nside of the output map. We never want counts to be upsampled.
-            if counts_n_pix < hp.nside2npix(self.nside):
+            if counts_n_pix < pset_n_pix:
                 raise ValueError(
                     f"Counts in the input PSET are sampled at nside "
                     f"{hp.npix2nside(counts_n_pix)}, and the pset is {self.nside}. "
                     f"This would require upsampling the counts, which we do not want."
                 )
             counts_nside = hp.npix2nside(counts_n_pix)
-            pset_n_pix = hp.nside2npix(self.nside)
             n_energy_bins = pset_data.sizes["energy_bin_geometric_mean"]
             order_diff = int(np.log2(counts_nside // self.nside))
             counts = pset_data["counts"].values[
diff --git a/imap_processing/tests/ena_maps/test_ena_maps.py b/imap_processing/tests/ena_maps/test_ena_maps.py
index 6b138450d..dddcb9d89 100644
--- a/imap_processing/tests/ena_maps/test_ena_maps.py
+++ b/imap_processing/tests/ena_maps/test_ena_maps.py
@@ -143,8 +143,8 @@ def test_downsample_counts(self):
         ultra_pset = self.l1c_pset_products[0]
 
         # First check that counts are at a finer resolution than exposure factor
-        counts_nside_before = hp.npix2nside(ultra_pset["counts"].shape[-1])
-        assert counts_nside_before != ultra_pset["exposure_factor"].shape[-1]
+        counts_npix_before = ultra_pset["counts"].shape[-1]
+        assert counts_npix_before != ultra_pset["exposure_factor"].shape[-1]
         pset = ena_maps.UltraPointingSet(ultra_pset)
 
         # Verify counts are now at the same resolution as pset

From c40b7cef618221915dbc2ad1c28f0ed0b8afe0b2 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Thu, 26 Mar 2026 09:48:45 -0600
Subject: [PATCH 08/10] trigger tests

---
 imap_processing/tests/ultra/unit/test_ultra_l2.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/imap_processing/tests/ultra/unit/test_ultra_l2.py b/imap_processing/tests/ultra/unit/test_ultra_l2.py
index f1bb18899..0aa84e058 100644
--- a/imap_processing/tests/ultra/unit/test_ultra_l2.py
+++ b/imap_processing/tests/ultra/unit/test_ultra_l2.py
@@ -139,7 +139,6 @@ def test_generate_ultra_healpix_skymap_single_pset(
             pset["energy_bin_delta"] = pset["energy_bin_delta"].expand_dims(
                 {CoordNames.TIME.value: pset["epoch"].values}
             )
-
         # Create the Healpix skymap in the desired frame.
         with furnish_kernels(self.required_kernel_names):
             hp_skymap, _ = ultra_l2.generate_ultra_healpix_skymap(

From d0594fb7cb7715df2848233d943ac8496ae7f122 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Thu, 26 Mar 2026 09:57:23 -0600
Subject: [PATCH 09/10] logging

---
 imap_processing/ena_maps/ena_maps.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/imap_processing/ena_maps/ena_maps.py b/imap_processing/ena_maps/ena_maps.py
index 80900dca0..2a7a9be82 100644
--- a/imap_processing/ena_maps/ena_maps.py
+++ b/imap_processing/ena_maps/ena_maps.py
@@ -663,6 +663,10 @@ def downsample_counts(self) -> None:
                     CoordNames.HEALPIX_INDEX.value,
                 ),
             )
+            logger.info(
+                f"Counts variable with nside = {counts_nside} downsampled to "
+                f"nside {self.nside}."
+            )
         else:
             # Update the counts variable with the correct dims
             self.data["counts"] = self.data["counts"].rename(

From 17734a0b22fcd0355d55c0292bf0396dc71766e7 Mon Sep 17 00:00:00 2001
From: Luisa <luisa.coakley@lasp.colorado.edu>
Date: Thu, 26 Mar 2026 10:21:44 -0600
Subject: [PATCH 10/10] reduce data in ultra tests

---
 imap_processing/tests/ena_maps/conftest.py      |  4 +++-
 .../tests/ultra/unit/test_ultra_l2.py           | 17 +++++++++++++----
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/imap_processing/tests/ena_maps/conftest.py b/imap_processing/tests/ena_maps/conftest.py
index 905f0b793..e6391c34a 100644
--- a/imap_processing/tests/ena_maps/conftest.py
+++ b/imap_processing/tests/ena_maps/conftest.py
@@ -12,12 +12,14 @@
 @pytest.fixture(scope="module")
 def ultra_l1c_pset_datasets():
     """Make fake L1C Ultra PSET products on a HEALPix tiling for testing"""
-    l1c_nside = 32
+    l1c_nside = 16
+    counts_nside = 32
     return {
         "nside": l1c_nside,
         "products": [
             mock_l1c_pset_product_healpix(
                 nside=l1c_nside,
+                counts_nside=counts_nside,
                 stripe_center_lat=mid_latitude,
                 width_scale=5,
                 counts_scaling_params=(50, 0.5),
diff --git a/imap_processing/tests/ultra/unit/test_ultra_l2.py b/imap_processing/tests/ultra/unit/test_ultra_l2.py
index 0aa84e058..c0ba0de7f 100644
--- a/imap_processing/tests/ultra/unit/test_ultra_l2.py
+++ b/imap_processing/tests/ultra/unit/test_ultra_l2.py
@@ -46,7 +46,8 @@ def _setup_spice_kernels_list(self, spice_test_data_path, furnish_kernels):
     def _mock_single_pset(self, _setup_spice_kernels_list, furnish_kernels):
         with furnish_kernels(self.required_kernel_names):
             self.ultra_pset = mock_l1c_pset_product_healpix(
-                nside=32,
+                nside=16,
+                counts_nside=32,
                 stripe_center_lat=0,
                 timestr="2025-05-15T12:00:00",
                 energy_dependent_exposure=True,
@@ -66,6 +67,7 @@ def _mock_multiple_psets(self, _setup_spice_kernels_list, furnish_kernels):
             self.ultra_psets = [
                 mock_l1c_pset_product_healpix(
                     nside=16,
+                    counts_nside=32,
                     stripe_center_lat=mid_latitude,
                     width_scale=5,
                     counts_scaling_params=(50, 0.5),
@@ -740,8 +742,13 @@ def test_ultra_l2_descriptor_hpmap(self, mock_data_dict, furnish_kernels):
     @pytest.mark.usefixtures("_mock_single_pset")
     def test_bin_pset_energy_bins_default(self):
         """Test binning with default bin sizes."""
-        # Avoid modifying the original pset
-        pset = self.ultra_pset.copy(deep=True)
+        pset = mock_l1c_pset_product_healpix(
+            nside=16,
+            counts_nside=16,
+            stripe_center_lat=0,
+            timestr="2025-05-15T12:00:00",
+            energy_dependent_exposure=True,
+        )
         # Set the values in the single input PSET
         # Create a mock array with known values to test binning
         # e.g., 0,0,0,0,1,1,1,1,2,2,2,2,...11,11
@@ -751,7 +758,9 @@ def test_bin_pset_energy_bins_default(self):
         mock_array = (
             np.ones_like(pset["exposure_factor"]) * mock_vals[np.newaxis, :, np.newaxis]
         )
-        pset["counts"].values = mock_array
+        pset["counts"].values = (
+            np.ones_like(pset["counts"]) * mock_vals[np.newaxis, :, np.newaxis]
+        )
         pset["exposure_factor"].values = mock_array
         pset["sensitivity"].values = mock_array[0]
         pset["geometric_function"].values = mock_array[0]