diff --git a/imap_processing/ena_maps/utils/corrections.py b/imap_processing/ena_maps/utils/corrections.py index 2e88ce961..5c4c9fd18 100644 --- a/imap_processing/ena_maps/utils/corrections.py +++ b/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) diff --git a/imap_processing/hi/hi_l2.py b/imap_processing/hi/hi_l2.py index ce880445b..a58a711e1 100644 --- a/imap_processing/hi/hi_l2.py +++ b/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 diff --git a/imap_processing/tests/ena_maps/test_corrections.py b/imap_processing/tests/ena_maps/test_corrections.py index e3dfe848d..6edcc0729 100644 --- a/imap_processing/tests/ena_maps/test_corrections.py +++ b/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.""" diff --git a/imap_processing/tests/hi/test_hi_l2.py b/imap_processing/tests/hi/test_hi_l2.py index e84a4e352..b1573929e 100644 --- a/imap_processing/tests/hi/test_hi_l2.py +++ b/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])