3165 bug hi l2 problems propegating error and uncertainty in maps

imap_processing/ena_maps/utils/corrections.py

@@ -957,7 +957,7 @@ def interpolate_map_flux_to_helio_frame(
             sys_err_helio = sys_err_sc * energy_ratio
 
         # Clamp negative fluxes to zero (can occur from linear extrapolation)
-        flux_helio = flux_helio.where(flux_helio >= 0, 0.0)
+        flux_helio = xr.where(flux_helio < 0, 0.0, flux_helio)
 
         # Set any location where the value is not finite to NaN (converts +/-inf to NaN)
         flux_helio = flux_helio.where(np.isfinite(flux_helio), np.nan)

imap_processing/hi/hi_l2.py

@@ -463,11 +463,13 @@ def calculate_ena_intensity(
     map_ds["ena_intensity_stat_uncert"] = (
         map_ds["ena_signal_rate_stat_unc"] / flux_conversion_divisor
     )
-    map_ds["ena_intensity_sys_err"] = (
-        np.sqrt(map_ds["bg_rate"] * map_ds["exposure_factor"])
-        / map_ds["exposure_factor"]
-        / flux_conversion_divisor
-    )
+
+    with np.errstate(divide="ignore", invalid="ignore"):
+        map_ds["ena_intensity_sys_err"] = (
+            np.sqrt(map_ds["bg_rate"] * map_ds["exposure_factor"])
+            / map_ds["exposure_factor"]
+            / flux_conversion_divisor
+        )
 
     # Combine calibration products using proper weighted averaging
     # as described in Hi Algorithm Document Section 3.1.2
@@ -540,11 +542,16 @@ def combine_calibration_products(
     map_ds["ena_intensity"] = combined_flux
     # Statistical uncertainty
     map_ds["ena_intensity_stat_uncert"] = np.sqrt(
-        1 / (1 / (map_ds["ena_intensity_stat_uncert"] ** 2)).sum(dim="calibration_prod")
+        1
+        / (1 / (map_ds["ena_intensity_stat_uncert"] ** 2)).sum(
+            dim="calibration_prod", skipna=True, min_count=1
+        )
     )
     # For systematic error, just do quadrature sum over the systematic error for
     # each calibration product.
-    map_ds["ena_intensity_sys_err"] = np.sqrt((sys_err**2).sum(dim="calibration_prod"))
+    map_ds["ena_intensity_sys_err"] = np.sqrt(
+        (sys_err**2).sum(dim="calibration_prod", skipna=True, min_count=1)
+    )
 
     return map_ds
 

imap_processing/tests/ena_maps/test_corrections.py

@@ -1382,6 +1382,33 @@ def test_infinite_values_converted_to_nan(self):
         assert not np.any(np.isinf(stat_unc))
         assert not np.any(np.isinf(sys_err))
 
+    def test_nan_input_propagation(self):
+        """Test that NaN inputs properly propagate to NaN outputs."""
+        map_ds, esa_energies, helio_energies = self.create_test_map_dataset(
+            n_energy=3, n_spatial=3
+        )
+
+        # Set some input values to NaN
+        map_ds["ena_intensity"].values[1, 0] = np.nan
+        map_ds["ena_intensity_stat_uncert"].values[1, 1] = np.nan
+        map_ds["ena_intensity_sys_err"].values[1, 2] = np.nan
+
+        result_ds = interpolate_map_flux_to_helio_frame(
+            map_ds, esa_energies, helio_energies, ["ena_intensity"]
+        )
+
+        # NaN in flux should propagate to flux, stat_uncert, and sys_err
+        # at that location
+        assert np.isnan(result_ds["ena_intensity"].values[1, 0])
+        assert np.isnan(result_ds["ena_intensity_stat_uncert"].values[1, 0])
+        assert np.isnan(result_ds["ena_intensity_sys_err"].values[1, 0])
+
+        # NaN in stat_uncert input should result in NaN stat_uncert output
+        assert np.isnan(result_ds["ena_intensity_stat_uncert"].values[1, 1])
+
+        # NaN in sys_err input should result in NaN sys_err output
+        assert np.isnan(result_ds["ena_intensity_sys_err"].values[1, 2])
+
     def test_multidimensional_spatial_coords(self):
         """Test that interpolation works with multi-dimensional spatial coordinates."""
 

imap_processing/tests/hi/test_hi_l2.py

@@ -733,6 +733,86 @@ def test_combine_calibration_products_edge_cases():
         assert "calibration_prod" not in result_ds[var].dims
 
 
+def test_combine_calibration_products_nan_handling():
+    """Test that combine_calibration_products handles NaN values correctly.
+
+    Tests the skipna=True, min_count=1 behavior:
+    - When one calibration product has NaN, the result should use the valid values
+    - When ALL calibration products have NaN at a position, the result should be NaN
+    """
+    # Create test dataset with 2 calibration products
+    coords = {
+        "epoch": 1,
+        "esa_energy_step": 2,
+        "calibration_prod": 2,
+        "longitude": 2,
+        "latitude": 1,
+    }
+
+    # Shape: (1, 2, 2, 2, 1) - epoch, energy, cal_prod, lon, lat
+    shape = (1, 2, 2, 2, 1)
+
+    # Create base arrays with valid values
+    intensity = np.full(shape, 100.0)
+    stat_uncert = np.full(shape, 10.0)
+    sys_err = np.full(shape, 5.0)
+    signal_rates = np.full(shape, 500.0)
+    bg_rate = np.full(shape, 20.0)
+    bg_rate_sys_err = np.full(shape, 2.0)
+    exposure_factor = np.full(shape, 1.0)
+
+    # Set NaN in one calibration product's uncertainty at position [0,0,0,0,0]
+    # The other calibration product is valid, so result should be finite
+    stat_uncert[0, 0, 0, 0, 0] = np.nan
+
+    # Set NaN in both calibration products at position [0,1,:,0,0]
+    # Since ALL products have NaN, result should be NaN (due to min_count=1)
+    stat_uncert[0, 1, 0, 0, 0] = np.nan
+    stat_uncert[0, 1, 1, 0, 0] = np.nan
+
+    # Set NaN in sys_err for one product at position [0,0,0,1,0]
+    sys_err[0, 0, 0, 1, 0] = np.nan
+
+    # Set NaN in sys_err for both products at position [0,1,:,1,0]
+    sys_err[0, 1, 0, 1, 0] = np.nan
+    sys_err[0, 1, 1, 1, 0] = np.nan
+
+    dim_names = list(coords.keys())
+    test_ds = xr.Dataset(
+        {
+            "ena_intensity": xr.DataArray(intensity, dims=dim_names),
+            "ena_intensity_stat_uncert": xr.DataArray(stat_uncert, dims=dim_names),
+            "ena_intensity_sys_err": xr.DataArray(sys_err, dims=dim_names),
+            "ena_signal_rates": xr.DataArray(signal_rates, dims=dim_names),
+            "bg_rate": xr.DataArray(bg_rate, dims=dim_names),
+            "bg_rate_sys_err": xr.DataArray(bg_rate_sys_err, dims=dim_names),
+            "exposure_factor": xr.DataArray(exposure_factor, dims=dim_names),
+        }
+    )
+
+    geom_factors = xr.DataArray(
+        np.ones((2, 2)), dims=["esa_energy_step", "calibration_prod"]
+    )
+
+    esa_energies = xr.DataArray(np.array([1.0, 2.0]), dims=["esa_energy_step"])
+
+    result_ds = combine_calibration_products(test_ds, geom_factors, esa_energies)
+
+    # Test 1: When one product has NaN stat_uncert, result should be finite
+    # (skipna=True uses the valid value from the other product)
+    assert np.isfinite(result_ds["ena_intensity_stat_uncert"].values[0, 0, 0, 0])
+
+    # Test 2: When ALL products have NaN stat_uncert, result should be NaN
+    # (min_count=1 ensures at least one valid value is needed)
+    assert np.isnan(result_ds["ena_intensity_stat_uncert"].values[0, 1, 0, 0])
+
+    # Test 3: When one product has NaN sys_err, result should be finite
+    assert np.isfinite(result_ds["ena_intensity_sys_err"].values[0, 0, 1, 0])
+
+    # Test 4: When ALL products have NaN sys_err, result should be NaN
+    assert np.isnan(result_ds["ena_intensity_sys_err"].values[0, 1, 1, 0])
+
+
 # =============================================================================
 # PSET PROCESSING TESTS
 # =============================================================================
@@ -1480,3 +1560,38 @@ def test_combine_maps_handles_nan_intensity(mock_sky_map_for_combine):
         expected_normal,
         decimal=5,
     )
+
+
+def test_combine_maps_handles_nan_uncertainties(mock_sky_map_for_combine):
+    """Test that combine_maps handles NaN values in uncertainties correctly.
+
+    When one map has NaN values in ena_intensity_stat_uncert, the weight for
+    that map should be zero, so the combined result uses only the other map's
+    value.
+    """
+    ram_map = mock_sky_map_for_combine()
+    anti_map = mock_sky_map_for_combine(intensity_offset=20)
+
+    # Set NaN in stat_uncert - should result in zero weight for that map
+    ram_stat_unc = ram_map.data_1d["ena_intensity_stat_uncert"].values.copy()
+    ram_stat_unc[0, 0, 0, 0] = np.nan
+    ram_map.data_1d["ena_intensity_stat_uncert"] = xr.DataArray(
+        ram_stat_unc, dims=ram_map.data_1d["ena_intensity_stat_uncert"].dims
+    )
+
+    sky_maps = {"ram": ram_map, "anti": anti_map}
+    result = combine_maps(sky_maps)
+
+    # Combined intensity at NaN uncertainty position should use only anti's
+    # value since ram has zero weight (due to NaN uncertainty)
+    assert np.isfinite(result.data_1d["ena_intensity"].values[0, 0, 0, 0])
+    # The result should be anti's intensity since ram has zero weight
+    expected_intensity = 70.0  # anti's intensity
+    np.testing.assert_almost_equal(
+        result.data_1d["ena_intensity"].values[0, 0, 0, 0],
+        expected_intensity,
+        decimal=5,
+    )
+
+    # Combined stat_uncert should also be finite (from anti's value)
+    assert np.isfinite(result.data_1d["ena_intensity_stat_uncert"].values[0, 0, 0, 0])