From 8f56b85251974c4487501e44a18ae7589b8e1e7f Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 29 Jul 2026 12:14:52 -0400 Subject: [PATCH 01/14] removed geometric factor as a constant, pulling it from ancillary --- imap_processing/lo/constants.py | 10 - imap_processing/lo/l2/lo_l2.py | 388 +++++++++++++++++++++++-- imap_processing/tests/lo/test_lo_l2.py | 85 ++++-- 3 files changed, 430 insertions(+), 53 deletions(-) diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index 40de8b1f9..0f173bb6c 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -81,16 +81,6 @@ class LoConstants: # The following are indexed by ESA level (0-indexed, ESA level = index + 1). # The 8th entry is the virtual E8 channel, unused by the map. - ESA_ENERGY: ClassVar[list[float]] = [ - 0.016, - 0.030, - 0.056, - 0.106, - 0.200, - 0.405, - 0.787, - 1.821, - ] GEO_FACTOR: ClassVar[list[float]] = [ 7.0e-5, 7.9e-5, diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 4821350a3..a60b8a2f9 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -1,10 +1,13 @@ """IMAP-Lo L2 data processing.""" import logging +from pathlib import Path import numpy as np +import pandas as pd import xarray as xr +from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.ena_maps.ena_maps import ( PointingSet, RectangularSkyMap, @@ -12,6 +15,7 @@ ) from imap_processing.ena_maps.utils.coordinates import CoordNames from imap_processing.ena_maps.utils.naming import MapDescriptor +from imap_processing.lo import lo_ancillary from imap_processing.lo.constants import LoConstants as c # noqa: N813 from imap_processing.lo.l1c.lo_l1c import compute_pointing_directions from imap_processing.spice.geometry import ( @@ -33,6 +37,9 @@ # or intensities are derived from them. ACCUMULATED_VARIABLES = ("ena_count", "exposure_factor", "bg_rate_exposure") +# The calibration ancillaries shipped with the package. +ANCILLARY_DATA_DIR = Path(__file__).parent.parent / "ancillary_data" + # ============================================================================= # MAIN ENTRY POINT # ============================================================================= @@ -84,6 +91,21 @@ def lo_l2( map_descriptor = MapDescriptor.from_string(descriptor) logger.info(f"Processing map for species: {map_descriptor.species}") + # Determine if corrections are needed and prepare oxygen data if required + ( + _sputtering_correction, + _bootstrap_correction, + _flux_correction, + _o_map_dataset, + _flux_factors, + _cg_correction, + ) = _prepare_corrections( + map_descriptor, descriptor, sci_dependencies, anc_dependencies + ) + + logger.info("Step 1: Loading ancillary data") + _efficiency_data = load_efficiency_data(anc_dependencies) + # The geometric factors in LoConstants are hydrogen only. if map_descriptor.species != "h": raise NotImplementedError( @@ -98,15 +120,20 @@ def lo_l2( pointings = _complete_pointings(sci_dependencies) logger.info(f"Building {descriptor} from {len(pointings)} pointings") - _initialize_accumulators(sky_map) - esa_mode = 0 + # Every pointing of a map is taken in the same ESA mode, so the last one + # sets the energies and passband widths the whole map is binned in. + esa_mode = _get_esa_mode(pointings[max(pointings)][2]) if pointings else 0 + energy = _esa_energy(map_descriptor.species, esa_mode) + + _initialize_accumulators(sky_map, energy) for repointing, (goodtimes, bgrates, histrates) in sorted(pointings.items()): logger.debug(f"Accumulating repoint{repointing:05d}") - esa_mode = _get_esa_mode(histrates) - _accumulate_pointing(goodtimes, bgrates, histrates, sky_map, map_descriptor) + _accumulate_pointing( + goodtimes, bgrates, histrates, sky_map, map_descriptor, energy + ) - variables = _calculate_rates_and_intensities(sky_map, esa_mode) + variables = _calculate_rates_and_intensities(sky_map, esa_mode, energy) dataset = _build_map_dataset(sky_map, variables, esa_mode) logger.info("IMAP-Lo L2 processing pipeline completed successfully") @@ -120,6 +147,224 @@ def lo_l2( ] +def _prepare_corrections( + map_descriptor: MapDescriptor, + descriptor: str, + sci_dependencies: dict, + anc_dependencies: list, +) -> tuple[bool, bool, bool, xr.Dataset | None, Path | None, bool]: + """ + Determine what corrections are needed and prepare oxygen dataset if required. + + This helper function encapsulates the logic for determining when sputtering + and bootstrap corrections should be applied, and handles the creation of + the oxygen dataset needed for sputtering corrections. + + Parameters + ---------- + map_descriptor : MapDescriptor + The parsed map descriptor containing species and data type information. + descriptor : str + The original descriptor string for creating the oxygen variant. + sci_dependencies : dict + Dictionary of datasets needed for L2 data product creation. + anc_dependencies : list + List of ancillary file paths. + + Returns + ------- + tuple[bool, bool, bool, xr.Dataset | None, Path | None, bool] + A tuple containing: + - sputtering_correction: Whether to apply sputtering corrections + - bootstrap_correction: Whether to apply bootstrap corrections + - flux_correction: Whether to apply flux corrections + - o_map_dataset: Oxygen dataset if needed, None otherwise + - flux_factors: Path to flux factors ancillary file if needed, + None otherwise + - cg_correction: Whether to apply CG correction to the dataset. + """ + # Default values - no corrections needed + sputtering_correction = False + bootstrap_correction = False + flux_correction = False + o_map_dataset = None + flux_factors: None | Path = None + + # Sputtering and bootstrap corrections are only applied to hydrogen ENA data + # Guard against recursion: don't process oxygen for oxygen maps + if ( + map_descriptor.species == "h" + and map_descriptor.principal_data == "ena" + and "-o-" not in descriptor + ): # Safety check to prevent infinite recursion + logger.info("Creating map for oxygen for sputtering corrections") + o_descriptor = descriptor.replace("-h-", "-o-") + o_map_dataset = lo_l2(sci_dependencies, anc_dependencies, o_descriptor)[0] + sputtering_correction = True + bootstrap_correction = True + + if "raw" not in map_descriptor.principal_data: + flux_correction = True + try: + flux_factors = next( + x for x in anc_dependencies if "esa-eta-fit-factors" in str(x) + ) + except StopIteration: + raise ValueError( + "No flux correction factor file found in ancillary dependencies" + ) from None + + cg_correction = True if map_descriptor.frame_descriptor == "hf" else False + + return ( + sputtering_correction, + bootstrap_correction, + flux_correction, + o_map_dataset, + flux_factors, + cg_correction, + ) + + +# ============================================================================= +# SETUP AND INITIALIZATION HELPERS +# ============================================================================= + + +def load_efficiency_data(anc_dependencies: list) -> pd.DataFrame: + """ + Load efficiency factor data from ancillary files. + + Parameters + ---------- + anc_dependencies : list + List of ancillary file paths to search for efficiency factor files. + + Returns + ------- + pd.DataFrame + Concatenated efficiency factor data from all matching files. + Returns empty DataFrame if no efficiency files found. + """ + efficiency_files = [ + anc_file + for anc_file in anc_dependencies + if "efficiency-factor" in str(anc_file) + ] + + if not efficiency_files: + logger.warning("No efficiency factor files found in ancillary dependencies") + return pd.DataFrame() + + logger.debug(f"Loading {len(efficiency_files)} efficiency factor files") + return pd.concat( + [lo_ancillary.read_ancillary_file(anc_file) for anc_file in efficiency_files], + ignore_index=True, + ) + + +def load_sputter_correction_data( + source_species: str, target_species: str +) -> pd.DataFrame: + """ + Load sputter correction factors from an ancillary file. + + Parameters + ---------- + source_species : str + The species doing the sputtering (e.g. "o" for oxygen). + target_species : str + The species being corrected (e.g. "h" for hydrogen). + + Returns + ------- + pd.DataFrame + Rows matching the given species pair, sorted ascending by esa_step, + with columns: source_species, target_species, esa_step, + sputter_factor, sputter_factor_uncertainty. + """ + sputter_files = sorted(ANCILLARY_DATA_DIR.glob("*sputter-correction-factors*")) + + if not sputter_files: + raise ValueError("No sputter correction files found") + + df = pd.concat( + [lo_ancillary.read_ancillary_file(f) for f in sputter_files], + ignore_index=True, + ) + mask = (df["source_species"] == source_species) & ( + df["target_species"] == target_species + ) + result = df[mask].sort_values("esa_step").reset_index(drop=True) + return result + + +def load_bootstrap_correction_data() -> pd.DataFrame: + """ + Load bootstrap correction factors from an ancillary file. + + Returns + ------- + pd.DataFrame + Bootstrap correction factors with columns: esa_step_i, esa_step_k, + bootstrap_factor. Indices are 1-based ESA step numbers where esa_step_k=8 + refers to the virtual E8 channel. + """ + bootstrap_files = sorted(ANCILLARY_DATA_DIR.glob("*bootstrap-correction-factors*")) + + if not bootstrap_files: + raise ValueError("No bootstrap correction factor files found") + + return pd.concat( + [lo_ancillary.read_ancillary_file(f) for f in bootstrap_files], + ignore_index=True, + ) + + +def finalize_dataset(dataset: xr.Dataset, descriptor: str) -> xr.Dataset: + """ + Add attributes and perform final dataset preparation. + + Parameters + ---------- + dataset : xr.Dataset + The dataset to finalize with attributes. + descriptor : str + The descriptor for this map dataset. + + Returns + ------- + xr.Dataset + The finalized dataset with all attributes added. + """ + # Initialize the attribute manager + attr_mgr = ImapCdfAttributes() + attr_mgr.add_instrument_global_attrs(instrument="lo") + attr_mgr.add_instrument_variable_attrs(instrument="enamaps", level="l2-common") + attr_mgr.add_instrument_variable_attrs(instrument="enamaps", level="l2-rectangular") + + # Add global and variable attributes + dataset.attrs.update(attr_mgr.get_global_attributes("imap_lo_l2_enamap")) + + # Our global attributes have placeholders for descriptor + # so iterate through here and fill that in with the map-specific descriptor + for key in ["Data_type", "Logical_source", "Logical_source_description"]: + dataset.attrs[key] = dataset.attrs[key].format(descriptor=descriptor) + for var in dataset.data_vars: + try: + dataset[var].attrs = attr_mgr.get_variable_attributes(var) + except KeyError: + # If no attributes found, try without schema validation + try: + dataset[var].attrs = attr_mgr.get_variable_attributes( + var, check_schema=False + ) + except KeyError: + logger.warning(f"No attributes found for variable {var}") + + return dataset + + # ============================================================================= # INPUT HANDLING # ============================================================================= @@ -166,18 +411,6 @@ def _complete_pointings( return pointings -def _esa_energy() -> np.ndarray: - """ - Get the energy of each ESA level the map is binned in. - - Returns - ------- - np.ndarray - The energy [keV] of each ESA level. - """ - return np.array(c.ESA_ENERGY[: c.N_ESA_LEVELS]) - - def _get_esa_mode(histrates: xr.Dataset) -> int: """ Read the ESA mode of a pointing, defaulting to HiRes. @@ -222,6 +455,8 @@ class LoSpinAnglePointingSet(PointingSet): The values of the pointing, each of shape (esa level, spin angle). frame : SpiceFrame The frame to compute the sky directions in, i.e. the map's frame. + energy : np.ndarray + The energy [keV] of each ESA level. """ tiling_type: SkyTilingType = SkyTilingType.RECTANGULAR @@ -233,6 +468,7 @@ def __init__( spin_angles: np.ndarray, values: dict[str, np.ndarray], frame: SpiceFrame, + energy: np.ndarray, ): dims = [CoordNames.TIME.value, CoordNames.ENERGY_L2.value, "spin_angle"] super().__init__( @@ -243,7 +479,7 @@ def __init__( }, coords={ CoordNames.TIME.value: [epoch], - CoordNames.ENERGY_L2.value: _esa_energy(), + CoordNames.ENERGY_L2.value: energy, }, ), spice_reference_frame=frame, @@ -278,7 +514,7 @@ def midpoint_j2000_et(self) -> float: return float(ttj2000ns_to_et(self.epoch)) -def _initialize_accumulators(sky_map: RectangularSkyMap) -> None: +def _initialize_accumulators(sky_map: RectangularSkyMap, energy: np.ndarray) -> None: """ Seed the map with the empty accumulators each pointing is added into. @@ -291,6 +527,8 @@ def _initialize_accumulators(sky_map: RectangularSkyMap) -> None: ---------- sky_map : RectangularSkyMap The map being built, modified in place. + energy : np.ndarray + The energy [keV] of each ESA level. """ for name in ACCUMULATED_VARIABLES: sky_map.data_1d[name] = xr.DataArray( @@ -300,7 +538,7 @@ def _initialize_accumulators(sky_map: RectangularSkyMap) -> None: CoordNames.ENERGY_L2.value, CoordNames.GENERIC_PIXEL.value, ], - coords={CoordNames.ENERGY_L2.value: _esa_energy()}, + coords={CoordNames.ENERGY_L2.value: energy}, ) @@ -310,6 +548,7 @@ def _accumulate_pointing( histrates: xr.Dataset, sky_map: RectangularSkyMap, map_descriptor: MapDescriptor, + energy: np.ndarray, ) -> None: """ Add one pointing's counts and exposure to the map. @@ -328,6 +567,8 @@ def _accumulate_pointing( The map being built, modified in place. map_descriptor : MapDescriptor The parsed descriptor of the map being made. + energy : np.ndarray + The energy [keV] of each ESA level. """ species = map_descriptor.species pivot_angle = float(np.atleast_1d(goodtimes["pivot"].values)[0]) @@ -369,6 +610,7 @@ def _accumulate_pointing( "bg_rate_exposure": background_rates[:, np.newaxis] * pointing_exposure, }, sky_map.spice_reference_frame, + energy, ) # The projection sums the spin-angle bins that land in the same map pixel, # and adds this pointing on top of what the earlier pointings left there. @@ -444,7 +686,105 @@ def _spin_phase_mask( # ============================================================================= -# RATES AND INTENSITIES +# GEOMETRIC FACTORS +# ============================================================================= + + +def load_geometric_factor_data(species: str) -> pd.DataFrame: + """ + Load geometric factor data for the specified species. + + Parameters + ---------- + species : str + The species to load geometric factors for ("h" or "o"). + + Returns + ------- + pd.DataFrame + Geometric factor dataframe for the specified species. + + Raises + ------ + ValueError + If species is not "h" or "o". + """ + if species not in ["h", "o"]: + raise ValueError( + f"Geometric factors only available for 'h' and 'o', got '{species}'" + ) + + if species == "h": + gf_file = sorted(ANCILLARY_DATA_DIR.glob("*hydrogen-geometric-factor*"))[-1] + else: # species == "o" + gf_file = sorted(ANCILLARY_DATA_DIR.glob("*oxygen-geometric-factor*"))[-1] + + return lo_ancillary.read_ancillary_file(gf_file) + + +def reduce_geometric_factor_data(species: str, esa_mode: int) -> pd.DataFrame: + """ + Get geometric factor data for a specific species and ESA mode. + + This helper function loads geometric factor data, filters by ESA mode, and + selects the row of each of the 7 energy steps, in ascending step order. + + Parameters + ---------- + species : str + The species to load geometric factors for ("h" or "o"). + esa_mode : int + ESA mode (0 for HiRes, 1 for HiThr). + + Returns + ------- + pd.DataFrame + Geometric factor data indexed by Observed_E-Step (1-7), containing all + columns from the geometric factor CSV file. + """ + # Load geometric factor data for this species + gf_data = load_geometric_factor_data(species) + + # Filter for the specific ESA mode + if "esa_mode" in gf_data.columns: + gf_data = gf_data[gf_data["esa_mode"] == esa_mode] + + # Lo Instrument team: Use only geometric factors where + # incident_E-Step == Observed_E-Step + diagonal = gf_data["incident_E-Step"] == gf_data["Observed_E-Step"] + gf_data = gf_data[diagonal].set_index("Observed_E-Step") + + # Select the energy steps, in order. Raises if the file is missing one. + return gf_data.loc[list(range(1, c.N_ESA_LEVELS + 1))] + + +def _esa_energy(species: str, esa_mode: int) -> np.ndarray: + """ + Get the energy of each ESA level the map is binned in. + + The energies are the passband centers the geometric factors were measured + at, so they are read from the same ancillary file, for the same species and + ESA mode. + + Parameters + ---------- + species : str + The species of the map ("h" or "o"). + esa_mode : int + The ESA mode, 0 for HiRes and 1 for HiThr. + + Returns + ------- + np.ndarray + The energy [keV] of each ESA level. + """ + return reduce_geometric_factor_data(species, esa_mode)["Cntr_E"].to_numpy( + dtype=float + ) + + +# ============================================================================= +# RATES AND INTENSITIES CALCULATIONS # ============================================================================= @@ -475,7 +815,7 @@ def _geometric_factors(esa_mode: int) -> tuple[np.ndarray, np.ndarray, np.ndarra def _calculate_rates_and_intensities( - sky_map: RectangularSkyMap, esa_mode: int + sky_map: RectangularSkyMap, esa_mode: int, energy: np.ndarray ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -488,6 +828,8 @@ def _calculate_rates_and_intensities( The map the pointings were projected onto, read for its accumulators. esa_mode : int The ESA mode, 0 for HiRes and 1 for HiThr. + energy : np.ndarray + The energy [keV] of each ESA level. Returns ------- @@ -498,7 +840,7 @@ def _calculate_rates_and_intensities( exposure = sky_map.data_1d["exposure_factor"].values bg_rate_exposure = sky_map.data_1d["bg_rate_exposure"].values - energy = _esa_energy()[:, np.newaxis] + energy = energy[:, np.newaxis] geometric_factor, error_upper, error_lower = _geometric_factors(esa_mode) geometric_factor = geometric_factor[:, np.newaxis] error_upper = error_upper[:, np.newaxis] diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index d8ed50d7f..c7ebe643a 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -6,6 +6,7 @@ import pytest import xarray as xr +from imap_processing import imap_module_directory from imap_processing.cdf.utils import load_cdf, write_cdf from imap_processing.ena_maps.ena_maps import match_coords_to_indices from imap_processing.ena_maps.utils.naming import MapDescriptor @@ -19,9 +20,11 @@ ) from imap_processing.spice.time import met_to_ttj2000ns +ANCILLARY_DIR = imap_module_directory / "tests/lo/test_anc" + # A full-spin map, so that every spin-angle bin lands on it. -FULL_DESCRIPTOR = "l090-ena-h-sf-nsp-full-hae-6deg-3mo" -RAM_DESCRIPTOR = "l090-ena-h-sf-nsp-ram-hae-6deg-3mo" +FULL_DESCRIPTOR = "l090-enansnbs-h-sf-nsp-full-hae-6deg-3mo" +RAM_DESCRIPTOR = "l090-enansnbs-h-sf-nsp-ram-hae-6deg-3mo" N_ESA = LoConstants.N_ESA_LEVELS N_SPIN_BINS = LoConstants.N_SPIN_ANGLE_BINS @@ -33,6 +36,26 @@ IN_METS = [511_000_150.0, 511_000_200.0, 511_000_250.0] OUT_METS = [510_990_000.0, 511_010_000.0] +# The ESA level energies [keV] of imap_lo_hydrogen-geometric-factor-small, by +# ESA mode, which the map takes its energy binning from. +ESA_ENERGIES = { + 0: np.array([0.010, 0.020, 0.040, 0.080, 0.160, 0.320, 0.640]), + 1: np.array([0.011, 0.022, 0.044, 0.088, 0.176, 0.352, 0.704]), +} + + +@pytest.fixture(autouse=True) +def use_test_geometric_factors(): + """Point the map at the small geometric factor ancillary in ``test_anc``. + + The map reads its geometric factors, and the ESA level energies they were + measured at, straight out of the ancillary shipped with the package. The + test file stands in for it so the tests do not have to track the flight + calibration. + """ + with patch("imap_processing.lo.l2.lo_l2.ANCILLARY_DATA_DIR", ANCILLARY_DIR): + yield + def product_attrs(repointing, product): """The global attributes an L1B input of a pointing is written with.""" @@ -122,6 +145,7 @@ def identity_pointing(et, az_el, *args, **kwargs): def make_pointing_set(sky_map, spin_angles, pivot=PIVOT): """Build the in-memory pointing set of one pointing, sky pointing mocked.""" + energy = np.arange(1.0, N_ESA + 1.0) values = { name: np.ones((N_ESA, spin_angles.size)) for name in ("ena_count", "exposure_factor", "bg_rate_exposure") @@ -136,6 +160,7 @@ def make_pointing_set(sky_map, spin_angles, pivot=PIVOT): spin_angles, values, sky_map.spice_reference_frame, + energy, ) @@ -146,13 +171,25 @@ def one_pointing(): @pytest.fixture -def full_map(one_pointing): +def anc_dependencies(): + """The ancillary files a map takes as a dependency. + + Every map but a raw one is flux corrected, so the ESA eta fit factors are + required; without them the map cannot be made. + """ + return [ANCILLARY_DIR / "imap_lo_esa-eta-fit-factors_20240101_v001.csv"] + + +@pytest.fixture +def full_map(one_pointing, anc_dependencies): """The full-spin map of one pointing, with the sky pointing mocked.""" with patch( "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", side_effect=identity_pointing, ): - (dataset,) = lo_l2(as_dependencies(one_pointing), [], FULL_DESCRIPTOR) + (dataset,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR + ) return dataset, one_pointing @@ -201,12 +238,12 @@ def test_map_variables(self, full_map): ) def test_energy_coordinate(self, full_map): - """The energy coordinate and its widths come from the ESA constants.""" + """The energy coordinate comes from the ancillary, its widths from the + ESA constants.""" dataset, _ = full_map - np.testing.assert_allclose( - dataset["energy"].values, LoConstants.ESA_ENERGY[:N_ESA] - ) + # The pointings are all in ESA mode 0. + np.testing.assert_allclose(dataset["energy"].values, ESA_ENERGIES[0]) np.testing.assert_allclose( dataset["energy_delta_plus"].values, LoConstants.ESA_ENERGY_DELTA[0] ) @@ -249,7 +286,7 @@ def test_out_of_goodtime_epochs_are_excluded(self, full_map): assert per_energy.max() < 999.0 * N_SPIN_BINS np.testing.assert_allclose(per_energy, pointing["expected_counts"]) - def test_pointings_accumulate(self, one_pointing): + def test_pointings_accumulate(self, one_pointing, anc_dependencies): """Two pointings contribute twice the counts of one.""" other = make_pointing(repointing=101, seed=7) @@ -257,22 +294,30 @@ def test_pointings_accumulate(self, one_pointing): "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", side_effect=identity_pointing, ): - (one,) = lo_l2(as_dependencies(one_pointing), [], FULL_DESCRIPTOR) - (both,) = lo_l2(as_dependencies(one_pointing, other), [], FULL_DESCRIPTOR) + (one,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR + ) + (both,) = lo_l2( + as_dependencies(one_pointing, other), anc_dependencies, FULL_DESCRIPTOR + ) np.testing.assert_allclose( both["ena_count"].values.sum(axis=(0, 2, 3)), one["ena_count"].values.sum(axis=(0, 2, 3)) + other["expected_counts"], ) - def test_ram_map_keeps_half_the_spin(self, one_pointing): + def test_ram_map_keeps_half_the_spin(self, one_pointing, anc_dependencies): """A ram map takes fewer counts than the full spin it is cut from.""" with patch( "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", side_effect=identity_pointing, ): - (full,) = lo_l2(as_dependencies(one_pointing), [], FULL_DESCRIPTOR) - (ram,) = lo_l2(as_dependencies(one_pointing), [], RAM_DESCRIPTOR) + (full,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR + ) + (ram,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, RAM_DESCRIPTOR + ) full_counts = full["ena_count"].values.sum() ram_counts = ram["ena_count"].values.sum() @@ -301,7 +346,7 @@ def test_rate_and_intensity(self, full_map): geometric_factor = ( np.array(LoConstants.GEO_FACTOR[:N_ESA]) * LoConstants.GEO_FACTOR_SCALE ) - energy = np.array(LoConstants.ESA_ENERGY[:N_ESA]) + energy = ESA_ENERGIES[0] expected = dataset["ena_count_rate"] / xr.DataArray( geometric_factor * energy, dims=["energy"] ) @@ -454,20 +499,20 @@ def test_missing_product_raises(self, one_pointing): class TestUnsupported: """Map flavours the Lo pipeline does not make.""" - def test_oxygen_not_supported(self, one_pointing): + def test_oxygen_not_supported(self, one_pointing, anc_dependencies): """Only hydrogen geometric factors are defined.""" with pytest.raises(NotImplementedError, match="species o"): lo_l2( as_dependencies(one_pointing), - [], + anc_dependencies, "l090-ena-o-sf-nsp-full-hae-6deg-3mo", ) - def test_healpix_not_supported(self, one_pointing): + def test_healpix_not_supported(self, one_pointing, anc_dependencies): """Lo makes rectangular maps only.""" with pytest.raises(NotImplementedError, match="HEALPix"): lo_l2( as_dependencies(one_pointing), - [], - "l090-ena-h-sf-nsp-full-hae-nside8-3mo", + anc_dependencies, + "l090-enansnbs-h-sf-nsp-full-hae-nside8-3mo", ) From a7fe482fd647c2e053210202bcc4b5803119a319 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 30 Jul 2026 14:48:08 -0400 Subject: [PATCH 02/14] Removed GEO_FACTOR* constants - reading from ancillary file instead --- imap_processing/lo/constants.py | 30 ---------- imap_processing/lo/l2/lo_l2.py | 77 +++++++++++++++----------- imap_processing/tests/lo/test_lo_l2.py | 10 +++- 3 files changed, 51 insertions(+), 66 deletions(-) diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index 0f173bb6c..92fe4d521 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -79,36 +79,6 @@ class LoConstants: RAM_HISTOGRAM_BINS: tuple[slice, ...] = (slice(0, 20), slice(50, 60)) ANTI_RAM_HISTOGRAM_BINS: tuple[slice, ...] = (slice(20, 50),) - # The following are indexed by ESA level (0-indexed, ESA level = index + 1). - # The 8th entry is the virtual E8 channel, unused by the map. - GEO_FACTOR: ClassVar[list[float]] = [ - 7.0e-5, - 7.9e-5, - 9.7e-5, - 11.2e-5, - 14.0e-5, - 17.7e-5, - 22.5e-5, - 6.721e-5, - ] - GEO_FACTOR_ERR: ClassVar[list[float]] = [ - 4.9e-5, - 5.5e-5, - 6.8e-5, - 3.0e-5, - 4.5e-5, - 2.0e-5, - 1.4e-5, - 6.721e-5, - ] - - # GEO_FACTOR/GEO_FACTOR_ERR above are the raw, pre-recalibration values; the - # map multiplies them by GEO_FACTOR_SCALE, and derives the asymmetric - # upper/lower G-factor bounds using the two scale factors below. - GEO_FACTOR_SCALE: float = 0.63529412 - GEO_FACTOR_SCALE_UPPER: float = 1.57407407 - GEO_FACTOR_SCALE_LOWER: float = 0.36728395 - # Half-widths [keV] of the ESA energy passbands, by ESA level, for the two # ESA modes. NOTE: From an e-mail from Nathan on 2025-09-11 (converted to keV). ESA_ENERGY_DELTA: ClassVar[dict[int, list[float]]] = { diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index a60b8a2f9..d0223b722 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -106,11 +106,11 @@ def lo_l2( logger.info("Step 1: Loading ancillary data") _efficiency_data = load_efficiency_data(anc_dependencies) - # The geometric factors in LoConstants are hydrogen only. + # Only hydrogen maps are supported end to end for now. if map_descriptor.species != "h": raise NotImplementedError( f"Cannot make a map of species {map_descriptor.species} for " - f"{descriptor}. Only hydrogen geometric factors are defined." + f"{descriptor}. Only hydrogen maps are supported." ) sky_map = map_descriptor.to_empty_map() @@ -133,7 +133,9 @@ def lo_l2( goodtimes, bgrates, histrates, sky_map, map_descriptor, energy ) - variables = _calculate_rates_and_intensities(sky_map, esa_mode, energy) + variables = _calculate_rates_and_intensities( + sky_map, map_descriptor.species, esa_mode, energy + ) dataset = _build_map_dataset(sky_map, variables, esa_mode) logger.info("IMAP-Lo L2 processing pipeline completed successfully") @@ -788,34 +790,47 @@ def _esa_energy(species: str, esa_mode: int) -> np.ndarray: # ============================================================================= -def _geometric_factors(esa_mode: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]: +def _geometric_factors( + species: str, esa_mode: int +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ Get the recalibrated geometric factors and their asymmetric bounds. + The ancillary names its two uncertainty columns for the direction the + intensity derived from the factor moves in, which is the opposite of the + direction the factor itself moves in: intensity goes as 1/G, so a smaller + factor gives a larger intensity. Its ``_unc_plus`` is therefore the + downward excursion of the factor, and its ``_unc_minus`` the upward one. + Parameters ---------- + species : str + The species of the map ("h" or "o"). esa_mode : int - The ESA mode, 0 for HiRes and 1 for HiThr. Unused for now, the - geometric factors are not yet split by ESA mode. + The ESA mode, 0 for HiRes and 1 for HiThr. Returns ------- tuple[np.ndarray, np.ndarray, np.ndarray] - The geometric factor of each ESA level, and its upper and lower error - bounds. + The geometric factor of each ESA level, and its lower and upper + calibration bounds. """ - levels = slice(0, c.N_ESA_LEVELS) - geometric_factor = np.array(c.GEO_FACTOR[levels]) * c.GEO_FACTOR_SCALE - error = np.array(c.GEO_FACTOR_ERR[levels]) * c.GEO_FACTOR_SCALE + gf_data = reduce_geometric_factor_data(species, esa_mode) + column = f"GF_Trpl_{species.upper()}" - error_upper = np.hypot(geometric_factor * (c.GEO_FACTOR_SCALE_UPPER - 1.0), error) - error_lower = np.hypot(geometric_factor * (1.0 - c.GEO_FACTOR_SCALE_LOWER), error) + geometric_factor = gf_data[column].to_numpy(dtype=float) + excursion_down = gf_data[f"{column}_unc_plus"].to_numpy(dtype=float) + excursion_up = gf_data[f"{column}_unc_minus"].to_numpy(dtype=float) - return geometric_factor, error_upper, error_lower + return ( + geometric_factor, + geometric_factor - excursion_down, + geometric_factor + excursion_up, + ) def _calculate_rates_and_intensities( - sky_map: RectangularSkyMap, esa_mode: int, energy: np.ndarray + sky_map: RectangularSkyMap, species: str, esa_mode: int, energy: np.ndarray ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -826,6 +841,8 @@ def _calculate_rates_and_intensities( ---------- sky_map : RectangularSkyMap The map the pointings were projected onto, read for its accumulators. + species : str + The species of the map ("h" or "o"), which sets the geometric factors. esa_mode : int The ESA mode, 0 for HiRes and 1 for HiThr. energy : np.ndarray @@ -841,10 +858,10 @@ def _calculate_rates_and_intensities( bg_rate_exposure = sky_map.data_1d["bg_rate_exposure"].values energy = energy[:, np.newaxis] - geometric_factor, error_upper, error_lower = _geometric_factors(esa_mode) + geometric_factor, gf_low, gf_high = _geometric_factors(species, esa_mode) geometric_factor = geometric_factor[:, np.newaxis] - error_upper = error_upper[:, np.newaxis] - error_lower = error_lower[:, np.newaxis] + gf_low = gf_low[:, np.newaxis] + gf_high = gf_high[:, np.newaxis] exposed = exposure > 0 @@ -878,27 +895,21 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: intensity = _divide(count_rate, geometric_factor * energy) intensity_stat_uncert = _divide(count_rate_stat_uncert, geometric_factor * energy) - # The systematic error is the flux excursion from the recalibrated G-factor - # bounds: the upper/lower excursions come from the lower/upper G-factor - # bounds respectively, and the symmetric error is their geometric mean. It - # is undefined where the lower bound would drive the G-factor non-positive. - valid = geometric_factor > error_lower + # The systematic error is the intensity excursion from the recalibrated + # G-factor bounds, and the symmetric error is the geometric mean of the two. + # Intensity goes as 1/G, so the lower G-factor bound gives the upper + # intensity. It is undefined where that bound is not positive. + valid = gf_low > 0 if not valid.all(): logger.warning( "The geometric factor of ESA levels " f"{(np.flatnonzero(~valid[:, 0]) + 1).tolist()} is below its lower " f"error bound; their systematic errors are left at zero." ) - intensity_sys_err_plus = np.where( - valid, - intensity * geometric_factor / (geometric_factor - error_lower) - intensity, - 0.0, - ) - intensity_sys_err_minus = np.where( - valid, - intensity - intensity * geometric_factor / (geometric_factor + error_upper), - 0.0, - ) + intensity_upper = _divide(count_rate, np.where(valid, gf_low, 1.0) * energy) + intensity_lower = _divide(count_rate, gf_high * energy) + intensity_sys_err_plus = np.where(valid, intensity_upper - intensity, 0.0) + intensity_sys_err_minus = np.where(valid, intensity - intensity_lower, 0.0) bg_rate = _divide(bg_rate_exposure, exposure) bg_rate_stat_uncert = np.sqrt(_divide(bg_rate, exposure)) diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index c7ebe643a..5c8fb4ed6 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -43,6 +43,12 @@ 1: np.array([0.011, 0.022, 0.044, 0.088, 0.176, 0.352, 0.704]), } +# The hydrogen geometric factors of the same file, by ESA mode. +GEO_FACTORS = { + 0: np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]) * 1e-5, + 1: np.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7]) * 1e-5, +} + @pytest.fixture(autouse=True) def use_test_geometric_factors(): @@ -343,9 +349,7 @@ def test_rate_and_intensity(self, full_map): ) assert np.all(dataset["ena_count_rate"].values[~exposed] == 0) - geometric_factor = ( - np.array(LoConstants.GEO_FACTOR[:N_ESA]) * LoConstants.GEO_FACTOR_SCALE - ) + geometric_factor = GEO_FACTORS[0] energy = ESA_ENERGIES[0] expected = dataset["ena_count_rate"] / xr.DataArray( geometric_factor * energy, dims=["energy"] From 304a80862f84b577dfd391946d3e534f6f32b7c9 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 30 Jul 2026 15:57:24 -0400 Subject: [PATCH 03/14] reading energy band half-widths from the ancillary instead of in constants --- ...imap_lo_hydrogen-geometric-factor_v004.csv | 32 ++--- .../imap_lo_oxygen-geometric-factor_v004.csv | 32 ++--- imap_processing/lo/constants.py | 43 ++++-- imap_processing/lo/l2/lo_l2.py | 134 ++++++++---------- imap_processing/tests/lo/test_lo_l2.py | 11 +- 5 files changed, 131 insertions(+), 121 deletions(-) diff --git a/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv b/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv index f110e64c2..cb0ced20f 100644 --- a/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv +++ b/imap_processing/lo/ancillary_data/imap_lo_hydrogen-geometric-factor_v004.csv @@ -1,16 +1,16 @@ -esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_H,GF_Trpl_H_unc_minus,GF_Trpl_H_unc_plus -# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV] -0,1,1,0.01633,0.00028,4.45E-05,4.03E-05,4.20E-05 -0,2,2,0.03047,0.00043,5.02E-05,4.53E-05,4.72E-05 -0,3,3,0.05576,0.00089,6.16E-05,5.58E-05,5.82E-05 -0,4,4,0.10626,0.0019,7.12E-05,4.51E-05,4.89E-05 -0,5,5,0.20004,0.0034,8.89E-05,5.85E-05,6.31E-05 -0,6,6,0.40496,0.0073,1.12E-04,6.58E-05,7.23E-05 -0,7,7,0.78729,0.022,1.43E-04,8.25E-05,9.09E-05 -1,1,1,0.01719,0.00022,1.05E-04,8.82E-05,9.26E-05 -1,2,2,0.03236,0.00036,1.22E-04,1.02E-04,1.07E-04 -1,3,3,0.05948,0.00077,1.44E-04,1.21E-04,1.27E-04 -1,4,4,0.11441,0.0013,1.73E-04,1.17E-04,1.26E-04 -1,5,5,0.2137,0.003,2.15E-04,1.37E-04,1.49E-04 -1,6,6,0.43736,0.0053,2.82E-04,1.65E-04,1.81E-04 -1,7,7,0.83888,0.0117,3.61E-04,2.08E-04,2.29E-04 +esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_H,GF_Trpl_H_unc_minus,GF_Trpl_H_unc_plus,Cntr_E_delta_minus,Cntr_E_delta_plus +# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV] +0,1,1,0.01633,0.00028,4.45E-05,4.03E-05,4.20E-05,0.00543,0.00543 +0,2,2,0.03047,0.00043,5.02E-05,4.53E-05,4.72E-05,0.01002,0.01002 +0,3,3,0.05576,0.00089,6.16E-05,5.58E-05,5.82E-05,0.01861,0.01861 +0,4,4,0.10626,0.0019,7.12E-05,4.51E-05,4.89E-05,0.03331,0.03331 +0,5,5,0.20004,0.0034,8.89E-05,5.85E-05,6.31E-05,0.06498,0.06498 +0,6,6,0.40496,0.0073,1.12E-04,6.58E-05,7.23E-05,0.13164,0.13164 +0,7,7,0.78729,0.022,1.43E-04,8.25E-05,9.09E-05,0.26235,0.26235 +1,1,1,0.01719,0.00022,1.05E-04,8.82E-05,9.26E-05,0.00881,0.00881 +1,2,2,0.03236,0.00036,1.22E-04,1.02E-04,1.07E-04,0.01604,0.01604 +1,3,3,0.05948,0.00077,1.44E-04,1.21E-04,1.27E-04,0.02850,0.02850 +1,4,4,0.11441,0.0013,1.73E-04,1.17E-04,1.26E-04,0.05313,0.05313 +1,5,5,0.2137,0.003,2.15E-04,1.37E-04,1.49E-04,0.10560,0.10560 +1,6,6,0.43736,0.0053,2.82E-04,1.65E-04,1.81E-04,0.21967,0.21967 +1,7,7,0.83888,0.0117,3.61E-04,2.08E-04,2.29E-04,0.41360,0.41360 diff --git a/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv b/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv index 0b07d8199..9f2c8cff3 100644 --- a/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv +++ b/imap_processing/lo/ancillary_data/imap_lo_oxygen-geometric-factor_v004.csv @@ -1,16 +1,16 @@ -esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_O,GF_Trpl_O_unc_minus,GF_Trpl_O_unc_plus -# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV] -0,1,1,0.01919,0.001372344105,1.98E-05,1.74E-05,1.81E-05 -0,2,2,0.03675,0.002537963082,3.18E-05,2.39E-05,2.53E-05 -0,3,3,0.07121,0.006591339559,5.34E-05,4.05E-05,4.29E-05 -0,4,4,0.14141,0.01743990907,8.64E-05,6.14E-05,6.56E-05 -0,5,5,0.274,0.03583900763,1.32E-04,9.07E-05,9.72E-05 -0,6,6,0.58503,0.09806681395,2.46E-04,1.53E-04,1.67E-04 -0,7,7,1.13506,0.2381076888,3.81E-04,2.23E-04,2.45E-04 -1,1,1,0.02043303552,0.001575924599,5.02E-05,4.41E-05,4.61E-05 -1,2,2,0.03972,0.002953020557,8.13E-05,6.09E-05,6.47E-05 -1,3,3,0.07648,0.007476494389,1.33E-04,1.01E-04,1.07E-04 -1,4,4,0.15353,0.01487953071,2.38E-04,1.69E-04,1.80E-04 -1,5,5,0.29846,0.04148950349,3.87E-04,2.67E-04,2.86E-04 -1,6,6,0.61524,0.08357629594,6.61E-04,4.11E-04,4.47E-04 -1,7,7,1.24282,0.2120806857,1.02E-03,6.10E-04,6.67E-04 +esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_O,GF_Trpl_O_unc_minus,GF_Trpl_O_unc_plus,Cntr_E_delta_minus,Cntr_E_delta_plus +# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV] +0,1,1,0.01919,0.001372344105,1.98E-05,1.74E-05,1.81E-05,0.00582,0.00582 +0,2,2,0.03675,0.002537963082,3.18E-05,2.39E-05,2.53E-05,0.01110,0.01110 +0,3,3,0.07121,0.006591339559,5.34E-05,4.05E-05,4.29E-05,0.02178,0.02178 +0,4,4,0.14141,0.01743990907,8.64E-05,6.14E-05,6.56E-05,0.04147,0.04147 +0,5,5,0.274,0.03583900763,1.32E-04,9.07E-05,9.72E-05,0.08561,0.08561 +0,6,6,0.58503,0.09806681395,2.46E-04,1.53E-04,1.67E-04,0.18067,0.18067 +0,7,7,1.13506,0.2381076888,3.81E-04,2.23E-04,2.45E-04,0.36193,0.36193 +1,1,1,0.02043303552,0.001575924599,5.02E-05,4.41E-05,4.61E-05,0.00945,0.00945 +1,2,2,0.03972,0.002953020557,8.13E-05,6.09E-05,6.47E-05,0.01784,0.01784 +1,3,3,0.07648,0.007476494389,1.33E-04,1.01E-04,1.07E-04,0.03351,0.03351 +1,4,4,0.15353,0.01487953071,2.38E-04,1.69E-04,1.80E-04,0.06661,0.06661 +1,5,5,0.29846,0.04148950349,3.87E-04,2.67E-04,2.86E-04,0.13995,0.13995 +1,6,6,0.61524,0.08357629594,6.61E-04,4.11E-04,4.47E-04,0.30224,0.30224 +1,7,7,1.24282,0.2120806857,1.02E-03,6.10E-04,6.67E-04,0.56948,0.56948 diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index 92fe4d521..87119c436 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -3,6 +3,40 @@ from dataclasses import dataclass from typing import ClassVar, NamedTuple +import numpy as np + + +class EsaCalibration(NamedTuple): + """ + The instrument's calibration settings for each ESA level. + + Every field holds one value per ESA level, in ascending level order, read + from the ancillary of one species and ESA mode. + + Attributes + ---------- + energy : np.ndarray + The energy [keV] of each level, the passband center its geometric + factor was measured at. + energy_delta_minus : np.ndarray + The half-width [keV] of each passband below its center. + energy_delta_plus : np.ndarray + The half-width [keV] of each passband above its center. + geometric_factor : np.ndarray + The recalibrated geometric factor [cm^2 sr keV/keV] of each level. + geometric_factor_low : np.ndarray + The lower calibration bound of each geometric factor. + geometric_factor_high : np.ndarray + The upper calibration bound of each geometric factor. + """ + + energy: np.ndarray + energy_delta_minus: np.ndarray + energy_delta_plus: np.ndarray + geometric_factor: np.ndarray + geometric_factor_low: np.ndarray + geometric_factor_high: np.ndarray + class PivotAngleSpec(NamedTuple): """ @@ -79,15 +113,6 @@ class LoConstants: RAM_HISTOGRAM_BINS: tuple[slice, ...] = (slice(0, 20), slice(50, 60)) ANTI_RAM_HISTOGRAM_BINS: tuple[slice, ...] = (slice(20, 50),) - # Half-widths [keV] of the ESA energy passbands, by ESA level, for the two - # ESA modes. NOTE: From an e-mail from Nathan on 2025-09-11 (converted to keV). - ESA_ENERGY_DELTA: ClassVar[dict[int, list[float]]] = { - # esa_mode 0, HiRes - 0: [0.00543, 0.01002, 0.01861, 0.03331, 0.06498, 0.13164, 0.26235], - # esa_mode 1, HiThr - 1: [0.00881, 0.01604, 0.02850, 0.05313, 0.10560, 0.21967, 0.41360], - } - # Nominal background rates [counts/s] for each species BG_RATES: ClassVar[dict[str, float]] = {"H": 0.0014925, "O": 0.000136635} # When no exposure is available, scale the nominal rate down as a conservative diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index d0223b722..e5882b15d 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -16,6 +16,7 @@ from imap_processing.ena_maps.utils.coordinates import CoordNames from imap_processing.ena_maps.utils.naming import MapDescriptor from imap_processing.lo import lo_ancillary +from imap_processing.lo.constants import EsaCalibration from imap_processing.lo.constants import LoConstants as c # noqa: N813 from imap_processing.lo.l1c.lo_l1c import compute_pointing_directions from imap_processing.spice.geometry import ( @@ -40,6 +41,7 @@ # The calibration ancillaries shipped with the package. ANCILLARY_DATA_DIR = Path(__file__).parent.parent / "ancillary_data" + # ============================================================================= # MAIN ENTRY POINT # ============================================================================= @@ -69,8 +71,9 @@ def lo_l2( The input datasets covering the pointings of the map window, keyed by repointing and then by product descriptor. anc_dependencies : list - List of ancillary file paths. Unused, the calibration constants of the - map live in ``LoConstants``. + List of ancillary file paths, read for the efficiency and correction + factors. The geometric factors and ESA level energies come from the + ancillaries shipped with the package, in ``ANCILLARY_DATA_DIR``. descriptor : str The map descriptor to be produced (e.g., "l090-ena-h-sf-nsp-ram-hae-6deg-3mo"). @@ -121,22 +124,20 @@ def lo_l2( logger.info(f"Building {descriptor} from {len(pointings)} pointings") # Every pointing of a map is taken in the same ESA mode, so the last one - # sets the energies and passband widths the whole map is binned in. + # sets the energy response the whole map is binned in. esa_mode = _get_esa_mode(pointings[max(pointings)][2]) if pointings else 0 - energy = _esa_energy(map_descriptor.species, esa_mode) + calibration = _esa_calibration(map_descriptor.species, esa_mode) - _initialize_accumulators(sky_map, energy) + _initialize_accumulators(sky_map, calibration.energy) for repointing, (goodtimes, bgrates, histrates) in sorted(pointings.items()): logger.debug(f"Accumulating repoint{repointing:05d}") _accumulate_pointing( - goodtimes, bgrates, histrates, sky_map, map_descriptor, energy + goodtimes, bgrates, histrates, sky_map, map_descriptor, calibration.energy ) - variables = _calculate_rates_and_intensities( - sky_map, map_descriptor.species, esa_mode, energy - ) - dataset = _build_map_dataset(sky_map, variables, esa_mode) + variables = _calculate_rates_and_intensities(sky_map, calibration) + dataset = _build_map_dataset(sky_map, variables, calibration) logger.info("IMAP-Lo L2 processing pipeline completed successfully") return [ @@ -760,13 +761,16 @@ def reduce_geometric_factor_data(species: str, esa_mode: int) -> pd.DataFrame: return gf_data.loc[list(range(1, c.N_ESA_LEVELS + 1))] -def _esa_energy(species: str, esa_mode: int) -> np.ndarray: +def _esa_calibration(species: str, esa_mode: int) -> EsaCalibration: """ - Get the energy of each ESA level the map is binned in. + Get the ESA level calibration one map is built from. - The energies are the passband centers the geometric factors were measured - at, so they are read from the same ancillary file, for the same species and - ESA mode. + The ancillary names its two geometric factor uncertainty columns for the + direction the intensity derived from the factor moves in, which is the + opposite of the direction the factor itself moves in: intensity goes as + 1/G, so a smaller factor gives a larger intensity. Its ``_unc_plus`` is + therefore the downward excursion of the factor, and its ``_unc_minus`` the + upward one. Parameters ---------- @@ -777,11 +781,24 @@ def _esa_energy(species: str, esa_mode: int) -> np.ndarray: Returns ------- - np.ndarray - The energy [keV] of each ESA level. - """ - return reduce_geometric_factor_data(species, esa_mode)["Cntr_E"].to_numpy( - dtype=float + EsaCalibration + The energies, passband half-widths and geometric factors of every ESA + level, in ascending level order. + """ + gf_data = reduce_geometric_factor_data(species, esa_mode).astype(float) + + factor = f"GF_Trpl_{species.upper()}" + geometric_factor = gf_data[factor].to_numpy() + + return EsaCalibration( + energy=gf_data["Cntr_E"].to_numpy(), + energy_delta_minus=gf_data["Cntr_E_delta_minus"].to_numpy(), + energy_delta_plus=gf_data["Cntr_E_delta_plus"].to_numpy(), + geometric_factor=geometric_factor, + geometric_factor_low=geometric_factor + - gf_data[f"{factor}_unc_plus"].to_numpy(), + geometric_factor_high=geometric_factor + + gf_data[f"{factor}_unc_minus"].to_numpy(), ) @@ -790,47 +807,8 @@ def _esa_energy(species: str, esa_mode: int) -> np.ndarray: # ============================================================================= -def _geometric_factors( - species: str, esa_mode: int -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - Get the recalibrated geometric factors and their asymmetric bounds. - - The ancillary names its two uncertainty columns for the direction the - intensity derived from the factor moves in, which is the opposite of the - direction the factor itself moves in: intensity goes as 1/G, so a smaller - factor gives a larger intensity. Its ``_unc_plus`` is therefore the - downward excursion of the factor, and its ``_unc_minus`` the upward one. - - Parameters - ---------- - species : str - The species of the map ("h" or "o"). - esa_mode : int - The ESA mode, 0 for HiRes and 1 for HiThr. - - Returns - ------- - tuple[np.ndarray, np.ndarray, np.ndarray] - The geometric factor of each ESA level, and its lower and upper - calibration bounds. - """ - gf_data = reduce_geometric_factor_data(species, esa_mode) - column = f"GF_Trpl_{species.upper()}" - - geometric_factor = gf_data[column].to_numpy(dtype=float) - excursion_down = gf_data[f"{column}_unc_plus"].to_numpy(dtype=float) - excursion_up = gf_data[f"{column}_unc_minus"].to_numpy(dtype=float) - - return ( - geometric_factor, - geometric_factor - excursion_down, - geometric_factor + excursion_up, - ) - - def _calculate_rates_and_intensities( - sky_map: RectangularSkyMap, species: str, esa_mode: int, energy: np.ndarray + sky_map: RectangularSkyMap, calibration: EsaCalibration ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -841,12 +819,9 @@ def _calculate_rates_and_intensities( ---------- sky_map : RectangularSkyMap The map the pointings were projected onto, read for its accumulators. - species : str - The species of the map ("h" or "o"), which sets the geometric factors. - esa_mode : int - The ESA mode, 0 for HiRes and 1 for HiThr. - energy : np.ndarray - The energy [keV] of each ESA level. + calibration : EsaCalibration + The energy response the map is binned in, read for the energies and + geometric factors the intensities are derived with. Returns ------- @@ -857,11 +832,11 @@ def _calculate_rates_and_intensities( exposure = sky_map.data_1d["exposure_factor"].values bg_rate_exposure = sky_map.data_1d["bg_rate_exposure"].values - energy = energy[:, np.newaxis] - geometric_factor, gf_low, gf_high = _geometric_factors(species, esa_mode) - geometric_factor = geometric_factor[:, np.newaxis] - gf_low = gf_low[:, np.newaxis] - gf_high = gf_high[:, np.newaxis] + # Every ESA level quantity gets a pixel axis to broadcast over the map. + energy = calibration.energy[:, np.newaxis] + geometric_factor = calibration.geometric_factor[:, np.newaxis] + gf_low = calibration.geometric_factor_low[:, np.newaxis] + gf_high = calibration.geometric_factor_high[:, np.newaxis] exposed = exposure > 0 @@ -936,7 +911,9 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: def _build_map_dataset( - sky_map: RectangularSkyMap, variables: dict[str, np.ndarray], esa_mode: int + sky_map: RectangularSkyMap, + variables: dict[str, np.ndarray], + calibration: EsaCalibration, ) -> xr.Dataset: """ Lay the map variables out on the map's sky grid. @@ -950,8 +927,8 @@ def _build_map_dataset( The map being built. variables : dict[str, np.ndarray] The map variables, each of shape (epoch, esa level, pixel). - esa_mode : int - The ESA mode, 0 for HiRes and 1 for HiThr, which sets the widths of the + calibration : EsaCalibration + The energy response the map is binned in, read for the widths of the ESA energy passbands. Returns @@ -968,8 +945,11 @@ def _build_map_dataset( dataset = sky_map.to_dataset() - energy_delta = np.array(c.ESA_ENERGY_DELTA[esa_mode]) - dataset["energy_delta_minus"] = xr.DataArray(energy_delta, dims=["energy"]) - dataset["energy_delta_plus"] = xr.DataArray(energy_delta, dims=["energy"]) + dataset["energy_delta_minus"] = xr.DataArray( + calibration.energy_delta_minus, dims=["energy"] + ) + dataset["energy_delta_plus"] = xr.DataArray( + calibration.energy_delta_plus, dims=["energy"] + ) return dataset diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 5c8fb4ed6..c8ac878fb 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -49,6 +49,9 @@ 1: np.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7]) * 1e-5, } +# The ESA passband half-widths of the same file, a tenth of each center energy. +ESA_ENERGY_DELTAS = {mode: energies / 10 for mode, energies in ESA_ENERGIES.items()} + @pytest.fixture(autouse=True) def use_test_geometric_factors(): @@ -244,14 +247,16 @@ def test_map_variables(self, full_map): ) def test_energy_coordinate(self, full_map): - """The energy coordinate comes from the ancillary, its widths from the - ESA constants.""" + """The energy coordinate and its widths both come from the ancillary.""" dataset, _ = full_map # The pointings are all in ESA mode 0. np.testing.assert_allclose(dataset["energy"].values, ESA_ENERGIES[0]) np.testing.assert_allclose( - dataset["energy_delta_plus"].values, LoConstants.ESA_ENERGY_DELTA[0] + dataset["energy_delta_minus"].values, ESA_ENERGY_DELTAS[0] + ) + np.testing.assert_allclose( + dataset["energy_delta_plus"].values, ESA_ENERGY_DELTAS[0] ) def test_map_writes_to_cdf(self, full_map): From 42a1b2736bc53b619ce70125dfe7c99f123d2f20 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 30 Jul 2026 16:17:37 -0400 Subject: [PATCH 04/14] added missing test ancillaries --- ...p_lo_esa-eta-fit-factors_20240101_v001.csv | 8 ++++++++ ...o_hydrogen-geometric-factor-small_v001.csv | 19 +++++++++++++++++++ 2 files changed, 27 insertions(+) create mode 100644 imap_processing/tests/lo/test_anc/imap_lo_esa-eta-fit-factors_20240101_v001.csv create mode 100644 imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv diff --git a/imap_processing/tests/lo/test_anc/imap_lo_esa-eta-fit-factors_20240101_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_esa-eta-fit-factors_20240101_v001.csv new file mode 100644 index 000000000..4fda625ba --- /dev/null +++ b/imap_processing/tests/lo/test_anc/imap_lo_esa-eta-fit-factors_20240101_v001.csv @@ -0,0 +1,8 @@ +esa_step,M0,M1,M2,M3,M4,M5 +1,1.0152,-0.047723,0.0304,-0.001817,0.0023649,-0.00032519 +2,1.0142,-0.045778,0.030061,-0.0020424,0.0021796,-0.00029036 +3,1.013,-0.045334,0.030649,-0.0021426,0.0021755,-0.0002869 +4,1.0109,-0.043671,0.029741,-0.0014197,0.0016756,-0.0002298 +5,1.0145,-0.045219,0.029705,-0.0021726,0.0020739,-0.0002652 +6,1.0116,-0.043433,0.030755,-0.0020747,0.001899,-0.0002476 +7,1.0156,-0.048728,0.029868,-0.0016762,0.0022885,-0.0003183 \ No newline at end of file diff --git a/imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv new file mode 100644 index 000000000..0a6733a9e --- /dev/null +++ b/imap_processing/tests/lo/test_anc/imap_lo_hydrogen-geometric-factor-small_v001.csv @@ -0,0 +1,19 @@ +esa_mode,incident_E-Step,Observed_E-Step,Cntr_E,Cntr_E_unc,GF_Trpl_H,GF_Trpl_H_unc_minus,GF_Trpl_H_unc_plus,Cntr_E_delta_minus,Cntr_E_delta_plus +# [1],[1],[1],[keV],[keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[cm^2sr keV/keV],[keV],[keV] +0,1,1,0.01000,0.00010,1.00E-05,1.00E-06,1.00E-06,0.00100,0.00100 +0,2,1,0.01000,0.00010,9.99E-05,9.99E-06,9.99E-06,0.00100,0.00100 +0,2,2,0.02000,0.00020,2.00E-05,2.00E-06,2.00E-06,0.00200,0.00200 +0,3,2,0.02000,0.00020,9.99E-05,9.99E-06,9.99E-06,0.00200,0.00200 +0,3,3,0.04000,0.00040,3.00E-05,3.00E-06,3.00E-06,0.00400,0.00400 +0,4,4,0.08000,0.00080,4.00E-05,4.00E-06,4.00E-06,0.00800,0.00800 +0,5,5,0.16000,0.00160,5.00E-05,5.00E-06,5.00E-06,0.01600,0.01600 +0,6,6,0.32000,0.00320,6.00E-05,6.00E-06,6.00E-06,0.03200,0.03200 +0,7,7,0.64000,0.00640,7.00E-05,7.00E-06,7.00E-06,0.06400,0.06400 +1,1,1,0.01100,0.00011,1.10E-05,1.10E-06,1.10E-06,0.00110,0.00110 +1,2,1,0.01100,0.00011,9.99E-05,9.99E-06,9.99E-06,0.00110,0.00110 +1,2,2,0.02200,0.00022,2.20E-05,2.20E-06,2.20E-06,0.00220,0.00220 +1,3,3,0.04400,0.00044,3.30E-05,3.30E-06,3.30E-06,0.00440,0.00440 +1,4,4,0.08800,0.00088,4.40E-05,4.40E-06,4.40E-06,0.00880,0.00880 +1,5,5,0.17600,0.00176,5.50E-05,5.50E-06,5.50E-06,0.01760,0.01760 +1,6,6,0.35200,0.00352,6.60E-05,6.60E-06,6.60E-06,0.03520,0.03520 +1,7,7,0.70400,0.00704,7.70E-05,7.70E-06,7.70E-06,0.07040,0.07040 From bda2ff51a1087f3e445138a2197a61a76c67f694 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 3 Aug 2026 10:32:29 -0400 Subject: [PATCH 05/14] sputter correction logic --- imap_processing/ena_maps/utils/naming.py | 44 +++- ...map_lo_sputter-correction-factors_v002.csv | 9 + imap_processing/lo/l2/lo_l2.py | 247 ++++++++++++------ imap_processing/tests/ena_maps/test_naming.py | 42 +++ ..._sputter-correction-factors-small_v001.csv | 4 + imap_processing/tests/lo/test_lo_l2.py | 140 +++++++++- 6 files changed, 398 insertions(+), 88 deletions(-) create mode 100644 imap_processing/lo/ancillary_data/imap_lo_sputter-correction-factors_v002.csv create mode 100644 imap_processing/tests/lo/test_anc/imap_lo_sputter-correction-factors-small_v001.csv diff --git a/imap_processing/ena_maps/utils/naming.py b/imap_processing/ena_maps/utils/naming.py index 191239d2c..6bada65c8 100644 --- a/imap_processing/ena_maps/utils/naming.py +++ b/imap_processing/ena_maps/utils/naming.py @@ -42,6 +42,12 @@ class MappableInstrumentShortName(Enum): "hk": "heliocentric kinetic", } +# The principal data part of a descriptor is a stem naming the quantity mapped, +# optionally followed by modifier codes saying which corrections were made. +PRINCIPAL_DATA_PATTERN = re.compile( + r"^(drt|ena|int|isn|spx)(?:(?<=spx)\d+)?([^-_\s]*)$" +) + @dataclass class MapDescriptor: @@ -251,9 +257,7 @@ def build_catdesc(self, quantity_text: str) -> str: instrument = instrument.title() sensor = " Combined" if self.sensor == "combined" else self.sensor species = "UV" if self.species == "uv" else self.species.title() - m = re.match( - r"^(drt|ena|int|isn|spx)(?:(?<=spx)\d+)?([^-_\s]*)$", self.principal_data - ) + m = PRINCIPAL_DATA_PATTERN.match(self.principal_data) if m.group(1) == "isn": species = "ISN " + species extras = m.group(2) @@ -310,6 +314,40 @@ def principal_data_var(self) -> str: "spx": "ena_spectral_index", }[self.principal_data[:3]] + @property + def principal_data_extras(self) -> str: + """ + The modifier codes following the principal data stem. + + These say which corrections a map was made with, e.g. the "nbs" of + "enanbs". Empty if the principal data is a bare stem, or if it does not + parse as a stem followed by modifiers at all. + + Returns + ------- + principal_data_extras : str + The modifier codes, run together in the order they appear. + """ + m = PRINCIPAL_DATA_PATTERN.match(self.principal_data) + return m.group(2) if m else "" + + @property + def sputter_corrected(self) -> bool: + """ + Whether the map has the counts sputtered in from a heavier species removed. + + Returns + ------- + sputter_corrected : bool + True if this map is sputter corrected. + """ + if not self.principal_data.startswith("ena"): + return False + extras = self.principal_data_extras + # "ns" is the code for a map made with no sputter correction, "nbs" the + # older code for one made with neither sputter nor bootstrap correction. + return not extras.startswith(("ns", "nbs")) + # Methods for parsing and building parts of the map descriptor string @staticmethod def get_instrument_descriptor( diff --git a/imap_processing/lo/ancillary_data/imap_lo_sputter-correction-factors_v002.csv b/imap_processing/lo/ancillary_data/imap_lo_sputter-correction-factors_v002.csv new file mode 100644 index 000000000..7455a2e3a --- /dev/null +++ b/imap_processing/lo/ancillary_data/imap_lo_sputter-correction-factors_v002.csv @@ -0,0 +1,9 @@ +source_species,target_species,target_esa,source_esa,sputter_factor +o,h,1,4,0.236 +o,h,2,4,0.372 +o,h,3,4,0.898 +o,h,4,4,0.891 +o,h,5,4,0.037 +o,h,5,6,0.32 +o,h,6,6,0.32 +o,h,7,6,0.22 diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index e5882b15d..b967ae059 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -94,17 +94,13 @@ def lo_l2( map_descriptor = MapDescriptor.from_string(descriptor) logger.info(f"Processing map for species: {map_descriptor.species}") - # Determine if corrections are needed and prepare oxygen data if required + # Determine which of the corrections the descriptor asks for are needed ( - _sputtering_correction, - _bootstrap_correction, _flux_correction, - _o_map_dataset, _flux_factors, + sputter_correction, _cg_correction, - ) = _prepare_corrections( - map_descriptor, descriptor, sci_dependencies, anc_dependencies - ) + ) = _prepare_corrections(map_descriptor, anc_dependencies) logger.info("Step 1: Loading ancillary data") _efficiency_data = load_efficiency_data(anc_dependencies) @@ -128,15 +124,33 @@ def lo_l2( esa_mode = _get_esa_mode(pointings[max(pointings)][2]) if pointings else 0 calibration = _esa_calibration(map_descriptor.species, esa_mode) - _initialize_accumulators(sky_map, calibration.energy) + # The species sputtering into this map, if it is to be sputter corrected, + # and the ESA levels it sputters into. Its counts are accumulated on the + # same grid, alongside the map's own. + sputter_source, sputter_matrix = ( + _sputter_correction(map_descriptor.species) + if sputter_correction + else (None, None) + ) + accumulators = ACCUMULATED_VARIABLES + ( + ("sputter_source_count",) if sputter_source else () + ) + + _initialize_accumulators(sky_map, calibration.energy, accumulators) for repointing, (goodtimes, bgrates, histrates) in sorted(pointings.items()): logger.debug(f"Accumulating repoint{repointing:05d}") _accumulate_pointing( - goodtimes, bgrates, histrates, sky_map, map_descriptor, calibration.energy + goodtimes, + bgrates, + histrates, + sky_map, + map_descriptor, + calibration.energy, + sputter_source, ) - variables = _calculate_rates_and_intensities(sky_map, calibration) + variables = _calculate_rates_and_intensities(sky_map, calibration, sputter_matrix) dataset = _build_map_dataset(sky_map, variables, calibration) logger.info("IMAP-Lo L2 processing pipeline completed successfully") @@ -152,60 +166,33 @@ def lo_l2( def _prepare_corrections( map_descriptor: MapDescriptor, - descriptor: str, - sci_dependencies: dict, anc_dependencies: list, -) -> tuple[bool, bool, bool, xr.Dataset | None, Path | None, bool]: +) -> tuple[bool, Path | None, bool, bool]: """ - Determine what corrections are needed and prepare oxygen dataset if required. - - This helper function encapsulates the logic for determining when sputtering - and bootstrap corrections should be applied, and handles the creation of - the oxygen dataset needed for sputtering corrections. + Determine which of the corrections the map descriptor asks for are needed. Parameters ---------- map_descriptor : MapDescriptor The parsed map descriptor containing species and data type information. - descriptor : str - The original descriptor string for creating the oxygen variant. - sci_dependencies : dict - Dictionary of datasets needed for L2 data product creation. anc_dependencies : list List of ancillary file paths. Returns ------- - tuple[bool, bool, bool, xr.Dataset | None, Path | None, bool] + tuple[bool, Path | None, bool, bool] A tuple containing: - - sputtering_correction: Whether to apply sputtering corrections - - bootstrap_correction: Whether to apply bootstrap corrections - flux_correction: Whether to apply flux corrections - - o_map_dataset: Oxygen dataset if needed, None otherwise - flux_factors: Path to flux factors ancillary file if needed, None otherwise + - sputter_correction: Whether to remove the counts sputtered into the + mapped species from a heavier one. - cg_correction: Whether to apply CG correction to the dataset. """ # Default values - no corrections needed - sputtering_correction = False - bootstrap_correction = False flux_correction = False - o_map_dataset = None flux_factors: None | Path = None - # Sputtering and bootstrap corrections are only applied to hydrogen ENA data - # Guard against recursion: don't process oxygen for oxygen maps - if ( - map_descriptor.species == "h" - and map_descriptor.principal_data == "ena" - and "-o-" not in descriptor - ): # Safety check to prevent infinite recursion - logger.info("Creating map for oxygen for sputtering corrections") - o_descriptor = descriptor.replace("-h-", "-o-") - o_map_dataset = lo_l2(sci_dependencies, anc_dependencies, o_descriptor)[0] - sputtering_correction = True - bootstrap_correction = True - if "raw" not in map_descriptor.principal_data: flux_correction = True try: @@ -217,14 +204,13 @@ def _prepare_corrections( "No flux correction factor file found in ancillary dependencies" ) from None + sputter_correction = map_descriptor.sputter_corrected cg_correction = True if map_descriptor.frame_descriptor == "hf" else False return ( - sputtering_correction, - bootstrap_correction, flux_correction, - o_map_dataset, flux_factors, + sputter_correction, cg_correction, ) @@ -266,40 +252,72 @@ def load_efficiency_data(anc_dependencies: list) -> pd.DataFrame: ) -def load_sputter_correction_data( - source_species: str, target_species: str -) -> pd.DataFrame: +def load_sputter_correction_data() -> pd.DataFrame: """ - Load sputter correction factors from an ancillary file. - - Parameters - ---------- - source_species : str - The species doing the sputtering (e.g. "o" for oxygen). - target_species : str - The species being corrected (e.g. "h" for hydrogen). + Load the sputter correction factors shipped with the package. Returns ------- pd.DataFrame - Rows matching the given species pair, sorted ascending by esa_step, - with columns: source_species, target_species, esa_step, - sputter_factor, sputter_factor_uncertainty. + The ancillary data, with columns: source_species, target_species, + target_esa, source_esa, sputter_factor. + + Raises + ------ + ValueError + If no sputter correction ancillary is shipped with the package. """ sputter_files = sorted(ANCILLARY_DATA_DIR.glob("*sputter-correction-factors*")) if not sputter_files: raise ValueError("No sputter correction files found") - df = pd.concat( - [lo_ancillary.read_ancillary_file(f) for f in sputter_files], - ignore_index=True, - ) - mask = (df["source_species"] == source_species) & ( - df["target_species"] == target_species - ) - result = df[mask].sort_values("esa_step").reset_index(drop=True) - return result + return lo_ancillary.read_ancillary_file(sputter_files[-1]) + + +def _sputter_correction(species: str) -> tuple[str, np.ndarray]: + """ + Load the sputter correction of a species from the ancillary. + + Parameters + ---------- + species : str + The species being mapped, whose counts are to be corrected. + + Returns + ------- + tuple[str, np.ndarray] + The species sputtering into the mapped species, and the (target level, + source level) fraction of its counts to remove from the mapped + species' counts, zero where a source level does not sputter into a + target level. + + Raises + ------ + ValueError + If the ancillary has nothing sputtering into the mapped species, or + names more than one species doing it, which the correction cannot + choose between. + """ + factors = load_sputter_correction_data() + factors = factors[factors["target_species"] == species] + sources = factors["source_species"].unique() + + if len(sources) == 0: + raise ValueError( + f"The map asks for a sputter correction, but the ancillary has no " + f"factors for {species}" + ) + if len(sources) > 1: + raise ValueError( + f"More than one species sputters into {species}: {sorted(sources)}" + ) + + matrix = np.zeros((c.N_ESA_LEVELS, c.N_ESA_LEVELS)) + matrix[ + factors["target_esa"].to_numpy() - 1, factors["source_esa"].to_numpy() - 1 + ] = factors["sputter_factor"].to_numpy() + return str(sources[0]), matrix def load_bootstrap_correction_data() -> pd.DataFrame: @@ -517,7 +535,9 @@ def midpoint_j2000_et(self) -> float: return float(ttj2000ns_to_et(self.epoch)) -def _initialize_accumulators(sky_map: RectangularSkyMap, energy: np.ndarray) -> None: +def _initialize_accumulators( + sky_map: RectangularSkyMap, energy: np.ndarray, names: tuple[str, ...] +) -> None: """ Seed the map with the empty accumulators each pointing is added into. @@ -532,8 +552,10 @@ def _initialize_accumulators(sky_map: RectangularSkyMap, energy: np.ndarray) -> The map being built, modified in place. energy : np.ndarray The energy [keV] of each ESA level. + names : tuple[str, ...] + The accumulators to seed. """ - for name in ACCUMULATED_VARIABLES: + for name in names: sky_map.data_1d[name] = xr.DataArray( np.zeros((1, c.N_ESA_LEVELS, sky_map.num_points)), dims=[ @@ -552,6 +574,7 @@ def _accumulate_pointing( sky_map: RectangularSkyMap, map_descriptor: MapDescriptor, energy: np.ndarray, + sputter_source: str | None = None, ) -> None: """ Add one pointing's counts and exposure to the map. @@ -572,6 +595,9 @@ def _accumulate_pointing( The parsed descriptor of the map being made. energy : np.ndarray The energy [keV] of each ESA level. + sputter_source : str | None + The species sputtering into the mapped species, whose counts are + accumulated alongside. None if this map is not sputter corrected. """ species = map_descriptor.species pivot_angle = float(np.atleast_1d(goodtimes["pivot"].values)[0]) @@ -600,18 +626,24 @@ def _accumulate_pointing( # The whole pointing is projected from the middle of its good times, which # is where the despun frame is sampled. epoch = met_to_ttj2000ns((gt_start.min() + gt_end.max()) / 2.0) + values = { + "ena_count": pointing_counts, + "exposure_factor": pointing_exposure, + # Background is a rate per ESA level per pointing, so it is + # accumulated weighted by exposure and divided by the total + # exposure at the end. + "bg_rate_exposure": background_rates[:, np.newaxis] * pointing_exposure, + } + if sputter_source: + values["sputter_source_count"] = ( + histrates[f"{sputter_source}_counts"].values[in_goodtime].sum(axis=0) + ) + pointing_set = LoSpinAnglePointingSet( epoch, pivot_angle, spin_angles, - { - "ena_count": pointing_counts, - "exposure_factor": pointing_exposure, - # Background is a rate per ESA level per pointing, so it is - # accumulated weighted by exposure and divided by the total - # exposure at the end. - "bg_rate_exposure": background_rates[:, np.newaxis] * pointing_exposure, - }, + values, sky_map.spice_reference_frame, energy, ) @@ -619,7 +651,7 @@ def _accumulate_pointing( # and adds this pointing on top of what the earlier pointings left there. sky_map.project_pset_values_to_map( pointing_set, - value_keys=list(ACCUMULATED_VARIABLES), + value_keys=list(values), pset_valid_mask=keep, ) @@ -807,8 +839,40 @@ def _esa_calibration(species: str, esa_mode: int) -> EsaCalibration: # ============================================================================= +def _sputter_correct_counts( + counts: np.ndarray, source_counts: np.ndarray, sputter_matrix: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """ + Remove the counts sputtered into the mapped species from another species. + + Parameters + ---------- + counts : np.ndarray + The accumulated counts of the mapped species, of shape + (epoch, esa level, pixel). + source_counts : np.ndarray + The accumulated counts of the sputtering species, same shape and grid. + sputter_matrix : np.ndarray + The (target level, source level) fraction of the source counts to + remove from the target counts. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + The corrected counts and their variance. Counting in both species is + Poisson, so each source term contributes its counts scaled by the + square of its factor, and the variance only ever grows. + """ + logger.info("Applying the sputter correction to the accumulated counts") + corrected = counts - np.einsum("ts,esp->etp", sputter_matrix, source_counts) + variance = counts + np.einsum("ts,esp->etp", sputter_matrix**2, source_counts) + return corrected, variance + + def _calculate_rates_and_intensities( - sky_map: RectangularSkyMap, calibration: EsaCalibration + sky_map: RectangularSkyMap, + calibration: EsaCalibration, + sputter_matrix: np.ndarray | None = None, ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -822,6 +886,10 @@ def _calculate_rates_and_intensities( calibration : EsaCalibration The energy response the map is binned in, read for the energies and geometric factors the intensities are derived with. + sputter_matrix : np.ndarray | None + The (target level, source level) sputter correction factors, applied + to the counts before the rate is taken. None leaves the counts as + they were observed. Returns ------- @@ -832,6 +900,13 @@ def _calculate_rates_and_intensities( exposure = sky_map.data_1d["exposure_factor"].values bg_rate_exposure = sky_map.data_1d["bg_rate_exposure"].values + if sputter_matrix is None: + rate_counts, rate_counts_var = counts, counts + else: + rate_counts, rate_counts_var = _sputter_correct_counts( + counts, sky_map.data_1d["sputter_source_count"].values, sputter_matrix + ) + # Every ESA level quantity gets a pixel axis to broadcast over the map. energy = calibration.energy[:, np.newaxis] geometric_factor = calibration.geometric_factor[:, np.newaxis] @@ -863,9 +938,11 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: where=exposed, ) - count_rate = _divide(counts, exposure) + # Removing the sputtered counts can take a low-count pixel below zero, + # which is not a rate the instrument can have observed. + count_rate = np.maximum(_divide(rate_counts, exposure), 0.0) # Poisson uncertainty on the counts, propagated to the rate - count_rate_stat_uncert = _divide(np.sqrt(counts), exposure) + count_rate_stat_uncert = _divide(np.sqrt(rate_counts_var), exposure) intensity = _divide(count_rate, geometric_factor * energy) intensity_stat_uncert = _divide(count_rate_stat_uncert, geometric_factor * energy) @@ -940,8 +1017,10 @@ def _build_map_dataset( dims = sky_map.data_1d["ena_count"].dims for name, values in variables.items(): sky_map.data_1d[name] = xr.DataArray(values.astype(np.float32), dims=dims) - # `bg_rate_exposure` is an accumulator, not a map variable. - sky_map.data_1d = sky_map.data_1d.drop_vars("bg_rate_exposure") + # These are accumulators, not map variables. + sky_map.data_1d = sky_map.data_1d.drop_vars( + ["bg_rate_exposure", "sputter_source_count"], errors="ignore" + ) dataset = sky_map.to_dataset() diff --git a/imap_processing/tests/ena_maps/test_naming.py b/imap_processing/tests/ena_maps/test_naming.py index cc32debf5..b1707a2b4 100644 --- a/imap_processing/tests/ena_maps/test_naming.py +++ b/imap_processing/tests/ena_maps/test_naming.py @@ -505,3 +505,45 @@ def test_build_map_var_catdesc_more_vars( def test_principal_data_var(self, descriptor_str, expected_principal_data_var): md = MapDescriptor.from_string(descriptor_str) assert md.principal_data_var == expected_principal_data_var + + @pytest.mark.parametrize( + "principal_data, expected_extras", + [ + ("ena", ""), + ("enas", "s"), + ("enans", "ns"), + ("enanbs", "nbs"), + ("enansnbs", "nsnbs"), + ("isnnbkgnd", "nbkgnd"), + # The digits of a spx stem belong to the stem, not the modifiers + ("spx0305", ""), + ], + ) + def test_principal_data_extras(self, principal_data, expected_extras): + md = MapDescriptor.from_string( + f"l090-{principal_data}-h-sf-nsp-ram-hae-6deg-1yr" + ) + assert md.principal_data_extras == expected_extras + + @pytest.mark.parametrize( + "principal_data, expected", + [ + # "s" and "ns" after the stem say whether the correction was made + ("enas", True), + ("enans", False), + # The codes mark the exceptions, so an unmarked ENA map is corrected + ("ena", True), + # "nbs" (no sputter/bootstrap) suppresses it as well + ("enanbs", False), + ("enansnbs", False), + # Only ENA maps are sputter corrected + ("spx0305", False), + ("isn", False), + ("drt", False), + ], + ) + def test_sputter_corrected(self, principal_data, expected): + md = MapDescriptor.from_string( + f"l090-{principal_data}-h-sf-nsp-ram-hae-6deg-1yr" + ) + assert md.sputter_corrected is expected diff --git a/imap_processing/tests/lo/test_anc/imap_lo_sputter-correction-factors-small_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_sputter-correction-factors-small_v001.csv new file mode 100644 index 000000000..15a2db61f --- /dev/null +++ b/imap_processing/tests/lo/test_anc/imap_lo_sputter-correction-factors-small_v001.csv @@ -0,0 +1,4 @@ +source_species,target_species,target_esa,source_esa,sputter_factor +o,h,2,3,0.5 +o,h,5,3,0.25 +o,h,5,6,0.1 diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index c8ac878fb..06b6faea8 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -22,10 +22,18 @@ ANCILLARY_DIR = imap_module_directory / "tests/lo/test_anc" -# A full-spin map, so that every spin-angle bin lands on it. +# A full-spin map, so that every spin-angle bin lands on it. The "ns" after +# "ena" asks for no sputter correction, so these are the uncorrected maps. FULL_DESCRIPTOR = "l090-enansnbs-h-sf-nsp-full-hae-6deg-3mo" RAM_DESCRIPTOR = "l090-enansnbs-h-sf-nsp-ram-hae-6deg-3mo" +# The same full-spin map, sputter corrected. +SPUTTER_DESCRIPTOR = "l090-enas-h-sf-nsp-full-hae-6deg-3mo" + +# The contents of imap_lo_sputter-correction-factors-small, as +# {target ESA step: {source ESA step: factor}}, 1-based as in the ancillary. +SPUTTER_FACTORS = {2: {3: 0.5}, 5: {3: 0.25, 6: 0.1}} + N_ESA = LoConstants.N_ESA_LEVELS N_SPIN_BINS = LoConstants.N_SPIN_ANGLE_BINS PIVOT = 90.0 @@ -98,6 +106,10 @@ def make_pointing(repointing=100, pivot=PIVOT, seed=42): histrates = xr.Dataset( { "h_counts": (["epoch", "esa_step", "spin_bin_6"], counts), + # The sputter correction reads the oxygen counts of the same + # pointing. Making them the hydrogen counts lets a test predict the + # correction from the hydrogen counts the uncorrected map reports. + "o_counts": (["epoch", "esa_step", "spin_bin_6"], counts), "exposure_time_6deg": (["epoch", "esa_step", "spin_bin_6"], exposure), "esa_mode": ("epoch", np.zeros(mets.size, dtype=int)), }, @@ -202,6 +214,24 @@ def full_map(one_pointing, anc_dependencies): return dataset, one_pointing +@pytest.fixture +def sputter_maps(one_pointing, anc_dependencies): + """The sputter corrected and uncorrected maps of the same pointing. + + The two differ only in whether the correction was applied, so the + uncorrected map supplies the counts the correction is predicted from. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (corrected,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, SPUTTER_DESCRIPTOR + ) + (raw,) = lo_l2(as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR) + return corrected, raw + + class TestMapStructure: """The shape and contents of the produced map.""" @@ -412,6 +442,114 @@ def test_systematic_error_bounds(self, full_map): assert np.all(plus[lit] >= minus[lit]) +class TestSputterCorrection: + """Removing the counts oxygen sputters into the hydrogen channels.""" + + @staticmethod + def sputtered(counts, target_esa, power=1): + """The counts sputtered into a target ESA step, from the source steps. + + ``power`` is 1 for the counts themselves and 2 for their variance, + which each source term contributes to scaled by the square of its + factor. + """ + return sum( + factor**power * counts[:, source - 1] + for source, factor in SPUTTER_FACTORS.get(target_esa, {}).items() + ) + + def test_correction_removes_the_sputtered_counts(self, sputter_maps): + """The rate is the counts less the sputtered ones, over the exposure.""" + corrected, raw = sputter_maps + + counts = raw["ena_count"].values + exposure = raw["exposure_factor"].values + + for target_esa in range(1, N_ESA + 1): + exposed = exposure[:, target_esa - 1] > 0 + expected = np.maximum( + counts[:, target_esa - 1] - self.sputtered(counts, target_esa), 0.0 + ) + np.testing.assert_allclose( + corrected["ena_count_rate"].values[:, target_esa - 1][exposed], + expected[exposed] / exposure[:, target_esa - 1][exposed], + rtol=1e-5, + ) + + def test_uncorrected_steps_are_untouched(self, sputter_maps): + """A step that nothing sputters into keeps the rate it already had.""" + corrected, raw = sputter_maps + + untouched = [esa for esa in range(1, N_ESA + 1) if esa not in SPUTTER_FACTORS] + assert untouched, "the test factors must leave some steps uncorrected" + + for target_esa in untouched: + np.testing.assert_allclose( + corrected["ena_count_rate"].values[:, target_esa - 1], + raw["ena_count_rate"].values[:, target_esa - 1], + rtol=1e-6, + ) + + def test_corrected_steps_lose_intensity(self, sputter_maps): + """The corrected steps come out below the uncorrected ones somewhere.""" + corrected, raw = sputter_maps + + for target_esa in SPUTTER_FACTORS: + correction = ( + raw["ena_intensity"].values[:, target_esa - 1] + - corrected["ena_intensity"].values[:, target_esa - 1] + ) + assert np.all(correction >= 0) + assert np.any(correction > 0), f"ESA {target_esa} was not corrected" + + def test_uncertainty_gains_the_source_counts(self, sputter_maps): + """Subtracting a measured quantity can only add to the variance.""" + corrected, raw = sputter_maps + + counts = raw["ena_count"].values + exposure = raw["exposure_factor"].values + + for target_esa in range(1, N_ESA + 1): + exposed = exposure[:, target_esa - 1] > 0 + variance = counts[:, target_esa - 1] + self.sputtered( + counts, target_esa, power=2 + ) + np.testing.assert_allclose( + corrected["ena_count_rate_stat_uncert"].values[:, target_esa - 1][ + exposed + ], + np.sqrt(variance[exposed]) / exposure[:, target_esa - 1][exposed], + rtol=1e-5, + ) + + def test_rate_is_never_negative(self, sputter_maps): + """Over-subtracting a low-count pixel floors it rather than going below 0.""" + corrected, raw = sputter_maps + + counts = raw["ena_count"].values + # The test pointing is sparse enough that some pixel is over-subtracted, + # which is the case this floor exists for. + over_subtracted = [ + counts[:, esa - 1] - self.sputtered(counts, esa) < 0 + for esa in SPUTTER_FACTORS + ] + assert np.any(over_subtracted), "no pixel exercised the floor" + + assert np.all(corrected["ena_count_rate"].values >= 0) + assert np.all(corrected["ena_intensity"].values >= 0) + + def test_counts_stay_as_observed(self, sputter_maps): + """The correction applies from the rate onward, not to the raw counts.""" + corrected, raw = sputter_maps + + np.testing.assert_array_equal( + corrected["ena_count"].values, raw["ena_count"].values + ) + np.testing.assert_array_equal( + corrected["exposure_factor"].values, raw["exposure_factor"].values + ) + + class TestGeometry: """The spin-angle to sky-pixel geometry.""" From 958a8a2bbbb714238e92a68370b189a508bf64da Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 3 Aug 2026 13:50:40 -0400 Subject: [PATCH 06/14] bootstrap correction logic --- imap_processing/ena_maps/utils/naming.py | 36 ++- imap_processing/lo/constants.py | 19 ++ imap_processing/lo/l2/lo_l2.py | 282 +++++++++++++++++- imap_processing/tests/ena_maps/test_naming.py | 37 ++- ...ootstrap-correction-factors-small_v001.csv | 5 + imap_processing/tests/lo/test_lo_l2.py | 243 ++++++++++++++- 6 files changed, 598 insertions(+), 24 deletions(-) create mode 100644 imap_processing/tests/lo/test_anc/imap_lo_bootstrap-correction-factors-small_v001.csv diff --git a/imap_processing/ena_maps/utils/naming.py b/imap_processing/ena_maps/utils/naming.py index 6bada65c8..f05aec6fe 100644 --- a/imap_processing/ena_maps/utils/naming.py +++ b/imap_processing/ena_maps/utils/naming.py @@ -331,6 +331,18 @@ def principal_data_extras(self) -> str: m = PRINCIPAL_DATA_PATTERN.match(self.principal_data) return m.group(2) if m else "" + @property + def raw(self) -> bool: + """ + Whether the map is of the raw data, made with none of the corrections. + + Returns + ------- + raw : bool + True if this map asks for no corrections at all. + """ + return "raw" in self.principal_data + @property def sputter_corrected(self) -> bool: """ @@ -341,12 +353,26 @@ def sputter_corrected(self) -> bool: sputter_corrected : bool True if this map is sputter corrected. """ - if not self.principal_data.startswith("ena"): + if not self.principal_data.startswith("ena") or self.raw: + return False + # The stem is followed by the correction codes, in order: "s" or "ns" + # for the sputter correction, then "bs" or "nbs" for the bootstrap one. + return self.principal_data_extras.startswith("s") + + @property + def bootstrap_corrected(self) -> bool: + """ + Whether the map has the intensity bled in from higher ESA levels removed. + + Returns + ------- + bootstrap_corrected : bool + True if this map is bootstrap corrected. + """ + if not self.principal_data.startswith("ena") or self.raw: return False - extras = self.principal_data_extras - # "ns" is the code for a map made with no sputter correction, "nbs" the - # older code for one made with neither sputter nor bootstrap correction. - return not extras.startswith(("ns", "nbs")) + # The bootstrap code follows the sputter one, see sputter_corrected. + return re.match(r"(?:n?s)?bs", self.principal_data_extras) is not None # Methods for parsing and building parts of the map descriptor string @staticmethod diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index 87119c436..612da05d9 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -96,6 +96,25 @@ class LoConstants: N_SPINS_PER_ESA_LEVEL: int = 4 # Spins per ESA step within one histogram cycle N_SPIN_ANGLE_BINS: int = 60 # Number of angular bins within a spin + # Bootstrap correction settings. The nominal coefficients of the ancillary + # are scaled by BOOTSTRAP_SCALE before they are applied; the low and high + # scalings bracket that choice and become the systematic error on the + # corrected intensity. + BOOTSTRAP_SCALE: float = 0.5 + BOOTSTRAP_SCALE_INTENSITY_HIGH: float = 0.25 + BOOTSTRAP_SCALE_INTENSITY_LOW: float = 1.0 + + # The bootstrap correction of the highest ESA levels needs an ESA level + # above them to subtract. That virtual "ESA 8" channel has no geometric + # factor of its own; its intensity is extrapolated from the top two levels + # with a power law, at this multiple of the top level's center energy. + ESA_8_ENERGY_RATIO: float = 2.1 + # Width [pixels] of the neighborhood the spectral index of a pixel that has + # no measurable one is taken from. + BOOTSTRAP_SPECTRAL_INDEX_FILTER_SIZE: int = 3 + # The spectral index to extrapolate with when the map has none to offer. + BOOTSTRAP_DEFAULT_SPECTRAL_INDEX: float = 1.6 + # Nominal spin period [s]. True spin duration is NOT 15 seconds. NOMINAL_SPIN_PERIOD_SEC: float = 15.0 diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index b967ae059..f4a9b2e72 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -6,6 +6,7 @@ import numpy as np import pandas as pd import xarray as xr +from scipy.ndimage import generic_filter from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.ena_maps.ena_maps import ( @@ -99,6 +100,7 @@ def lo_l2( _flux_correction, _flux_factors, sputter_correction, + bootstrap_correction, _cg_correction, ) = _prepare_corrections(map_descriptor, anc_dependencies) @@ -150,7 +152,11 @@ def lo_l2( sputter_source, ) - variables = _calculate_rates_and_intensities(sky_map, calibration, sputter_matrix) + bootstrap_matrix = _bootstrap_correction() if bootstrap_correction else None + + variables = _calculate_rates_and_intensities( + sky_map, calibration, sputter_matrix, bootstrap_matrix + ) dataset = _build_map_dataset(sky_map, variables, calibration) logger.info("IMAP-Lo L2 processing pipeline completed successfully") @@ -167,7 +173,7 @@ def lo_l2( def _prepare_corrections( map_descriptor: MapDescriptor, anc_dependencies: list, -) -> tuple[bool, Path | None, bool, bool]: +) -> tuple[bool, Path | None, bool, bool, bool]: """ Determine which of the corrections the map descriptor asks for are needed. @@ -180,20 +186,22 @@ def _prepare_corrections( Returns ------- - tuple[bool, Path | None, bool, bool] + tuple[bool, Path | None, bool, bool, bool] A tuple containing: - flux_correction: Whether to apply flux corrections - flux_factors: Path to flux factors ancillary file if needed, None otherwise - sputter_correction: Whether to remove the counts sputtered into the mapped species from a heavier one. + - bootstrap_correction: Whether to remove the intensity that bled into + each ESA level from the levels above it. - cg_correction: Whether to apply CG correction to the dataset. """ # Default values - no corrections needed flux_correction = False flux_factors: None | Path = None - if "raw" not in map_descriptor.principal_data: + if not map_descriptor.raw: flux_correction = True try: flux_factors = next( @@ -205,12 +213,14 @@ def _prepare_corrections( ) from None sputter_correction = map_descriptor.sputter_corrected + bootstrap_correction = map_descriptor.bootstrap_corrected cg_correction = True if map_descriptor.frame_descriptor == "hf" else False return ( flux_correction, flux_factors, sputter_correction, + bootstrap_correction, cg_correction, ) @@ -330,16 +340,41 @@ def load_bootstrap_correction_data() -> pd.DataFrame: Bootstrap correction factors with columns: esa_step_i, esa_step_k, bootstrap_factor. Indices are 1-based ESA step numbers where esa_step_k=8 refers to the virtual E8 channel. + + Raises + ------ + ValueError + If no bootstrap correction ancillary is shipped with the package. """ bootstrap_files = sorted(ANCILLARY_DATA_DIR.glob("*bootstrap-correction-factors*")) if not bootstrap_files: raise ValueError("No bootstrap correction factor files found") - return pd.concat( - [lo_ancillary.read_ancillary_file(f) for f in bootstrap_files], - ignore_index=True, - ) + return lo_ancillary.read_ancillary_file(bootstrap_files[-1]) + + +def _bootstrap_correction() -> np.ndarray: + """ + Load the bootstrap coefficients from the ancillary. + + Returns + ------- + np.ndarray + The (target level, source level) fraction of a source level's intensity + to remove from a target level below it, zero where a source level does + not bleed into a target level. The source axis carries one level more + than the target axis, for the virtual ESA level above the top of the + map. These are the nominal coefficients; the correction scales them + itself, once for the correction and once for each of its bounds. + """ + factors = load_bootstrap_correction_data() + + matrix = np.zeros((c.N_ESA_LEVELS, c.N_ESA_LEVELS + 1)) + matrix[ + factors["esa_step_i"].to_numpy() - 1, factors["esa_step_k"].to_numpy() - 1 + ] = factors["bootstrap_factor"].to_numpy() + return matrix def finalize_dataset(dataset: xr.Dataset, descriptor: str) -> xr.Dataset: @@ -869,10 +904,222 @@ def _sputter_correct_counts( return corrected, variance +def _local_median_spectral_index( + spectral_index: np.ndarray, grid_shape: tuple[int, ...] +) -> np.ndarray: + """ + Fill in the spectral index of a pixel from the pixels around it. + + Parameters + ---------- + spectral_index : np.ndarray + The measured spectral index of every pixel, of shape (epoch, pixel), + NaN where the pixel has none. + grid_shape : tuple[int, ...] + The (azimuth, elevation) shape the pixel axis unwraps to, which is what + makes two pixels neighbors. + + Returns + ------- + np.ndarray + The median of the measured indices in each pixel's neighborhood, of the + same shape, NaN where the whole neighborhood is unmeasured. The map + does not wrap around in azimuth, so the pixels at its edges take the + median of the neighbors they have. + """ + + def median_of_measured(neighborhood: np.ndarray) -> float: + """ + Take the median of the pixels of a neighborhood that have an index. + + Parameters + ---------- + neighborhood : np.ndarray + The spectral indices of the pixels around (and including) one + pixel, NaN where a pixel has none. + + Returns + ------- + float + The median, or NaN if no pixel of the neighborhood has an index. + """ + measured = neighborhood[~np.isnan(neighborhood)] + return float(np.median(measured)) if measured.size else float(np.nan) + + return np.stack( + [ + generic_filter( + epoch_index.reshape(grid_shape), + median_of_measured, + size=c.BOOTSTRAP_SPECTRAL_INDEX_FILTER_SIZE, + mode="constant", + cval=np.nan, + ).ravel() + for epoch_index in spectral_index + ] + ) + + +def _extrapolate_top_intensity( + intensity: np.ndarray, energy: np.ndarray, grid_shape: tuple[int, ...] +) -> np.ndarray: + """ + Extrapolate the intensity of the virtual ESA level above the top of the map. + + The top ESA levels have nothing above them in the map to be bootstrap + corrected against, so a virtual level is extrapolated from the top two + levels of each pixel, taking the spectrum between them as a power law. + + Parameters + ---------- + intensity : np.ndarray + The intensity of every ESA level, of shape (epoch, esa level, pixel). + energy : np.ndarray + The energy [keV] of each ESA level. + grid_shape : tuple[int, ...] + The (azimuth, elevation) shape the pixel axis unwraps to. + + Returns + ------- + np.ndarray + The intensity of the virtual level, of shape (epoch, pixel), zero in + the pixels the top level saw nothing in. + """ + second, top = energy[-2], energy[-1] + virtual = top * c.ESA_8_ENERGY_RATIO + second_intensity, top_intensity = intensity[:, -2], intensity[:, -1] + + # The spectral index the two levels of a pixel imply, where it has both. + measured = (second_intensity > 0) & (top_intensity > 0) + spectral_index = np.zeros_like(top_intensity) + spectral_index[measured] = -np.log( + top_intensity[measured] / second_intensity[measured] + ) / np.log(top / second) + + extrapolated = np.zeros_like(top_intensity) + extrapolated[measured] = top_intensity[measured] * (virtual / top) ** ( + -spectral_index[measured] + ) + + # A pixel the second level saw nothing in has no spectrum of its own to + # extrapolate along, so it borrows one from its neighbors, falling back to + # the whole map and then to a nominal index. + borrowing = (top_intensity > 0) & ~measured + if borrowing.any(): + local = _local_median_spectral_index( + np.where(measured, spectral_index, np.nan), grid_shape + ) + global_index = ( + float(np.median(spectral_index[measured])) + if measured.any() + else c.BOOTSTRAP_DEFAULT_SPECTRAL_INDEX + ) + borrowed = np.where(np.isfinite(local), local, global_index) + extrapolated[borrowing] = top_intensity[borrowing] * (virtual / top) ** ( + -borrowed[borrowing] + ) + + return extrapolated + + +def _bootstrap_correct_intensity( + intensity: np.ndarray, + variance: np.ndarray, + calibration: EsaCalibration, + bootstrap_matrix: np.ndarray, + grid_shape: tuple[int, ...], + valid_gf_bounds: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Remove the intensity that bled into each ESA level from the levels above it. + + Every level is corrected against the intensities as they were measured, so + the correction of one level does not feed the correction of the next. + + Parameters + ---------- + intensity : np.ndarray + The intensity of every ESA level, of shape (epoch, esa level, pixel). + variance : np.ndarray + The statistical variance of those intensities, same shape. + calibration : EsaCalibration + The energy response the map is binned in, read for the energies the + virtual level is extrapolated along and the geometric factor bounds the + systematic error is taken from. + bootstrap_matrix : np.ndarray + The (target level, source level) nominal bootstrap coefficients. + grid_shape : tuple[int, ...] + The (azimuth, elevation) shape the pixel axis unwraps to. + valid_gf_bounds : np.ndarray + Whether the lower geometric factor bound of each ESA level is usable, + of shape (esa level, 1). The systematic error is zero where it is not. + + Returns + ------- + tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray] + The corrected intensity, its statistical uncertainty, and the upward + and downward systematic excursions of the correction. + """ + logger.info("Applying the bootstrap correction to the intensities") + + # The correction subtracts the level above the top of the map as well, so + # the intensities are extended by the virtual level it is extrapolated to. + # Its variance is approximated by that of the level it was extrapolated + # from, which contributes little to the levels it is subtracted from. + top_intensity = _extrapolate_top_intensity( + intensity, calibration.energy, grid_shape + ) + extended = np.concatenate([intensity, top_intensity[:, np.newaxis]], axis=1) + extended_variance = np.concatenate([variance, variance[:, -1:]], axis=1) + + def subtract(scale: float) -> np.ndarray: + """ + Subtract the bled intensity, at one scaling of the coefficients. + + Parameters + ---------- + scale : float + The scaling of the nominal coefficients. + + Returns + ------- + np.ndarray + The corrected intensity, which the subtraction can take below zero + in a faint pixel, floored there. + """ + return np.maximum( + intensity - np.einsum("ik,ekp->eip", scale * bootstrap_matrix, extended), + 0.0, + ) + + corrected = subtract(c.BOOTSTRAP_SCALE) + corrected_variance = variance + np.einsum( + "ik,ekp->eip", (c.BOOTSTRAP_SCALE * bootstrap_matrix) ** 2, extended_variance + ) + + # The systematic error spans the two scalings the correction is bracketed + # by, each moved on by the geometric factor bound of the same direction. + geometric_factor = calibration.geometric_factor[:, np.newaxis] + gf_high = calibration.geometric_factor_high[:, np.newaxis] + gf_low = np.where( + valid_gf_bounds, calibration.geometric_factor_low[:, np.newaxis], 1.0 + ) + lower = subtract(c.BOOTSTRAP_SCALE_INTENSITY_LOW) * geometric_factor / gf_high + upper = subtract(c.BOOTSTRAP_SCALE_INTENSITY_HIGH) * geometric_factor / gf_low + + return ( + corrected, + np.sqrt(corrected_variance), + np.where(valid_gf_bounds, upper - corrected, 0.0), + np.where(valid_gf_bounds, corrected - lower, 0.0), + ) + + def _calculate_rates_and_intensities( sky_map: RectangularSkyMap, calibration: EsaCalibration, sputter_matrix: np.ndarray | None = None, + bootstrap_matrix: np.ndarray | None = None, ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -890,6 +1137,10 @@ def _calculate_rates_and_intensities( The (target level, source level) sputter correction factors, applied to the counts before the rate is taken. None leaves the counts as they were observed. + bootstrap_matrix : np.ndarray | None + The (target level, source level) bootstrap correction coefficients, + applied to the intensities once they are derived. None leaves the + intensities as the counts gave them. Returns ------- @@ -963,6 +1214,21 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: intensity_sys_err_plus = np.where(valid, intensity_upper - intensity, 0.0) intensity_sys_err_minus = np.where(valid, intensity - intensity_lower, 0.0) + if bootstrap_matrix is not None: + ( + intensity, + intensity_stat_uncert, + intensity_sys_err_plus, + intensity_sys_err_minus, + ) = _bootstrap_correct_intensity( + intensity, + intensity_stat_uncert**2, + calibration, + bootstrap_matrix, + sky_map.binning_grid_shape, + valid, + ) + bg_rate = _divide(bg_rate_exposure, exposure) bg_rate_stat_uncert = np.sqrt(_divide(bg_rate, exposure)) bg_intensity = _divide(bg_rate, geometric_factor * energy) diff --git a/imap_processing/tests/ena_maps/test_naming.py b/imap_processing/tests/ena_maps/test_naming.py index b1707a2b4..195c23876 100644 --- a/imap_processing/tests/ena_maps/test_naming.py +++ b/imap_processing/tests/ena_maps/test_naming.py @@ -528,14 +528,14 @@ def test_principal_data_extras(self, principal_data, expected_extras): @pytest.mark.parametrize( "principal_data, expected", [ - # "s" and "ns" after the stem say whether the correction was made - ("enas", True), - ("enans", False), - # The codes mark the exceptions, so an unmarked ENA map is corrected - ("ena", True), - # "nbs" (no sputter/bootstrap) suppresses it as well - ("enanbs", False), + # "s" and "ns" after the stem say whether the correction was made, + # ahead of the "bs" or "nbs" of the bootstrap correction + ("enasbs", True), + ("enasnbs", True), + ("enansbs", False), ("enansnbs", False), + # A raw map asks for none of the corrections + ("enaraw", False), # Only ENA maps are sputter corrected ("spx0305", False), ("isn", False), @@ -547,3 +547,26 @@ def test_sputter_corrected(self, principal_data, expected): f"l090-{principal_data}-h-sf-nsp-ram-hae-6deg-1yr" ) assert md.sputter_corrected is expected + + @pytest.mark.parametrize( + "principal_data, expected", + [ + # "bs" and "nbs" say whether the correction was made, following the + # "s" or "ns" of the sputter correction + ("enasbs", True), + ("enasnbs", False), + ("enansbs", True), + ("enansnbs", False), + # A raw map asks for none of the corrections + ("enaraw", False), + # Only ENA maps are bootstrap corrected + ("spx0305", False), + ("isn", False), + ("drt", False), + ], + ) + def test_bootstrap_corrected(self, principal_data, expected): + md = MapDescriptor.from_string( + f"l090-{principal_data}-h-sf-nsp-ram-hae-6deg-1yr" + ) + assert md.bootstrap_corrected is expected diff --git a/imap_processing/tests/lo/test_anc/imap_lo_bootstrap-correction-factors-small_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_bootstrap-correction-factors-small_v001.csv new file mode 100644 index 000000000..c7148d13d --- /dev/null +++ b/imap_processing/tests/lo/test_anc/imap_lo_bootstrap-correction-factors-small_v001.csv @@ -0,0 +1,5 @@ +esa_step_i,esa_step_k,bootstrap_factor +2,3,0.4 +2,5,0.2 +6,7,0.5 +7,8,0.6 diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 06b6faea8..92fabbe5e 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -13,8 +13,11 @@ from imap_processing.lo.constants import LoConstants from imap_processing.lo.l2.lo_l2 import ( LoSpinAnglePointingSet, + _bootstrap_correct_intensity, _complete_pointings, _dps_spin_angles, + _esa_calibration, + _extrapolate_top_intensity, _spin_phase_mask, lo_l2, ) @@ -22,18 +25,27 @@ ANCILLARY_DIR = imap_module_directory / "tests/lo/test_anc" -# A full-spin map, so that every spin-angle bin lands on it. The "ns" after -# "ena" asks for no sputter correction, so these are the uncorrected maps. +# A full-spin map, so that every spin-angle bin lands on it. The "ns" and +# "nbs" after "ena" ask for neither the sputter nor the bootstrap correction, +# so these are the uncorrected maps. FULL_DESCRIPTOR = "l090-enansnbs-h-sf-nsp-full-hae-6deg-3mo" RAM_DESCRIPTOR = "l090-enansnbs-h-sf-nsp-ram-hae-6deg-3mo" -# The same full-spin map, sputter corrected. -SPUTTER_DESCRIPTOR = "l090-enas-h-sf-nsp-full-hae-6deg-3mo" +# The same full-spin map, sputter corrected only. +SPUTTER_DESCRIPTOR = "l090-enasnbs-h-sf-nsp-full-hae-6deg-3mo" + +# The same full-spin map, bootstrap corrected only. +BOOTSTRAP_DESCRIPTOR = "l090-enansbs-h-sf-nsp-full-hae-6deg-3mo" # The contents of imap_lo_sputter-correction-factors-small, as # {target ESA step: {source ESA step: factor}}, 1-based as in the ancillary. SPUTTER_FACTORS = {2: {3: 0.5}, 5: {3: 0.25, 6: 0.1}} +# The contents of imap_lo_bootstrap-correction-factors-small, as +# {target ESA step: {source ESA step: coefficient}}, 1-based as in the +# ancillary. Step 8 is the virtual ESA step above the top of the map. +BOOTSTRAP_FACTORS = {2: {3: 0.4, 5: 0.2}, 6: {7: 0.5}, 7: {8: 0.6}} + N_ESA = LoConstants.N_ESA_LEVELS N_SPIN_BINS = LoConstants.N_SPIN_ANGLE_BINS PIVOT = 90.0 @@ -232,6 +244,24 @@ def sputter_maps(one_pointing, anc_dependencies): return corrected, raw +@pytest.fixture +def bootstrap_maps(one_pointing, anc_dependencies): + """The bootstrap corrected and uncorrected maps of the same pointing. + + The two differ only in whether the correction was applied, so the + uncorrected map supplies the intensities the correction is predicted from. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (corrected,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, BOOTSTRAP_DESCRIPTOR + ) + (raw,) = lo_l2(as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR) + return corrected, raw + + class TestMapStructure: """The shape and contents of the produced map.""" @@ -550,6 +580,211 @@ def test_counts_stay_as_observed(self, sputter_maps): ) +class TestBootstrapCorrection: + """Removing the intensity that bled down from the higher ESA steps.""" + + # The ESA steps the test ancillary corrects from steps of the map itself, + # rather than from the virtual step above the top of it. + MAPPED_SOURCE_STEPS = (2, 6) + + @staticmethod + def bled(intensity, target_esa, power=1): + """The intensity bled into a target ESA step, from the source steps. + + ``power`` is 1 for the intensity itself and 2 for its variance, which + each source term contributes to scaled by the square of its + coefficient. + """ + return sum( + (LoConstants.BOOTSTRAP_SCALE * coefficient) ** power + * intensity[:, source - 1] + for source, coefficient in BOOTSTRAP_FACTORS.get(target_esa, {}).items() + ) + + def test_correction_removes_the_bled_intensity(self, bootstrap_maps): + """A corrected step loses the scaled intensity of the steps above it.""" + corrected, raw = bootstrap_maps + + intensity = raw["ena_intensity"].values + for target_esa in self.MAPPED_SOURCE_STEPS: + expected = np.maximum( + intensity[:, target_esa - 1] - self.bled(intensity, target_esa), 0.0 + ) + np.testing.assert_allclose( + corrected["ena_intensity"].values[:, target_esa - 1], + expected, + rtol=1e-4, + ) + + def test_uncorrected_steps_are_untouched(self, bootstrap_maps): + """A step that nothing bleeds into keeps the intensity it already had.""" + corrected, raw = bootstrap_maps + + untouched = [esa for esa in range(1, N_ESA + 1) if esa not in BOOTSTRAP_FACTORS] + assert untouched, "the test coefficients must leave some steps uncorrected" + + for target_esa in untouched: + np.testing.assert_allclose( + corrected["ena_intensity"].values[:, target_esa - 1], + raw["ena_intensity"].values[:, target_esa - 1], + rtol=1e-6, + ) + + def test_top_step_is_corrected_against_the_virtual_step(self, bootstrap_maps): + """The top step has only the extrapolated step above it to lose to.""" + corrected, raw = bootstrap_maps + + top = N_ESA + assert set(BOOTSTRAP_FACTORS[top]) == {N_ESA + 1}, ( + "the top step must be fed by the virtual step alone" + ) + + correction = ( + raw["ena_intensity"].values[:, top - 1] + - corrected["ena_intensity"].values[:, top - 1] + ) + assert np.all(correction >= 0) + assert np.any(correction > 0), "the virtual step corrected nothing" + + def test_uncertainty_gains_the_source_intensities(self, bootstrap_maps): + """Subtracting a measured quantity can only add to the variance.""" + corrected, raw = bootstrap_maps + + variance = raw["ena_intensity_stat_uncert"].values ** 2 + for target_esa in self.MAPPED_SOURCE_STEPS: + expected = variance[:, target_esa - 1] + self.bled( + variance, target_esa, power=2 + ) + np.testing.assert_allclose( + corrected["ena_intensity_stat_uncert"].values[:, target_esa - 1], + np.sqrt(expected), + rtol=1e-4, + ) + + def test_intensity_is_never_negative(self, bootstrap_maps): + """The corrected map holds no negative intensity.""" + corrected, _ = bootstrap_maps + + assert np.all(corrected["ena_intensity"].values >= 0) + + def test_over_subtraction_is_floored_at_zero(self): + """Over-subtracting a pixel floors it rather than going below zero.""" + calibration = _esa_calibration("h", 0) + intensity = np.ones((1, N_ESA, 4)) + # Every step loses twice its own intensity to the step above it. + coefficients = np.zeros((N_ESA, N_ESA + 1)) + coefficients[np.arange(N_ESA), np.arange(1, N_ESA + 1)] = ( + 2.0 / LoConstants.BOOTSTRAP_SCALE + ) + + corrected, _, plus, minus = _bootstrap_correct_intensity( + intensity, + np.ones_like(intensity), + calibration, + coefficients, + (2, 2), + calibration.geometric_factor_low[:, np.newaxis] > 0, + ) + + np.testing.assert_array_equal(corrected, np.zeros_like(corrected)) + assert np.all(plus >= 0) + assert np.all(minus >= 0) + + def test_systematic_error_brackets_the_correction(self, bootstrap_maps): + """The corrected steps keep a two-sided systematic error.""" + corrected, _ = bootstrap_maps + + plus = corrected["ena_intensity_sys_err_plus"].values + minus = corrected["ena_intensity_sys_err_minus"].values + symmetric = corrected["ena_intensity_sys_err"].values + lit = corrected["ena_intensity"].values > 0 + + assert np.all(plus >= 0) + assert np.all(minus >= 0) + assert np.all(plus[lit] > 0) + np.testing.assert_allclose( + symmetric[lit], np.sqrt(plus[lit] * minus[lit]), rtol=1e-4 + ) + # The correction is bracketed by a smaller and a larger subtraction, so + # its systematic error is wider than the G-factor one it starts from. + for target_esa in self.MAPPED_SOURCE_STEPS: + assert np.any(plus[:, target_esa - 1] > minus[:, target_esa - 1]) + + def test_counts_and_rates_stay_as_observed(self, bootstrap_maps): + """The correction applies to the intensities alone.""" + corrected, raw = bootstrap_maps + + for variable in ("ena_count", "exposure_factor", "ena_count_rate"): + np.testing.assert_array_equal( + corrected[variable].values, raw[variable].values + ) + + +class TestVirtualStepExtrapolation: + """Extrapolating the ESA step above the top of the map.""" + + energy = np.array([0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64]) + grid_shape = (4, 3) + + def intensity(self, second, top): + """A map of the two top ESA steps, one epoch, on a 4x3 sky grid.""" + values = np.zeros((1, N_ESA, np.prod(self.grid_shape))) + values[0, -2] = np.asarray(second, dtype=float).ravel() + values[0, -1] = np.asarray(top, dtype=float).ravel() + return values + + def test_power_law_between_the_top_two_steps(self): + """A pixel with both steps lit extrapolates along its own spectrum.""" + gamma = 1.4 + top = np.full(np.prod(self.grid_shape), 3.0) + second = top * (self.energy[-1] / self.energy[-2]) ** gamma + + extrapolated = _extrapolate_top_intensity( + self.intensity(second, top), self.energy, self.grid_shape + ) + + np.testing.assert_allclose( + extrapolated[0], + top * LoConstants.ESA_8_ENERGY_RATIO**-gamma, + rtol=1e-6, + ) + + def test_missing_spectrum_borrows_from_the_neighborhood(self): + """A pixel with no spectrum of its own uses its neighbors' median.""" + gamma = 1.4 + top = np.full(np.prod(self.grid_shape), 3.0) + second = top * (self.energy[-1] / self.energy[-2]) ** gamma + # One pixel is dark in the second step, so it has no spectrum, but its + # neighbors all share the same one. + second[0] = 0.0 + + extrapolated = _extrapolate_top_intensity( + self.intensity(second, top), self.energy, self.grid_shape + ) + + np.testing.assert_allclose( + extrapolated[0], + top * LoConstants.ESA_8_ENERGY_RATIO**-gamma, + rtol=1e-6, + ) + + def test_map_with_no_spectrum_falls_back_to_the_nominal_index(self): + """With no pixel to learn a spectrum from, a nominal one stands in.""" + top = np.full(np.prod(self.grid_shape), 3.0) + + extrapolated = _extrapolate_top_intensity( + self.intensity(np.zeros_like(top), top), self.energy, self.grid_shape + ) + + np.testing.assert_allclose( + extrapolated[0], + top + * LoConstants.ESA_8_ENERGY_RATIO + ** -LoConstants.BOOTSTRAP_DEFAULT_SPECTRAL_INDEX, + rtol=1e-6, + ) + + class TestGeometry: """The spin-angle to sky-pixel geometry.""" From 87b95f9410a0ca9dedc1253992dbe15cbd027e0d Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 3 Aug 2026 15:15:36 -0400 Subject: [PATCH 07/14] cg correction logic --- imap_processing/lo/constants.py | 11 + imap_processing/lo/l2/lo_l2.py | 310 ++++++++++++++++++++++++- imap_processing/tests/lo/test_lo_l2.py | 250 +++++++++++++++++++- 3 files changed, 562 insertions(+), 9 deletions(-) diff --git a/imap_processing/lo/constants.py b/imap_processing/lo/constants.py index 612da05d9..8d8c8b481 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -115,6 +115,17 @@ class LoConstants: # The spectral index to extrapolate with when the map has none to offer. BOOTSTRAP_DEFAULT_SPECTRAL_INDEX: float = 1.6 + # Compton-Getting correction settings. The energy [eV] a hydrogen ENA has + # in the spacecraft frame purely from the spacecraft's own motion, i.e. + # 1/2 m_H U^2 at the nominal spacecraft speed of ~30 km/s. The kinematics + # of the correction are scaled by it. + CG_ENA_ENERGY_AT_SPACECRAFT_SPEED_EV: float = 4.661 + # The predictor-corrector that estimates the source spectrum behind the + # observed one runs until the RMS change in the intensities falls below the + # tolerance, or the iterations run out. + CG_MAX_ITERATIONS: int = 20 + CG_CONVERGENCE_TOLERANCE: float = 0.005 + # Nominal spin period [s]. True spin duration is NOT 15 seconds. NOMINAL_SPIN_PERIOD_SEC: float = 15.0 diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index f4a9b2e72..5771e460d 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -15,6 +15,7 @@ SkyTilingType, ) from imap_processing.ena_maps.utils.coordinates import CoordNames +from imap_processing.ena_maps.utils.corrections import PowerLawFluxCorrector from imap_processing.ena_maps.utils.naming import MapDescriptor from imap_processing.lo import lo_ancillary from imap_processing.lo.constants import EsaCalibration @@ -89,6 +90,9 @@ def lo_l2( NotImplementedError If a HEALPix map is requested (only rectangular maps supported for Lo), or if the map is of a species other than hydrogen. + ValueError + If the map is to be Compton-Getting corrected but the ancillary + dependencies hold no ESA eta fit factors to correct it with. """ logger.info("Starting IMAP-Lo L2 processing pipeline") @@ -98,10 +102,10 @@ def lo_l2( # Determine which of the corrections the descriptor asks for are needed ( _flux_correction, - _flux_factors, + flux_factors, sputter_correction, bootstrap_correction, - _cg_correction, + cg_correction, ) = _prepare_corrections(map_descriptor, anc_dependencies) logger.info("Step 1: Loading ancillary data") @@ -134,8 +138,10 @@ def lo_l2( if sputter_correction else (None, None) ) - accumulators = ACCUMULATED_VARIABLES + ( - ("sputter_source_count",) if sputter_source else () + accumulators = ( + ACCUMULATED_VARIABLES + + (("sputter_source_count",) if sputter_source else ()) + + (("cos_alpha_exposure",) if cg_correction else ()) ) _initialize_accumulators(sky_map, calibration.energy, accumulators) @@ -150,12 +156,25 @@ def lo_l2( map_descriptor, calibration.energy, sputter_source, + cg_correction, ) bootstrap_matrix = _bootstrap_correction() if bootstrap_correction else None + # The Compton-Getting correction reads the source spectrum of each pixel + # through the ESA transmission factors of the eta fit ancillary. + flux_corrector = None + if cg_correction: + if flux_factors is None: + raise ValueError( + "A heliospheric frame map needs the ESA eta fit factors to be " + "Compton-Getting corrected, and none were found in the " + "ancillary dependencies" + ) + flux_corrector = PowerLawFluxCorrector(flux_factors) + variables = _calculate_rates_and_intensities( - sky_map, calibration, sputter_matrix, bootstrap_matrix + sky_map, calibration, sputter_matrix, bootstrap_matrix, flux_corrector ) dataset = _build_map_dataset(sky_map, variables, calibration) @@ -610,6 +629,7 @@ def _accumulate_pointing( map_descriptor: MapDescriptor, energy: np.ndarray, sputter_source: str | None = None, + accumulate_cos_alpha: bool = False, ) -> None: """ Add one pointing's counts and exposure to the map. @@ -633,6 +653,10 @@ def _accumulate_pointing( sputter_source : str | None The species sputtering into the mapped species, whose counts are accumulated alongside. None if this map is not sputter corrected. + accumulate_cos_alpha : bool + Whether to accumulate the RAM projection of each spin-angle bin, which + the Compton-Getting correction needs. False if this map is not CG + corrected. """ species = map_descriptor.species pivot_angle = float(np.atleast_1d(goodtimes["pivot"].values)[0]) @@ -673,6 +697,13 @@ def _accumulate_pointing( values["sputter_source_count"] = ( histrates[f"{sputter_source}_counts"].values[in_goodtime].sum(axis=0) ) + if accumulate_cos_alpha: + # The RAM projection is a property of the bin, not of the counts in it, + # so like the background it is accumulated weighted by exposure and + # divided by the total exposure at the end. + values["cos_alpha_exposure"] = ( + _ram_projection(spin_angles, pivot_angle) * pointing_exposure + ) pointing_set = LoSpinAnglePointingSet( epoch, @@ -749,10 +780,35 @@ def _spin_phase_mask( f"Invalid spin phase: {spin_phase}. Must be 'ram', 'anti' or 'full'." ) - ram_projection = np.sin(np.radians(pivot_angle + c.PIVOT_RAM_OFFSET)) * np.sin( + ram_projection = _ram_projection(spin_angles, pivot_angle) + return ram_projection > 0 if spin_phase == "ram" else ram_projection < 0 + + +def _ram_projection(spin_angles: np.ndarray, pivot_angle: float) -> np.ndarray: + """ + Project the look direction of each spin-angle bin onto the RAM direction. + + The projection is the cosine of the angle between the bin's look direction + and the spacecraft's velocity, which is what tells the RAM half of the spin + from the anti-RAM half, and what the Compton-Getting correction is a + function of. + + Parameters + ---------- + spin_angles : np.ndarray + The IMAP_DPS azimuth [degrees] of each spin-angle bin. + pivot_angle : float + The pivot angle [degrees] of the pointing. + + Returns + ------- + np.ndarray + The projection factor of each bin, positive looking into the RAM + direction and negative looking away from it. + """ + return np.sin(np.radians(pivot_angle + c.PIVOT_RAM_OFFSET)) * np.sin( np.radians(spin_angles) ) - return ram_projection > 0 if spin_phase == "ram" else ram_projection < 0 # ============================================================================= @@ -1115,11 +1171,211 @@ def subtract(scale: float) -> np.ndarray: ) +def _power_law_slopes(intensity: np.ndarray, energy: np.ndarray) -> np.ndarray: + """ + Estimate the spectral index of every ESA level of every pixel. + + Parameters + ---------- + intensity : np.ndarray + The intensity of each ESA level, of shape (epoch, esa level, pixel), + NaN where a pixel saw nothing at a level. + energy : np.ndarray + The energy of each ESA level, in any unit. + + Returns + ------- + np.ndarray + The index of the power law through each level and the one above it, of + the same shape, NaN where either level is unmeasured. The top level has + no level above it, so it keeps the index of the one below it. + """ + with np.errstate(divide="ignore", invalid="ignore"): + slopes = ( + np.log(intensity[:, 1:] / intensity[:, :-1]) + / np.log(energy[1:] / energy[:-1])[:, np.newaxis] + ) + return np.concatenate([slopes, slopes[:, -1:]], axis=1) + + +def _source_intensity( + intensity: np.ndarray, + energy: np.ndarray, + flux_corrector: PowerLawFluxCorrector, +) -> tuple[np.ndarray, np.ndarray]: + """ + Undo the ESA transmission bias in the observed intensities. + + An ESA level integrates over a passband rather than sampling a single + energy, so what it observes depends on the spectrum falling through it. The + transmission factor that relates the two is itself a function of the + spectral index, so the source spectrum is recovered by iterating: estimate + the index, undo the transmission, re-estimate the index, until the + intensities settle. + + Parameters + ---------- + intensity : np.ndarray + The observed intensity of each ESA level, of shape + (epoch, esa level, pixel). + energy : np.ndarray + The energy of each ESA level, in any unit. + flux_corrector : PowerLawFluxCorrector + The ESA transmission factors, read from the eta fit ancillary. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + The source intensity and its spectral index, both NaN in the pixels + that saw nothing. + """ + levels: np.ndarray = np.arange(c.N_ESA_LEVELS) + 1 + + def transmission(spectral_index: np.ndarray) -> np.ndarray: + """ + Get the ESA transmission factor of each level of each pixel. + + Parameters + ---------- + spectral_index : np.ndarray + The index of the power law through each ESA level of each pixel. + + Returns + ------- + np.ndarray + The transmission factor, of the same shape. + """ + # The shared corrector takes the ESA level on the leading axis. + return np.moveaxis( + flux_corrector.eta_esa(levels, np.moveaxis(spectral_index, 1, 0)), 0, 1 + ) + + # A pixel that saw nothing at a level has no spectrum through it; the NaN + # carries through the iteration and out to the corrected map. + observed = np.where(intensity > 0, intensity, np.nan) + + index = _power_law_slopes(observed, energy) + source = observed / transmission(index) + + for iteration in range(c.CG_MAX_ITERATIONS): + predicted = 0.5 * (index + _power_law_slopes(source, energy)) + corrected = 0.5 * ( + index + _power_law_slopes(observed / transmission(predicted), energy) + ) + + previous, source = source, observed / transmission(corrected) + index = corrected + + with np.errstate(divide="ignore", invalid="ignore"): + change = np.sqrt(np.nanmean((source / previous) ** 2)) - 1.0 + if np.isfinite(change) and abs(change) < c.CG_CONVERGENCE_TOLERANCE: + logger.debug(f"Source spectrum converged after {iteration + 1} iterations") + break + else: + logger.warning( + f"Source spectrum did not converge in {c.CG_MAX_ITERATIONS} iterations" + ) + + return source, index + + +def _spacecraft_frame_energy(cos_alpha: np.ndarray, energy: np.ndarray) -> np.ndarray: + """ + Get the energy an ENA of a given heliospheric energy arrives with. + + An ENA arriving at the spacecraft is seen at a different energy than it has + in the heliosphere, by the spacecraft's own motion through it: the same ENA + is faster in the spacecraft frame when the spacecraft is moving into it and + slower when it is moving away. + + Parameters + ---------- + cos_alpha : np.ndarray + The cosine of the angle between the look direction of each pixel and + the spacecraft's velocity, of shape (epoch, esa level, pixel). + energy : np.ndarray + The heliospheric-frame energy [eV] of each ESA level. + + Returns + ------- + np.ndarray + The spacecraft-frame energy [eV] of each ESA level of each pixel, of + the shape of ``cos_alpha``. + """ + energy_u = c.CG_ENA_ENERGY_AT_SPACECRAFT_SPEED_EV + cos_alpha = np.clip(cos_alpha, -1.0, 1.0) + + # The speed of the ENA in the spacecraft frame, in units of the + # spacecraft's own speed: x = cos(a) + sqrt(y^2 - sin^2(a)), y^2 = E / E_u. + ratio = (energy / energy_u)[:, np.newaxis] + speed = cos_alpha + np.sqrt(np.maximum(ratio + cos_alpha**2 - 1.0, 0.0)) + return speed**2 * energy_u + + +def _compton_getting_correct_intensity( + intensity: np.ndarray, + uncertainties: tuple[np.ndarray, ...], + cos_alpha: np.ndarray, + calibration: EsaCalibration, + flux_corrector: PowerLawFluxCorrector, +) -> tuple[np.ndarray, tuple[np.ndarray, ...]]: + """ + Move the intensities from the spacecraft frame into the heliospheric one. + + The intensity a pixel reports is of ENAs at their spacecraft-frame energy, + which the spacecraft's motion has shifted away from the heliospheric-frame + energy the map is binned in. The correction reads the source spectrum of + the pixel off its own power law at the shifted energy, and scales the + intensity back onto the map's energy. + + Parameters + ---------- + intensity : np.ndarray + The intensity of every ESA level, of shape (epoch, esa level, pixel). + uncertainties : tuple[np.ndarray, ...] + The uncertainties on that intensity, each of the same shape. They are + scaled by the same factor as the intensity they belong to. + cos_alpha : np.ndarray + The cosine of the angle between the look direction of each pixel and + the spacecraft's velocity, same shape. + calibration : EsaCalibration + The energy response the map is binned in, read for the energies the + correction shifts between. + flux_corrector : PowerLawFluxCorrector + The ESA transmission factors, read from the eta fit ancillary. + + Returns + ------- + tuple[np.ndarray, tuple[np.ndarray, ...]] + The corrected intensity and its correspondingly scaled uncertainties, + zero in the pixels the correction has nothing to say about. + """ + # The kinematics are in eV; the map is binned in keV. + energy = calibration.energy * 1e3 + + source, spectral_index = _source_intensity(intensity, energy, flux_corrector) + energy_sc = _spacecraft_frame_energy(cos_alpha, energy) + + with np.errstate(divide="ignore", invalid="ignore"): + corrected = source * (energy_sc / energy[:, np.newaxis]) ** ( + spectral_index + 1.0 + ) + # The uncertainties are fractionally unchanged, so they move with the + # intensity. This folds in the transmission factor as well, which the + # source intensity was already divided by. + scaling = np.where(intensity > 0, corrected / intensity, np.nan) + + return np.nan_to_num(corrected), tuple( + np.nan_to_num(uncertainty * scaling) for uncertainty in uncertainties + ) + + def _calculate_rates_and_intensities( sky_map: RectangularSkyMap, calibration: EsaCalibration, sputter_matrix: np.ndarray | None = None, bootstrap_matrix: np.ndarray | None = None, + flux_corrector: PowerLawFluxCorrector | None = None, ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -1141,6 +1397,10 @@ def _calculate_rates_and_intensities( The (target level, source level) bootstrap correction coefficients, applied to the intensities once they are derived. None leaves the intensities as the counts gave them. + flux_corrector : PowerLawFluxCorrector | None + The ESA transmission factors the Compton-Getting correction recovers + the source spectrum with. None leaves the intensities in the + spacecraft frame. Returns ------- @@ -1234,6 +1494,39 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: bg_intensity = _divide(bg_rate, geometric_factor * energy) bg_intensity_stat_uncert = _divide(bg_rate_stat_uncert, geometric_factor * energy) + # The Compton-Getting correction comes last, moving the intensities out of + # the frame the instrument observed them in. The background is a spectrum + # of its own, so it is corrected on its own terms rather than with the + # map's. + if flux_corrector is not None: + cos_alpha = _divide(sky_map.data_1d["cos_alpha_exposure"].values, exposure) + logger.info("Applying the Compton-Getting correction to the intensities") + ( + intensity, + ( + intensity_stat_uncert, + intensity_sys_err_plus, + intensity_sys_err_minus, + ), + ) = _compton_getting_correct_intensity( + intensity, + ( + intensity_stat_uncert, + intensity_sys_err_plus, + intensity_sys_err_minus, + ), + cos_alpha, + calibration, + flux_corrector, + ) + bg_intensity, (bg_intensity_stat_uncert,) = _compton_getting_correct_intensity( + bg_intensity, + (bg_intensity_stat_uncert,), + cos_alpha, + calibration, + flux_corrector, + ) + return { "ena_count": counts, "exposure_factor": exposure, @@ -1285,7 +1578,8 @@ def _build_map_dataset( sky_map.data_1d[name] = xr.DataArray(values.astype(np.float32), dims=dims) # These are accumulators, not map variables. sky_map.data_1d = sky_map.data_1d.drop_vars( - ["bg_rate_exposure", "sputter_source_count"], errors="ignore" + ["bg_rate_exposure", "sputter_source_count", "cos_alpha_exposure"], + errors="ignore", ) dataset = sky_map.to_dataset() diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 92fabbe5e..23729d70c 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -9,15 +9,18 @@ from imap_processing import imap_module_directory from imap_processing.cdf.utils import load_cdf, write_cdf from imap_processing.ena_maps.ena_maps import match_coords_to_indices +from imap_processing.ena_maps.utils.corrections import PowerLawFluxCorrector from imap_processing.ena_maps.utils.naming import MapDescriptor -from imap_processing.lo.constants import LoConstants +from imap_processing.lo.constants import EsaCalibration, LoConstants from imap_processing.lo.l2.lo_l2 import ( LoSpinAnglePointingSet, _bootstrap_correct_intensity, _complete_pointings, + _compton_getting_correct_intensity, _dps_spin_angles, _esa_calibration, _extrapolate_top_intensity, + _spacecraft_frame_energy, _spin_phase_mask, lo_l2, ) @@ -37,6 +40,10 @@ # The same full-spin map, bootstrap corrected only. BOOTSTRAP_DESCRIPTOR = "l090-enansbs-h-sf-nsp-full-hae-6deg-3mo" +# The same map in the heliospheric frame, which is what asks for the +# Compton-Getting correction. Neither of the other corrections is made. +CG_DESCRIPTOR = "l090-enansnbs-h-hf-nsp-full-hae-6deg-3mo" + # The contents of imap_lo_sputter-correction-factors-small, as # {target ESA step: {source ESA step: factor}}, 1-based as in the ancillary. SPUTTER_FACTORS = {2: {3: 0.5}, 5: {3: 0.25, 6: 0.1}} @@ -262,6 +269,25 @@ def bootstrap_maps(one_pointing, anc_dependencies): return corrected, raw +@pytest.fixture +def cg_maps(one_pointing, anc_dependencies): + """The Compton-Getting corrected and uncorrected maps of one pointing. + + The two differ only in the frame they are made in, which is what asks for + the correction, so the spacecraft frame map supplies the intensities the + correction is predicted from. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (corrected,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, CG_DESCRIPTOR + ) + (raw,) = lo_l2(as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR) + return corrected, raw + + class TestMapStructure: """The shape and contents of the produced map.""" @@ -785,6 +811,228 @@ def test_map_with_no_spectrum_falls_back_to_the_nominal_index(self): ) +class TestComptonGettingCorrection: + """Moving the intensities into the frame the heliosphere sees them in.""" + + def test_counts_and_rates_stay_as_observed(self, cg_maps): + """The correction applies to the intensities alone.""" + corrected, raw = cg_maps + + for variable in ("ena_count", "exposure_factor", "ena_count_rate"): + np.testing.assert_array_equal( + corrected[variable].values, raw[variable].values + ) + + def test_intensities_are_corrected(self, cg_maps): + """The intensities and the background come out shifted.""" + corrected, raw = cg_maps + + for variable in ("ena_intensity", "bg_intensity"): + lit = raw[variable].values > 0 + assert lit.any() + assert np.any( + ~np.isclose(corrected[variable].values[lit], raw[variable].values[lit]) + ), f"{variable} was not corrected" + + def test_map_holds_no_fill_values(self, cg_maps): + """A pixel the correction says nothing about comes out at zero.""" + corrected, _ = cg_maps + + for variable in ( + "ena_intensity", + "ena_intensity_stat_uncert", + "ena_intensity_sys_err_plus", + "ena_intensity_sys_err_minus", + "bg_intensity", + "bg_intensity_stat_uncert", + ): + values = corrected[variable].values + assert np.isfinite(values).all(), f"{variable} holds NaN or inf" + assert np.all(values >= 0), f"{variable} went negative" + + def test_dark_pixels_stay_dark(self, cg_maps): + """A pixel that saw nothing has no spectrum to correct.""" + corrected, raw = cg_maps + + dark = raw["ena_intensity"].values <= 0 + assert dark.any() + assert np.all(corrected["ena_intensity"].values[dark] == 0) + + def test_accumulator_is_not_written_out(self, cg_maps): + """The projected cos(alpha) is an accumulator, not a map variable.""" + corrected, _ = cg_maps + + assert "cos_alpha_exposure" not in corrected.data_vars + + +class TestSpacecraftFrameEnergy: + """The kinematics of an ENA seen from a moving spacecraft.""" + + # The energy [eV] an ENA of the top ESA level has in the helio frame. + energy = np.array([16.0, 30.0, 56.0, 106.0, 200.0, 404.0, 787.0]) + + def spacecraft_energy(self, cos_alpha): + """The spacecraft frame energies of one pixel at every ESA level.""" + alpha = np.full((1, N_ESA, 1), float(cos_alpha)) + return _spacecraft_frame_energy(alpha, self.energy)[0, :, 0] + + def test_head_on_ena_gains_the_spacecraft_speed(self): + """Looking into the ram direction, the speeds add.""" + energy_u = LoConstants.CG_ENA_ENERGY_AT_SPACECRAFT_SPEED_EV + + np.testing.assert_allclose( + self.spacecraft_energy(1.0), + (np.sqrt(self.energy) + np.sqrt(energy_u)) ** 2, + rtol=1e-12, + ) + + def test_overtaken_ena_loses_the_spacecraft_speed(self): + """Looking away from the ram direction, the speeds subtract.""" + energy_u = LoConstants.CG_ENA_ENERGY_AT_SPACECRAFT_SPEED_EV + + np.testing.assert_allclose( + self.spacecraft_energy(-1.0), + (np.sqrt(self.energy) - np.sqrt(energy_u)) ** 2, + rtol=1e-12, + ) + + def test_side_on_ena_loses_the_spacecraft_energy(self): + """Looking across the ram direction, the energies subtract.""" + energy_u = LoConstants.CG_ENA_ENERGY_AT_SPACECRAFT_SPEED_EV + + np.testing.assert_allclose( + self.spacecraft_energy(0.0), self.energy - energy_u, rtol=1e-12 + ) + + def test_ram_side_is_the_energetic_one(self): + """The shift grows monotonically towards the ram direction.""" + cosines = np.linspace(-1.0, 1.0, 21) + energies = np.array([self.spacecraft_energy(c) for c in cosines]) + + assert np.all(np.diff(energies, axis=0) > 0) + + +class TestComptonGettingMaths: + """The correction applied to a spectrum whose answer is known.""" + + energy_kev = np.array([0.016, 0.030, 0.056, 0.106, 0.200, 0.404, 0.787]) + spectral_index = -1.8 + + @pytest.fixture + def calibration(self): + """A calibration carrying only the energies the correction reads.""" + ones = np.ones(N_ESA) + return EsaCalibration( + energy=self.energy_kev, + energy_delta_minus=np.zeros(N_ESA), + energy_delta_plus=np.zeros(N_ESA), + geometric_factor=ones, + geometric_factor_low=ones, + geometric_factor_high=ones, + ) + + @pytest.fixture + def flux_corrector(self): + """The ESA transmission factors of the test eta fit ancillary.""" + return PowerLawFluxCorrector( + ANCILLARY_DIR / "imap_lo_esa-eta-fit-factors_20240101_v001.csv" + ) + + @pytest.fixture + def transparent_corrector(self, tmp_path): + """A corrector whose every ESA level transmits perfectly. + + With no transmission to undo, the source spectrum is the observed one + and the correction is the energy shift alone, which is what makes the + expected values here analytic. + """ + coefficients = tmp_path / "imap_lo_esa-eta-fit-factors_20240101_v999.csv" + rows = "\n".join(f"{step},1,0,0,0,0,0" for step in range(1, N_ESA + 1)) + coefficients.write_text(f"esa_step,M0,M1,M2,M3,M4,M5\n{rows}\n") + return PowerLawFluxCorrector(coefficients) + + def power_law(self, n_pixels): + """A spectrum of exactly the test spectral index, at every pixel.""" + spectrum = 1e5 * (self.energy_kev / self.energy_kev[0]) ** self.spectral_index + return np.tile(spectrum[np.newaxis, :, np.newaxis], (1, 1, n_pixels)) + + def test_correction_follows_the_source_power_law( + self, calibration, transparent_corrector + ): + """A power-law spectrum is scaled by the shift along its own slope.""" + cos_alpha = np.linspace(-1.0, 1.0, 12)[np.newaxis, np.newaxis, :] + cos_alpha = np.tile(cos_alpha, (1, N_ESA, 1)) + intensity = self.power_law(cos_alpha.shape[-1]) + + corrected, _ = _compton_getting_correct_intensity( + intensity, (), cos_alpha, calibration, transparent_corrector + ) + + # The index of a pure power law is recovered exactly, and a perfectly + # transmitting ESA leaves the source intensity as the observed one. + energy_ev = self.energy_kev * 1e3 + energy_sc = _spacecraft_frame_energy(cos_alpha, energy_ev) + expected = intensity * (energy_sc / energy_ev[:, np.newaxis]) ** ( + self.spectral_index + 1.0 + ) + + np.testing.assert_allclose(corrected, expected, rtol=1e-10) + + def test_transmission_is_divided_out( + self, calibration, flux_corrector, transparent_corrector + ): + """An ESA that over-transmits has to be corrected back down.""" + cos_alpha = np.full((1, N_ESA, 4), 0.5) + intensity = self.power_law(4) + + corrected, _ = _compton_getting_correct_intensity( + intensity, (), cos_alpha, calibration, flux_corrector + ) + transparent, _ = _compton_getting_correct_intensity( + intensity, (), cos_alpha, calibration, transparent_corrector + ) + + # The transmission of a falling spectrum is above one at every level, + # so dividing it out leaves less intensity than a perfect ESA would. + index = np.full((N_ESA, 1), self.spectral_index) + assert np.all(flux_corrector.eta_esa(np.arange(N_ESA) + 1, index) > 1.0) + assert np.all(corrected < transparent) + + def test_uncertainties_keep_their_fraction(self, calibration, flux_corrector): + """Every uncertainty moves by the same factor as its intensity.""" + cos_alpha = np.tile( + np.linspace(-1.0, 1.0, 12)[np.newaxis, np.newaxis, :], (1, N_ESA, 1) + ) + intensity = self.power_law(cos_alpha.shape[-1]) + stat_uncert = 0.1 * intensity + sys_err = 0.05 * intensity + + corrected, (corrected_stat, corrected_sys) = _compton_getting_correct_intensity( + intensity, + (stat_uncert, sys_err), + cos_alpha, + calibration, + flux_corrector, + ) + + np.testing.assert_allclose(corrected_stat, 0.1 * corrected, rtol=1e-10) + np.testing.assert_allclose(corrected_sys, 0.05 * corrected, rtol=1e-10) + + def test_unlit_levels_are_left_at_zero(self, calibration, flux_corrector): + """A level with no intensity has no spectrum, and stays empty.""" + cos_alpha = np.full((1, N_ESA, 4), 0.5) + intensity = self.power_law(4) + intensity[0, :, 2] = 0.0 + + corrected, (corrected_stat,) = _compton_getting_correct_intensity( + intensity, (intensity,), cos_alpha, calibration, flux_corrector + ) + + np.testing.assert_array_equal(corrected[0, :, 2], np.zeros(N_ESA)) + np.testing.assert_array_equal(corrected_stat[0, :, 2], np.zeros(N_ESA)) + assert np.all(corrected[0, :, [0, 1, 3]] > 0) + + class TestGeometry: """The spin-angle to sky-pixel geometry.""" From 1674dea78e8ec35d25449cd6670825d37445ca87 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Tue, 4 Aug 2026 09:46:02 -0400 Subject: [PATCH 08/14] added tests to bump coverage --- imap_processing/tests/lo/test_lo_l2.py | 94 ++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index c8ac878fb..51321bd05 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -11,12 +11,18 @@ from imap_processing.ena_maps.ena_maps import match_coords_to_indices from imap_processing.ena_maps.utils.naming import MapDescriptor from imap_processing.lo.constants import LoConstants +from imap_processing.lo.l2.lo_l2 import ( + ANCILLARY_DATA_DIR as PACKAGE_ANCILLARY_DIR, +) from imap_processing.lo.l2.lo_l2 import ( LoSpinAnglePointingSet, _complete_pointings, _dps_spin_angles, _spin_phase_mask, + finalize_dataset, lo_l2, + load_bootstrap_correction_data, + load_sputter_correction_data, ) from imap_processing.spice.time import met_to_ttj2000ns @@ -189,6 +195,17 @@ def anc_dependencies(): return [ANCILLARY_DIR / "imap_lo_esa-eta-fit-factors_20240101_v001.csv"] +@pytest.fixture +def shipped_ancillaries(): + """Read the calibration ancillaries shipped with the package. + + Undoes the autouse ``use_test_geometric_factors`` patch (a more common use-case) + for tests that need the shipped ancillaries. + """ + with patch("imap_processing.lo.l2.lo_l2.ANCILLARY_DATA_DIR", PACKAGE_ANCILLARY_DIR): + yield + + @pytest.fixture def full_map(one_pointing, anc_dependencies): """The full-spin map of one pointing, with the sky pointing mocked.""" @@ -505,6 +522,83 @@ def test_missing_product_raises(self, one_pointing): _complete_pointings(dependencies) +class TestCorrectionFactors: + """Reading the sputter and bootstrap correction ancillaries.""" + + def test_sputter_factors_are_selected_by_species_pair( + self, shipped_ancillaries, tmp_path + ): + """Only the rows of the requested pair come back, in ESA step order.""" + factors = load_sputter_correction_data("o", "h") + + assert list(factors.columns) == [ + "source_species", + "target_species", + "esa_step", + "sputter_factor", + "sputter_factor_uncertainty", + ] + assert not factors.empty + assert (factors["source_species"] == "o").all() + assert (factors["target_species"] == "h").all() + + with patch("imap_processing.lo.l2.lo_l2.ANCILLARY_DATA_DIR", tmp_path): + with pytest.raises(ValueError, match="No sputter correction files"): + load_sputter_correction_data("o", "h") + + def test_bootstrap_factors_relate_lower_steps_to_higher( + self, shipped_ancillaries, tmp_path + ): + """Each factor carries a step pair, the source below the target.""" + factors = load_bootstrap_correction_data() + + assert list(factors.columns) == [ + "esa_step_i", + "esa_step_k", + "bootstrap_factor", + ] + assert not factors.empty + assert (factors["esa_step_i"] < factors["esa_step_k"]).all() + # Steps are 1-based, and step 8 is the virtual E8 channel. + assert factors["esa_step_i"].min() >= 1 + assert factors["esa_step_k"].max() <= N_ESA + 1 + assert (factors["bootstrap_factor"] > 0).all() + + with patch("imap_processing.lo.l2.lo_l2.ANCILLARY_DATA_DIR", tmp_path): + with pytest.raises(ValueError, match="No bootstrap correction factor"): + load_bootstrap_correction_data() + + +class TestFinalizeDataset: + """Attaching the CDF attributes to a finished map.""" + + def test_attributes_are_filled_in_from_the_descriptor(self): + """The map is labelled with its descriptor and its variables described.""" + dataset = xr.Dataset( + { + "ena_intensity": ( + ["epoch", "energy"], + np.zeros((1, N_ESA)), + ), + "not_a_map_variable": ("epoch", np.zeros(1)), + }, + coords={"epoch": [0], "energy": ESA_ENERGIES[0]}, + ) + + finalized = finalize_dataset(dataset, FULL_DESCRIPTOR) + + assert finalized.attrs["Logical_source"] == f"imap_lo_l2_{FULL_DESCRIPTOR}" + assert FULL_DESCRIPTOR in finalized.attrs["Data_type"] + + # A known map variable picks up its attributes from the enamaps config. + assert finalized["ena_intensity"].attrs["FIELDNAM"] == "Intensity" + assert finalized["ena_intensity"].attrs["UNITS"] == "cm -2 s -1 sr -1 keV -1" + + # A variable the config says nothing about is left without attributes, + # rather than failing the map. + assert finalized["not_a_map_variable"].attrs == {} + + class TestUnsupported: """Map flavours the Lo pipeline does not make.""" From 45d20bc613d58f033176a9b5ffe6eacd96554d7c Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Tue, 4 Aug 2026 14:17:53 -0400 Subject: [PATCH 09/14] explanatory note to a docstring --- imap_processing/lo/l2/lo_l2.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 5771e460d..085db14ec 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -1040,6 +1040,29 @@ def _extrapolate_top_intensity( np.ndarray The intensity of the virtual level, of shape (epoch, pixel), zero in the pixels the top level saw nothing in. + + Notes + ----- + The spectrum is taken as a power law in energy, ``I(E) = A * E ** -gamma``, + where ``gamma`` is the spectral index. The normalization ``A`` never has to + be evaluated: writing the law at the top two levels and dividing one by the + other cancels it, leaving:: + + I_top / I_second = (E_top / E_second) ** -gamma + + Taking logs of both sides gives the index the code solves for:: + + log(I_top / I_second) = -gamma * log(E_top / E_second) + gamma = -log(I_top / I_second) / log(E_top / E_second) + + The virtual level is then read off the same law anchored on the top level, + which keeps ``A`` cancelled, so a pixel needs only a spectral index and one + measured point to be extrapolated:: + + I_virtual = I_top * (E_virtual / E_top) ** -gamma + + That is what lets a pixel with no spectrum of its own borrow an index from + its neighbors and still use its own top level as the anchor. """ second, top = energy[-2], energy[-1] virtual = top * c.ESA_8_ENERGY_RATIO From 6e4a22fb5c4ccd6c2a310378f404475914463523 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Tue, 4 Aug 2026 15:55:08 -0400 Subject: [PATCH 10/14] ISN masking --- imap_processing/ena_maps/utils/naming.py | 30 ++ .../imap_lo_isn-mask-parameters_v001.csv | 22 ++ imap_processing/lo/l2/lo_l2.py | 292 ++++++++++++------ imap_processing/tests/ena_maps/test_naming.py | 46 +++ ...imap_lo_isn-mask-parameters-small_v001.csv | 22 ++ imap_processing/tests/lo/test_lo_l2.py | 203 +++++++++++- 6 files changed, 521 insertions(+), 94 deletions(-) create mode 100644 imap_processing/lo/ancillary_data/imap_lo_isn-mask-parameters_v001.csv create mode 100644 imap_processing/tests/lo/test_anc/imap_lo_isn-mask-parameters-small_v001.csv diff --git a/imap_processing/ena_maps/utils/naming.py b/imap_processing/ena_maps/utils/naming.py index f05aec6fe..41fab6e5c 100644 --- a/imap_processing/ena_maps/utils/naming.py +++ b/imap_processing/ena_maps/utils/naming.py @@ -374,6 +374,36 @@ def bootstrap_corrected(self) -> bool: # The bootstrap code follows the sputter one, see sputter_corrected. return re.match(r"(?:n?s)?bs", self.principal_data_extras) is not None + @property + def cg_corrected(self) -> bool: + """ + Whether the map's intensities were moved into the heliospheric frame. + + Returns + ------- + cg_corrected : bool + True if this map is Compton-Getting corrected. + """ + return self.frame_descriptor == "hf" + + @property + def isn_masked(self) -> bool: + """ + Whether the map has the interstellar neutral band masked out. + + Returns + ------- + isn_masked : bool + True if this map is ISN masked. + """ + if not self.principal_data.startswith("ena") or self.raw: + return False + # The mask code follows the bootstrap one, see bootstrap_corrected. + return ( + re.match(r"(?:n?s)?n?bsmsk", self.principal_data_extras, re.IGNORECASE) + is not None + ) + # Methods for parsing and building parts of the map descriptor string @staticmethod def get_instrument_descriptor( diff --git a/imap_processing/lo/ancillary_data/imap_lo_isn-mask-parameters_v001.csv b/imap_processing/lo/ancillary_data/imap_lo_isn-mask-parameters_v001.csv new file mode 100644 index 000000000..e0eb565a0 --- /dev/null +++ b/imap_processing/lo/ancillary_data/imap_lo_isn-mask-parameters_v001.csv @@ -0,0 +1,22 @@ +pivot_angle,esa_step,intensity_threshold_fraction,angular_width_deg,outlier_percentile +75,1,0.015,80.0,99.999 +75,2,0.015,70.0,99.999 +75,3,0.015,70.0,99.999 +75,4,0.008,70.0,99.999 +75,5,0.02,25.0,96.999 +75,6,1.0,35.0,99.999 +75,7,1.0,35.0,99.999 +90,1,0.01,80.0,99.999 +90,2,0.01,70.0,99.999 +90,3,0.005,65.0,99.999 +90,4,0.005,50.0,99.999 +90,5,0.1,45.0,99.999 +90,6,1.0,35.0,99.999 +90,7,1.0,35.0,99.999 +105,1,0.01,80.0,99.999 +105,2,0.005,60.0,99.999 +105,3,0.001,60.0,99.999 +105,4,0.005,40.0,99.999 +105,5,0.15,35.0,99.999 +105,6,1.0,35.0,99.999 +105,7,1.0,35.0,99.999 diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 085db14ec..2c3a6d94b 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -40,6 +40,19 @@ # or intensities are derived from them. ACCUMULATED_VARIABLES = ("ena_count", "exposure_factor", "bg_rate_exposure") +# The map variables the ISN mask blanks out: everything derived from the counts +# of the mapped species. +ISN_MASKED_VARIABLES = ( + "ena_count", + "ena_count_rate", + "ena_count_rate_stat_uncert", + "ena_intensity", + "ena_intensity_stat_uncert", + "ena_intensity_sys_err", + "ena_intensity_sys_err_plus", + "ena_intensity_sys_err_minus", +) + # The calibration ancillaries shipped with the package. ANCILLARY_DATA_DIR = Path(__file__).parent.parent / "ancillary_data" @@ -92,21 +105,25 @@ def lo_l2( or if the map is of a species other than hydrogen. ValueError If the map is to be Compton-Getting corrected but the ancillary - dependencies hold no ESA eta fit factors to correct it with. + dependencies hold no ESA eta fit factors to correct it with, or if it + is to be ISN masked at a pivot angle the mask has no tuning for. """ logger.info("Starting IMAP-Lo L2 processing pipeline") map_descriptor = MapDescriptor.from_string(descriptor) logger.info(f"Processing map for species: {map_descriptor.species}") - # Determine which of the corrections the descriptor asks for are needed - ( - _flux_correction, - flux_factors, - sputter_correction, - bootstrap_correction, - cg_correction, - ) = _prepare_corrections(map_descriptor, anc_dependencies) + # Read the ancillaries the corrections asked for by the descriptor need up + # front, so that a map missing one fails before any of its pointings are + # accumulated. + flux_corrector = ( + _flux_corrector(anc_dependencies) if map_descriptor.cg_corrected else None + ) + isn_mask_parameters = ( + _isn_mask_parameters(int(map_descriptor.sensor)) + if map_descriptor.isn_masked + else None + ) logger.info("Step 1: Loading ancillary data") _efficiency_data = load_efficiency_data(anc_dependencies) @@ -135,13 +152,13 @@ def lo_l2( # same grid, alongside the map's own. sputter_source, sputter_matrix = ( _sputter_correction(map_descriptor.species) - if sputter_correction + if map_descriptor.sputter_corrected else (None, None) ) accumulators = ( ACCUMULATED_VARIABLES + (("sputter_source_count",) if sputter_source else ()) - + (("cos_alpha_exposure",) if cg_correction else ()) + + (("cos_alpha_exposure",) if map_descriptor.cg_corrected else ()) ) _initialize_accumulators(sky_map, calibration.energy, accumulators) @@ -156,25 +173,19 @@ def lo_l2( map_descriptor, calibration.energy, sputter_source, - cg_correction, ) - bootstrap_matrix = _bootstrap_correction() if bootstrap_correction else None - - # The Compton-Getting correction reads the source spectrum of each pixel - # through the ESA transmission factors of the eta fit ancillary. - flux_corrector = None - if cg_correction: - if flux_factors is None: - raise ValueError( - "A heliospheric frame map needs the ESA eta fit factors to be " - "Compton-Getting corrected, and none were found in the " - "ancillary dependencies" - ) - flux_corrector = PowerLawFluxCorrector(flux_factors) + bootstrap_matrix = ( + _bootstrap_correction() if map_descriptor.bootstrap_corrected else None + ) variables = _calculate_rates_and_intensities( - sky_map, calibration, sputter_matrix, bootstrap_matrix, flux_corrector + sky_map, + calibration, + sputter_matrix, + bootstrap_matrix, + flux_corrector, + isn_mask_parameters, ) dataset = _build_map_dataset(sky_map, variables, calibration) @@ -189,61 +200,6 @@ def lo_l2( ] -def _prepare_corrections( - map_descriptor: MapDescriptor, - anc_dependencies: list, -) -> tuple[bool, Path | None, bool, bool, bool]: - """ - Determine which of the corrections the map descriptor asks for are needed. - - Parameters - ---------- - map_descriptor : MapDescriptor - The parsed map descriptor containing species and data type information. - anc_dependencies : list - List of ancillary file paths. - - Returns - ------- - tuple[bool, Path | None, bool, bool, bool] - A tuple containing: - - flux_correction: Whether to apply flux corrections - - flux_factors: Path to flux factors ancillary file if needed, - None otherwise - - sputter_correction: Whether to remove the counts sputtered into the - mapped species from a heavier one. - - bootstrap_correction: Whether to remove the intensity that bled into - each ESA level from the levels above it. - - cg_correction: Whether to apply CG correction to the dataset. - """ - # Default values - no corrections needed - flux_correction = False - flux_factors: None | Path = None - - if not map_descriptor.raw: - flux_correction = True - try: - flux_factors = next( - x for x in anc_dependencies if "esa-eta-fit-factors" in str(x) - ) - except StopIteration: - raise ValueError( - "No flux correction factor file found in ancillary dependencies" - ) from None - - sputter_correction = map_descriptor.sputter_corrected - bootstrap_correction = map_descriptor.bootstrap_corrected - cg_correction = True if map_descriptor.frame_descriptor == "hf" else False - - return ( - flux_correction, - flux_factors, - sputter_correction, - bootstrap_correction, - cg_correction, - ) - - # ============================================================================= # SETUP AND INITIALIZATION HELPERS # ============================================================================= @@ -396,6 +352,147 @@ def _bootstrap_correction() -> np.ndarray: return matrix +def _flux_corrector(anc_dependencies: list) -> PowerLawFluxCorrector: + """ + Load the ESA transmission factors the Compton-Getting correction reads. + + Parameters + ---------- + anc_dependencies : list + List of ancillary file paths, searched for the ESA eta fit factors. + + Returns + ------- + PowerLawFluxCorrector + The transmission factors, which the correction recovers the source + spectrum of a pixel through. + + Raises + ------ + ValueError + If the ancillary dependencies hold no ESA eta fit factors. + """ + try: + flux_factors = next( + x for x in anc_dependencies if "esa-eta-fit-factors" in str(x) + ) + except StopIteration: + raise ValueError( + "A heliospheric frame map needs the ESA eta fit factors to be " + "Compton-Getting corrected, and none were found in the ancillary " + "dependencies" + ) from None + + return PowerLawFluxCorrector(flux_factors) + + +def load_isn_mask_parameters() -> pd.DataFrame: + """ + Load the ISN mask tuning parameters shipped with the package. + + Returns + ------- + pd.DataFrame + The ancillary data, with columns: pivot_angle, esa_step, + intensity_threshold_fraction, angular_width_deg, outlier_percentile. + + Raises + ------ + ValueError + If no ISN mask parameter ancillary is shipped with the package. + """ + mask_files = sorted(ANCILLARY_DATA_DIR.glob("*isn-mask-parameters*")) + + if not mask_files: + raise ValueError("No ISN mask parameter files found") + + return lo_ancillary.read_ancillary_file(mask_files[-1]) + + +def _isn_mask_parameters(pivot_angle: int) -> pd.DataFrame: + """ + Get the ISN mask tuning of one pivot angle, in ascending ESA level order. + + Parameters + ---------- + pivot_angle : int + The nominal pivot angle [degrees] of the map being made. + + Returns + ------- + pd.DataFrame + The tuning of each ESA level, indexed by 1-based ESA step. + + Raises + ------ + ValueError + If the ancillary has no tuning for the pivot angle, which is not a + pivot the mask has been tuned for. + """ + parameters = load_isn_mask_parameters() + parameters = parameters[parameters["pivot_angle"] == pivot_angle] + + if parameters.empty: + raise ValueError( + f"The map asks for the ISN band to be masked out, but the ancillary " + f"has no mask tuning for the {pivot_angle} degree pivot angle" + ) + + # Select the ESA levels, in order. Raises if the ancillary is missing one. + return parameters.set_index("esa_step").loc[list(range(1, c.N_ESA_LEVELS + 1))] + + +def _isn_mask( + intensity: np.ndarray, elevation: np.ndarray, parameters: pd.DataFrame +) -> np.ndarray: + """ + Find the pixels the interstellar neutral flow dominates the map in. + + The ISN hydrogen the instrument sees is not the heliospheric ENA signal the + map is of, and it is bright enough to swamp it. It arrives as a band along + the ecliptic plane, so the mask is the pixels of a level that are both + bright and close enough to the plane, plus the brightest few pixels of the + level wherever they are. + + Parameters + ---------- + intensity : np.ndarray + The uncorrected intensity of every ESA level, of shape + (epoch, esa level, pixel). + elevation : np.ndarray + The ecliptic latitude [degrees] of each pixel, of shape (pixel,). The + ISN band lies along zero elevation. + parameters : pd.DataFrame + The mask tuning of each ESA level, in ascending level order. + + Returns + ------- + np.ndarray + Whether each pixel of each level is masked, of the shape of + ``intensity``. + """ + threshold = parameters["intensity_threshold_fraction"].to_numpy()[:, np.newaxis] + angular_width = parameters["angular_width_deg"].to_numpy()[:, np.newaxis] + percentile = parameters["outlier_percentile"].to_numpy() + + # The band: the bright pixels of a level, near enough to the ecliptic. A + # level that saw nothing has no brightest pixel to take a fraction of. + peak = np.nanmax(intensity, axis=-1, keepdims=True) + bright = np.where(peak > 0, intensity >= threshold * peak, False) + mask = bright & (np.abs(elevation) <= angular_width) + + # The outliers: the top tail of a level's own intensity distribution, which + # catches the ISN pixels that sit off the plane. + cutoff = np.stack( + [ + np.nanpercentile(intensity[:, level], level_percentile, axis=-1) + for level, level_percentile in enumerate(percentile) + ], + axis=1, + ) + return mask | (intensity > cutoff[..., np.newaxis]) + + def finalize_dataset(dataset: xr.Dataset, descriptor: str) -> xr.Dataset: """ Add attributes and perform final dataset preparation. @@ -629,7 +726,6 @@ def _accumulate_pointing( map_descriptor: MapDescriptor, energy: np.ndarray, sputter_source: str | None = None, - accumulate_cos_alpha: bool = False, ) -> None: """ Add one pointing's counts and exposure to the map. @@ -653,10 +749,6 @@ def _accumulate_pointing( sputter_source : str | None The species sputtering into the mapped species, whose counts are accumulated alongside. None if this map is not sputter corrected. - accumulate_cos_alpha : bool - Whether to accumulate the RAM projection of each spin-angle bin, which - the Compton-Getting correction needs. False if this map is not CG - corrected. """ species = map_descriptor.species pivot_angle = float(np.atleast_1d(goodtimes["pivot"].values)[0]) @@ -697,7 +789,7 @@ def _accumulate_pointing( values["sputter_source_count"] = ( histrates[f"{sputter_source}_counts"].values[in_goodtime].sum(axis=0) ) - if accumulate_cos_alpha: + if map_descriptor.cg_corrected: # The RAM projection is a property of the bin, not of the counts in it, # so like the background it is accumulated weighted by exposure and # divided by the total exposure at the end. @@ -1399,6 +1491,7 @@ def _calculate_rates_and_intensities( sputter_matrix: np.ndarray | None = None, bootstrap_matrix: np.ndarray | None = None, flux_corrector: PowerLawFluxCorrector | None = None, + isn_mask_parameters: pd.DataFrame | None = None, ) -> dict[str, np.ndarray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -1424,6 +1517,10 @@ def _calculate_rates_and_intensities( The ESA transmission factors the Compton-Getting correction recovers the source spectrum with. None leaves the intensities in the spacecraft frame. + isn_mask_parameters : pd.DataFrame | None + The tuning of the ISN mask, which blanks out the pixels the + interstellar neutral flow dominates once every correction has been + made. None leaves the whole map in place. Returns ------- @@ -1481,6 +1578,18 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: intensity = _divide(count_rate, geometric_factor * energy) intensity_stat_uncert = _divide(count_rate_stat_uncert, geometric_factor * energy) + # The ISN mask reads the intensity as the instrument observed it, with none + # of the corrections below made, so that a map is masked the same way + # regardless of the map descriptor. + isn_mask = None + if isn_mask_parameters is not None: + logger.info("Masking the ISN band out of the map") + isn_mask = _isn_mask( + _divide(_divide(counts, exposure), geometric_factor * energy), + sky_map.az_el_points.values[:, 1], + isn_mask_parameters, + ) + # The systematic error is the intensity excursion from the recalibrated # G-factor bounds, and the symmetric error is the geometric mean of the two. # Intensity goes as 1/G, so the lower G-factor bound gives the upper @@ -1550,7 +1659,7 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: flux_corrector, ) - return { + variables = { "ena_count": counts, "exposure_factor": exposure, "ena_count_rate": count_rate, @@ -1568,6 +1677,13 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: "bg_intensity_stat_uncert": bg_intensity_stat_uncert, } + # A masked pixel has no ENA measurement to report. + if isn_mask is not None: + for name in ISN_MASKED_VARIABLES: + variables[name] = np.where(isn_mask, np.nan, variables[name]) + + return variables + def _build_map_dataset( sky_map: RectangularSkyMap, diff --git a/imap_processing/tests/ena_maps/test_naming.py b/imap_processing/tests/ena_maps/test_naming.py index 195c23876..fd797c99d 100644 --- a/imap_processing/tests/ena_maps/test_naming.py +++ b/imap_processing/tests/ena_maps/test_naming.py @@ -570,3 +570,49 @@ def test_bootstrap_corrected(self, principal_data, expected): f"l090-{principal_data}-h-sf-nsp-ram-hae-6deg-1yr" ) assert md.bootstrap_corrected is expected + + @pytest.mark.parametrize( + "frame, expected", + [ + # The heliospheric frame is the one the correction moves a map into + ("hf", True), + # A map left in the frame it was observed in is not corrected + ("sf", False), + ("hk", False), + ], + ) + def test_cg_corrected(self, frame, expected): + md = MapDescriptor.from_string(f"l090-enasbs-h-{frame}-nsp-ram-hae-6deg-1yr") + assert md.cg_corrected is expected + + @pytest.mark.parametrize( + "principal_data, expected", + [ + # "msk" follows the "bs" or "nbs" of the bootstrap correction, + # whichever of the two the map asks for + ("enasbsmsk", True), + ("enasnbsmsk", True), + ("enansbsmsk", True), + ("enansnbsmsk", True), + ("enabsmsk", True), + # The SOC writes the code capitalized; both spellings are accepted + ("enasbsMsk", True), + # Without the code the ISN band is left in the map + ("enasbs", False), + ("enansnbs", False), + # The code has to follow the bootstrap one to be a mask code + ("enamsk", False), + ("enamsksbs", False), + # A raw map asks for none of the corrections + ("enaraw", False), + # Only ENA maps are ISN masked + ("spx0305", False), + ("isn", False), + ("drt", False), + ], + ) + def test_isn_masked(self, principal_data, expected): + md = MapDescriptor.from_string( + f"l090-{principal_data}-h-sf-nsp-ram-hae-6deg-1yr" + ) + assert md.isn_masked is expected diff --git a/imap_processing/tests/lo/test_anc/imap_lo_isn-mask-parameters-small_v001.csv b/imap_processing/tests/lo/test_anc/imap_lo_isn-mask-parameters-small_v001.csv new file mode 100644 index 000000000..80e575973 --- /dev/null +++ b/imap_processing/tests/lo/test_anc/imap_lo_isn-mask-parameters-small_v001.csv @@ -0,0 +1,22 @@ +pivot_angle,esa_step,intensity_threshold_fraction,angular_width_deg,outlier_percentile +75,1,0.5,1.0,100.0 +75,2,0.5,1.0,100.0 +75,3,0.5,1.0,100.0 +75,4,0.5,1.0,100.0 +75,5,0.5,1.0,100.0 +75,6,0.5,1.0,100.0 +75,7,0.5,1.0,100.0 +90,1,0.5,90.0,100.0 +90,2,0.5,90.0,100.0 +90,3,0.5,90.0,100.0 +90,4,0.5,90.0,100.0 +90,5,0.5,90.0,100.0 +90,6,0.5,90.0,100.0 +90,7,0.5,90.0,100.0 +105,1,2.0,90.0,50.0 +105,2,2.0,90.0,50.0 +105,3,2.0,90.0,50.0 +105,4,2.0,90.0,50.0 +105,5,2.0,90.0,50.0 +105,6,2.0,90.0,50.0 +105,7,2.0,90.0,50.0 diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 79ab69a36..195c129ed 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -16,6 +16,7 @@ ANCILLARY_DATA_DIR as PACKAGE_ANCILLARY_DIR, ) from imap_processing.lo.l2.lo_l2 import ( + ISN_MASKED_VARIABLES, LoSpinAnglePointingSet, _bootstrap_correct_intensity, _complete_pointings, @@ -28,6 +29,7 @@ finalize_dataset, lo_l2, load_bootstrap_correction_data, + load_isn_mask_parameters, load_sputter_correction_data, ) from imap_processing.spice.time import met_to_ttj2000ns @@ -50,6 +52,22 @@ # Compton-Getting correction. Neither of the other corrections is made. CG_DESCRIPTOR = "l090-enansnbs-h-hf-nsp-full-hae-6deg-3mo" +# The same uncorrected full-spin map, with the ISN band masked out. The "msk" +# follows the bootstrap code. The mask is tuned per pivot angle, which is the +# sensor of the descriptor, so each pivot gets its own descriptor: 90 masks the +# bright pixels, 75 has too narrow a band to mask any, and 105 masks the top +# half of each ESA level. See imap_lo_isn-mask-parameters-small_v001.csv. +MASK_DESCRIPTOR = "l090-enansnbsmsk-h-sf-nsp-full-hae-6deg-3mo" +NARROW_MASK_DESCRIPTOR = "l075-enansnbsmsk-h-sf-nsp-full-hae-6deg-3mo" +OUTLIER_MASK_DESCRIPTOR = "l105-enansnbsmsk-h-sf-nsp-full-hae-6deg-3mo" + +# A pivot angle the mask ancillary carries no tuning for. +UNTUNED_MASK_DESCRIPTOR = "l060-enansnbsmsk-h-sf-nsp-full-hae-6deg-3mo" + +# The fraction of an ESA level's peak intensity that the pivot 90 tuning masks +# from, and the variables it blanks out. +MASK_THRESHOLD_FRACTION = 0.5 + # The contents of imap_lo_sputter-correction-factors-small, as # {target ESA step: {source ESA step: factor}}, 1-based as in the ancillary. SPUTTER_FACTORS = {2: {3: 0.5}, 5: {3: 0.25, 6: 0.1}} @@ -220,8 +238,8 @@ def one_pointing(): def anc_dependencies(): """The ancillary files a map takes as a dependency. - Every map but a raw one is flux corrected, so the ESA eta fit factors are - required; without them the map cannot be made. + A heliospheric frame map is Compton-Getting corrected, which is the only + thing the ESA eta fit factors are read for. """ return [ANCILLARY_DIR / "imap_lo_esa-eta-fit-factors_20240101_v001.csv"] @@ -305,6 +323,24 @@ def cg_maps(one_pointing, anc_dependencies): return corrected, raw +@pytest.fixture +def masked_maps(one_pointing, anc_dependencies): + """The ISN masked and unmasked maps of the same pointing. + + The two differ only in whether the band was masked out, so the unmasked map + supplies the intensities the mask is predicted from. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (masked,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, MASK_DESCRIPTOR + ) + (raw,) = lo_l2(as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR) + return masked, raw + + class TestMapStructure: """The shape and contents of the produced map.""" @@ -831,6 +867,21 @@ def test_map_with_no_spectrum_falls_back_to_the_nominal_index(self): class TestComptonGettingCorrection: """Moving the intensities into the frame the heliosphere sees them in.""" + def test_the_eta_fit_factors_are_required(self, one_pointing): + """A heliospheric frame map cannot be made without the transmission factors.""" + with pytest.raises(ValueError, match="ESA eta fit factors"): + lo_l2(as_dependencies(one_pointing), [], CG_DESCRIPTOR) + + def test_a_spacecraft_frame_map_needs_no_eta_fit_factors(self, one_pointing): + """The factors are read for the correction alone, so an sf map goes without.""" + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (dataset,) = lo_l2(as_dependencies(one_pointing), [], FULL_DESCRIPTOR) + + assert (dataset["ena_intensity"].values > 0).any() + def test_counts_and_rates_stay_as_observed(self, cg_maps): """The correction applies to the intensities alone.""" corrected, raw = cg_maps @@ -1150,14 +1201,14 @@ def test_sputter_factors_are_selected_by_species_pair( self, shipped_ancillaries, tmp_path ): """Only the rows of the requested pair come back, in ESA step order.""" - factors = load_sputter_correction_data("o", "h") + factors = load_sputter_correction_data() assert list(factors.columns) == [ "source_species", "target_species", - "esa_step", + "target_esa", + "source_esa", "sputter_factor", - "sputter_factor_uncertainty", ] assert not factors.empty assert (factors["source_species"] == "o").all() @@ -1165,7 +1216,7 @@ def test_sputter_factors_are_selected_by_species_pair( with patch("imap_processing.lo.l2.lo_l2.ANCILLARY_DATA_DIR", tmp_path): with pytest.raises(ValueError, match="No sputter correction files"): - load_sputter_correction_data("o", "h") + load_sputter_correction_data() def test_bootstrap_factors_relate_lower_steps_to_higher( self, shipped_ancillaries, tmp_path @@ -1190,6 +1241,146 @@ def test_bootstrap_factors_relate_lower_steps_to_higher( load_bootstrap_correction_data() +class TestIsnMask: + """Blanking out the pixels the interstellar neutral flow dominates.""" + + def test_the_bright_pixels_of_every_level_are_masked(self, masked_maps): + """A pixel above the level's threshold comes back undefined.""" + masked, raw = masked_maps + + intensity = raw["ena_intensity"].values + peak = np.nanmax(intensity, axis=(-2, -1), keepdims=True) + expected = intensity >= MASK_THRESHOLD_FRACTION * peak + assert expected.any(), "the tuning must mask something" + assert not expected.all(), "the tuning must leave something unmasked" + + np.testing.assert_array_equal( + np.isnan(masked["ena_intensity"].values), expected + ) + + def test_the_unmasked_pixels_keep_their_values(self, masked_maps): + """Masking changes nothing about the pixels it does not blank out.""" + masked, raw = masked_maps + + keep = ~np.isnan(masked["ena_intensity"].values) + np.testing.assert_allclose( + masked["ena_intensity"].values[keep], + raw["ena_intensity"].values[keep], + rtol=1e-6, + ) + + def test_every_species_variable_is_masked_together(self, masked_maps): + """The counts, rates and uncertainties are blanked with the intensity.""" + masked, _ = masked_maps + + expected = np.isnan(masked["ena_intensity"].values) + for name in ISN_MASKED_VARIABLES: + np.testing.assert_array_equal( + np.isnan(masked[name].values), expected, err_msg=name + ) + + def test_the_exposure_and_background_are_left_alone(self, masked_maps): + """The mask says nothing about what the map was pointed at.""" + masked, raw = masked_maps + + for name in ( + "exposure_factor", + "bg_rate", + "bg_rate_stat_uncert", + "bg_intensity", + "bg_intensity_stat_uncert", + ): + np.testing.assert_allclose( + masked[name].values, raw[name].values, rtol=1e-6, err_msg=name + ) + + def test_an_unmasked_map_is_left_whole(self, full_map): + """A descriptor without the mask code blanks nothing out.""" + dataset, _ = full_map + + for name in ISN_MASKED_VARIABLES: + assert not np.isnan(dataset[name].values).any(), name + + def test_a_narrow_band_masks_nothing_off_the_ecliptic( + self, one_pointing, anc_dependencies + ): + """A pixel outside the angular width of the band survives its brightness.""" + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (masked,) = lo_l2( + as_dependencies(one_pointing), + anc_dependencies, + NARROW_MASK_DESCRIPTOR, + ) + + # The mocked sky pointing puts every pixel of the test map on the + # ecliptic, which the pivot 75 tuning is too narrow to reach. + assert not np.isnan(masked["ena_intensity"].values).any() + + def test_the_outlier_tail_of_a_level_is_masked( + self, one_pointing, anc_dependencies + ): + """The pixels above the level's percentile are masked wherever they are.""" + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (masked,) = lo_l2( + as_dependencies(one_pointing), + anc_dependencies, + OUTLIER_MASK_DESCRIPTOR, + ) + (raw,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, FULL_DESCRIPTOR + ) + + # The pivot 105 tuning cannot mask a band, so only the top half of each + # level's own intensity distribution is masked. + intensity = raw["ena_intensity"].values + median = np.nanpercentile(intensity, 50, axis=(-2, -1), keepdims=True) + expected = intensity > median + assert expected.any(), "the tuning must mask something" + + np.testing.assert_array_equal( + np.isnan(masked["ena_intensity"].values), expected + ) + + def test_an_untuned_pivot_angle_is_refused(self, one_pointing, anc_dependencies): + """A map cannot be masked at a pivot the ancillary says nothing about.""" + with pytest.raises(ValueError, match="no mask tuning for the 60 degree"): + lo_l2( + as_dependencies(one_pointing), + anc_dependencies, + UNTUNED_MASK_DESCRIPTOR, + ) + + def test_mask_parameters_cover_every_level_of_every_tuned_pivot( + self, shipped_ancillaries, tmp_path + ): + """The shipped ancillary tunes all 7 levels of each pivot it carries.""" + parameters = load_isn_mask_parameters() + + assert list(parameters.columns) == [ + "pivot_angle", + "esa_step", + "intensity_threshold_fraction", + "angular_width_deg", + "outlier_percentile", + ] + assert not parameters.empty + for pivot, tuning in parameters.groupby("pivot_angle"): + assert sorted(tuning["esa_step"]) == list(range(1, N_ESA + 1)), pivot + assert (parameters["intensity_threshold_fraction"] > 0).all() + assert (parameters["angular_width_deg"] > 0).all() + assert parameters["outlier_percentile"].between(0, 100).all() + + with patch("imap_processing.lo.l2.lo_l2.ANCILLARY_DATA_DIR", tmp_path): + with pytest.raises(ValueError, match="No ISN mask parameter files"): + load_isn_mask_parameters() + + class TestFinalizeDataset: """Attaching the CDF attributes to a finished map.""" From 6fefb194767104d94a7c8a8279bd7f9929389c23 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Tue, 4 Aug 2026 17:01:25 -0400 Subject: [PATCH 11/14] combined maps --- imap_processing/lo/l2/lo_l2.py | 142 +++++++++-- imap_processing/tests/ena_maps/test_naming.py | 10 + imap_processing/tests/lo/test_lo_l2.py | 221 +++++++++++++++++- imap_processing/tests/test_cli.py | 48 ++++ 4 files changed, 401 insertions(+), 20 deletions(-) diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 2c3a6d94b..4bac0477f 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -78,7 +78,10 @@ def lo_l2( The inputs are expected to have already been filtered down to the pivot angle of the map being made, which is done in pre-processing (see ``cli.Lo.pre_processing``) so that the map records only the files it was - made from as its parents. + made from as its parents. A combined map, written "ilo" rather than with a + pivot angle of its own, is filtered by nothing and accumulates every + pointing it is given. Each pointing is projected from the pivot angle its + own goodtimes report. Parameters ---------- @@ -113,17 +116,12 @@ def lo_l2( map_descriptor = MapDescriptor.from_string(descriptor) logger.info(f"Processing map for species: {map_descriptor.species}") - # Read the ancillaries the corrections asked for by the descriptor need up - # front, so that a map missing one fails before any of its pointings are - # accumulated. + # The Compton-Getting correction reads the source spectrum of each pixel + # through the ESA transmission factors of the eta fit ancillary. Read up + # front, so that a map missing it fails before anything is accumulated. flux_corrector = ( _flux_corrector(anc_dependencies) if map_descriptor.cg_corrected else None ) - isn_mask_parameters = ( - _isn_mask_parameters(int(map_descriptor.sensor)) - if map_descriptor.isn_masked - else None - ) logger.info("Step 1: Loading ancillary data") _efficiency_data = load_efficiency_data(anc_dependencies) @@ -142,6 +140,15 @@ def lo_l2( pointings = _complete_pointings(sci_dependencies) logger.info(f"Building {descriptor} from {len(pointings)} pointings") + # The mask is tuned per pivot angle, which a combined map takes from the + # pointings themselves. Resolved before anything is accumulated, so that a + # map the mask cannot be tuned for fails before the work is done. + isn_mask_parameters = ( + _isn_mask_parameters(_map_pivot_angles(map_descriptor, pointings)) + if map_descriptor.isn_masked + else None + ) + # Every pointing of a map is taken in the same ESA mode, so the last one # sets the energy response the whole map is binned in. esa_mode = _get_esa_mode(pointings[max(pointings)][2]) if pointings else 0 @@ -409,14 +416,20 @@ def load_isn_mask_parameters() -> pd.DataFrame: return lo_ancillary.read_ancillary_file(mask_files[-1]) -def _isn_mask_parameters(pivot_angle: int) -> pd.DataFrame: +def _isn_mask_parameters(pivot_angles: list[int]) -> pd.DataFrame: """ - Get the ISN mask tuning of one pivot angle, in ascending ESA level order. + Get the ISN mask tuning of a map, in ascending ESA level order. + + A map made at one pivot angle is masked with that pivot's tuning. A map + combining several is masked with the most permissive tuning of the pivots + that went into it: a pixel of the combined map holds the interstellar + neutrals seen at every one of them, so it is masked if any of those pivots + would have masked it. Parameters ---------- - pivot_angle : int - The nominal pivot angle [degrees] of the map being made. + pivot_angles : list[int] + The nominal pivot angles [degrees] the map was built from. Returns ------- @@ -426,20 +439,30 @@ def _isn_mask_parameters(pivot_angle: int) -> pd.DataFrame: Raises ------ ValueError - If the ancillary has no tuning for the pivot angle, which is not a - pivot the mask has been tuned for. + If the ancillary has no tuning for one of the pivot angles. """ parameters = load_isn_mask_parameters() - parameters = parameters[parameters["pivot_angle"] == pivot_angle] - if parameters.empty: + untuned = sorted(set(pivot_angles) - set(parameters["pivot_angle"])) + if untuned: raise ValueError( f"The map asks for the ISN band to be masked out, but the ancillary " - f"has no mask tuning for the {pivot_angle} degree pivot angle" + f"has no mask tuning for the {untuned} degree pivot angle(s) it was " + f"built from" ) + parameters = parameters[parameters["pivot_angle"].isin(pivot_angles)] + + # The widest band, the faintest pixel taken as bright, and the shortest + # outlier tail, i.e. the union of what each contributing pivot would mask. + tuning = parameters.groupby("esa_step").agg( + intensity_threshold_fraction=("intensity_threshold_fraction", "min"), + angular_width_deg=("angular_width_deg", "max"), + outlier_percentile=("outlier_percentile", "min"), + ) + # Select the ESA levels, in order. Raises if the ancillary is missing one. - return parameters.set_index("esa_step").loc[list(range(1, c.N_ESA_LEVELS + 1))] + return tuning.loc[list(range(1, c.N_ESA_LEVELS + 1))] def _isn_mask( @@ -542,6 +565,87 @@ def finalize_dataset(dataset: xr.Dataset, descriptor: str) -> xr.Dataset: # ============================================================================= +def _nominal_pivot_angle(pivot_angle: float) -> int | None: + """ + Snap a measured pivot angle onto the nominal pivot angle it was flown at. + + Parameters + ---------- + pivot_angle : float + The pivot angle [degrees] a pointing's goodtimes report. + + Returns + ------- + int | None + The nominal pivot angle [degrees] whose range contains it, or None if + it falls in none of them. + """ + for nominal, spec in c.PIVOT_ANGLES.items(): + if spec.min <= pivot_angle <= spec.max: + return nominal + return None + + +def _map_pivot_angles( + map_descriptor: MapDescriptor, pointings: dict[int, tuple] +) -> list[int]: + """ + Get the nominal pivot angles a map is built from. + + A Lo map carries its pivot angle as its sensor, e.g. the 90 of "l090", and + its inputs were filtered down to that pivot in pre-processing. A map that + combines every pivot angle instead of selecting one is written without a + sensor, as "ilo", so the pivot angles it holds are the ones its pointings + were actually flown at, which their goodtimes report. + + Parameters + ---------- + map_descriptor : MapDescriptor + The parsed descriptor of the map being made. + pointings : dict[int, tuple] + The (goodtimes, bgrates, histrates) datasets of each mappable pointing. + + Returns + ------- + list[int] + The nominal pivot angles [degrees] of the map, in ascending order. + + Raises + ------ + ValueError + If the map combines pivot angles but none of its pointings reports one + that is recognisably nominal, leaving nothing to identify it by. + """ + if isinstance(map_descriptor.sensor, int): + return [map_descriptor.sensor] + + measured = { + float(np.atleast_1d(goodtimes["pivot"].values)[0]) + for goodtimes, _, _ in pointings.values() + } + nominal = {_nominal_pivot_angle(pivot) for pivot in measured} + + unrecognised = sorted( + pivot for pivot in measured if _nominal_pivot_angle(pivot) is None + ) + if unrecognised: + logger.warning( + f"Ignoring the pivot angles {unrecognised} of " + f"{map_descriptor.instrument_descriptor}, they match none of the " + f"nominal pivot angles." + ) + + pivot_angles = sorted(pivot for pivot in nominal if pivot is not None) + if not pivot_angles: + raise ValueError( + f"The map asks for the ISN band to be masked out, but none of the " + f"pointings of {map_descriptor.instrument_descriptor} reports a " + f"nominal pivot angle to look the mask tuning up by" + ) + + return pivot_angles + + def _complete_pointings( sci_dependencies: dict[int, dict[str, xr.Dataset]], ) -> dict[int, tuple]: diff --git a/imap_processing/tests/ena_maps/test_naming.py b/imap_processing/tests/ena_maps/test_naming.py index fd797c99d..dcd92205a 100644 --- a/imap_processing/tests/ena_maps/test_naming.py +++ b/imap_processing/tests/ena_maps/test_naming.py @@ -80,6 +80,16 @@ def test_parse_instrument_descriptor( MappableInstrumentShortName.ULTRA, "combined", ) + # A Lo map of no particular pivot angle carries no sensor, which is + # what tells it apart from the "l090" of a single pivot angle + assert MapDescriptor.parse_instrument_descriptor("ilo") == ( + MappableInstrumentShortName.LO, + "", + ) + assert MapDescriptor.parse_instrument_descriptor("l090") == ( + MappableInstrumentShortName.LO_HI_RES, + 90, + ) with pytest.raises( ValueError, match="'abc' is not a valid MappableInstrumentShortName" ): diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 195c129ed..5a10779d1 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -64,6 +64,15 @@ # A pivot angle the mask ancillary carries no tuning for. UNTUNED_MASK_DESCRIPTOR = "l060-enansnbsmsk-h-sf-nsp-full-hae-6deg-3mo" +# The combined maps, written "ilo" rather than with a pivot angle of their own. +# These select no pivot angle, so every pointing given to them is accumulated. +COMBINED_DESCRIPTOR = "ilo-enansnbs-h-sf-nsp-full-hae-6deg-3mo" +COMBINED_RAM_DESCRIPTOR = "ilo-enansnbs-h-sf-nsp-ram-hae-6deg-3mo" +COMBINED_MASK_DESCRIPTOR = "ilo-enansnbsmsk-h-sf-nsp-full-hae-6deg-3mo" + +# The pivot angles a combined map is built from in these tests. +COMBINED_PIVOTS = (75.0, 90.0, 105.0) + # The fraction of an ESA level's peak intensity that the pivot 90 tuning masks # from, and the variables it blanks out. MASK_THRESHOLD_FRACTION = 0.5 @@ -341,6 +350,216 @@ def masked_maps(one_pointing, anc_dependencies): return masked, raw +@pytest.fixture +def pointings_at_every_pivot(): + """One synthetic pointing at each of the pivot angles a map combines.""" + return [ + make_pointing(repointing=100 + i, pivot=pivot, seed=i + 1) + for i, pivot in enumerate(COMBINED_PIVOTS) + ] + + +class TestCombinedMap: + """A map that selects no pivot angle and takes every pointing it is given.""" + + def test_every_pivot_angle_is_accumulated( + self, pointings_at_every_pivot, anc_dependencies + ): + """The counts of all three pointings land on the one combined map.""" + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (combined,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + COMBINED_DESCRIPTOR, + ) + + expected = sum( + pointing["expected_counts"].sum() for pointing in pointings_at_every_pivot + ) + assert combined["ena_count"].values.sum() == pytest.approx(expected) + + def test_the_combined_map_is_the_sum_of_its_pointings( + self, pointings_at_every_pivot, anc_dependencies + ): + """Accumulating together matches accumulating each pointing alone.""" + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (combined,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + COMBINED_DESCRIPTOR, + ) + singly = [ + lo_l2(as_dependencies(pointing), anc_dependencies, COMBINED_DESCRIPTOR)[ + 0 + ] + for pointing in pointings_at_every_pivot + ] + + for variable in ("ena_count", "exposure_factor"): + np.testing.assert_allclose( + combined[variable].values, + sum(single[variable].values for single in singly), + rtol=1e-5, + err_msg=variable, + ) + + def test_a_pointing_is_projected_from_its_own_pivot_angle( + self, pointings_at_every_pivot, anc_dependencies + ): + """The pivot angles of a combined map need not agree with each other. + + A RAM map keeps the spin-angle bins looking into the RAM direction, + which is a function of the pointing's own pivot angle, so the three + pointings would not all survive one shared pivot angle. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (combined,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + COMBINED_RAM_DESCRIPTOR, + ) + singly = [ + lo_l2( + as_dependencies(pointing), + anc_dependencies, + COMBINED_RAM_DESCRIPTOR, + )[0] + for pointing in pointings_at_every_pivot + ] + + # Each pointing keeps its own RAM half of the spin, and the combined + # map is exactly those three halves added together. + assert (combined["ena_count"].values > 0).any() + np.testing.assert_allclose( + combined["ena_count"].values, + sum(single["ena_count"].values for single in singly), + rtol=1e-5, + ) + + def test_a_pivot_angle_map_of_the_same_pointings_is_no_larger( + self, pointings_at_every_pivot, anc_dependencies + ): + """The combined map holds at least what any single-pivot map holds. + + ``lo_l2`` itself does no pivot filtering, so a single-pivot descriptor + given these same pointings accumulates them all too; the point here is + that the combined descriptor loses nothing by naming no pivot angle. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (combined,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + COMBINED_DESCRIPTOR, + ) + (single_pivot,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + FULL_DESCRIPTOR, + ) + + np.testing.assert_allclose( + combined["ena_count"].values, single_pivot["ena_count"].values, rtol=1e-5 + ) + + def test_the_isn_mask_reads_the_pivot_angles_off_the_pointings( + self, pointings_at_every_pivot, anc_dependencies + ): + """A combined map is masked with the union of its pivots' tunings. + + The pivot 105 tuning of the test ancillary masks the top half of each + ESA level and the pivot 90 one the bright pixels, so a map holding both + is masked at least everywhere either of them would mask. + """ + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + (combined,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + COMBINED_MASK_DESCRIPTOR, + ) + (unmasked,) = lo_l2( + as_dependencies(*pointings_at_every_pivot), + anc_dependencies, + COMBINED_DESCRIPTOR, + ) + + masked = np.isnan(combined["ena_intensity"].values) + intensity = unmasked["ena_intensity"].values + + # The union of the three tunings: the widest band (pivot 75's 1 degree + # loses to the 90 degrees of the others), the faintest brightness taken + # as bright, and the shortest outlier tail. + peak = np.nanmax(intensity, axis=(-2, -1), keepdims=True) + median = np.nanpercentile(intensity, 50, axis=(-2, -1), keepdims=True) + expected = (intensity >= MASK_THRESHOLD_FRACTION * peak) | (intensity > median) + + assert expected.any(), "the tuning must mask something" + np.testing.assert_array_equal(masked, expected) + + def test_a_combined_map_of_one_pivot_angle_uses_that_pivots_tuning( + self, one_pointing, anc_dependencies + ): + """With a single pivot angle in the data the union is that pivot's row.""" + with patch( + "imap_processing.lo.l1c.lo_l1c.frame_transform_az_el", + side_effect=identity_pointing, + ): + # The pointing is at pivot 90, so this must match the map made with + # the pivot 90 descriptor. + (combined,) = lo_l2( + as_dependencies(one_pointing), + anc_dependencies, + COMBINED_MASK_DESCRIPTOR, + ) + (by_descriptor,) = lo_l2( + as_dependencies(one_pointing), anc_dependencies, MASK_DESCRIPTOR + ) + + assert np.isnan(combined["ena_intensity"].values).any() + np.testing.assert_array_equal( + np.isnan(combined["ena_intensity"].values), + np.isnan(by_descriptor["ena_intensity"].values), + ) + + def test_an_unrecognisable_pivot_angle_cannot_tune_the_mask( + self, anc_dependencies, caplog + ): + """A pivot angle matching no nominal one is dropped, and none is left.""" + off_nominal = make_pointing(repointing=100, pivot=42.0) + + with pytest.raises(ValueError, match="reports a nominal pivot angle"): + lo_l2( + as_dependencies(off_nominal), + anc_dependencies, + COMBINED_MASK_DESCRIPTOR, + ) + assert "match none of the nominal pivot angles" in caplog.text + + def test_a_pivot_angle_the_mask_has_no_tuning_for_is_refused( + self, anc_dependencies + ): + """A combined map holding an untuned pivot angle cannot be masked.""" + # 60 is a nominal pivot angle, but the mask ancillary has no row for it. + untuned = make_pointing(repointing=100, pivot=60.0) + + with pytest.raises(ValueError, match=r"no mask tuning for the \[60\]"): + lo_l2(as_dependencies(untuned), anc_dependencies, COMBINED_MASK_DESCRIPTOR) + + class TestMapStructure: """The shape and contents of the produced map.""" @@ -1349,7 +1568,7 @@ def test_the_outlier_tail_of_a_level_is_masked( def test_an_untuned_pivot_angle_is_refused(self, one_pointing, anc_dependencies): """A map cannot be masked at a pivot the ancillary says nothing about.""" - with pytest.raises(ValueError, match="no mask tuning for the 60 degree"): + with pytest.raises(ValueError, match=r"no mask tuning for the \[60\] degree"): lo_l2( as_dependencies(one_pointing), anc_dependencies, diff --git a/imap_processing/tests/test_cli.py b/imap_processing/tests/test_cli.py index 38b143c95..99640ca90 100644 --- a/imap_processing/tests/test_cli.py +++ b/imap_processing/tests/test_cli.py @@ -581,6 +581,54 @@ def test_lo_pre_processing_pivot_angle_filter(mock_super_pre_processing, mock_lo assert mock_load_cdf.call_count == 2 +@mock.patch("imap_processing.cli.load_cdf") +@mock.patch("imap_processing.cli.ProcessInstrument.pre_processing") +def test_lo_pre_processing_combined_map_keeps_every_pivot_angle( + mock_super_pre_processing, mock_load_cdf +): + """Test that a combined map takes the pointings of every pivot angle.""" + first = "-repoint00217_v001.cdf" + second = "-repoint00218_v001.cdf" + goodtimes = [ + f"imap_lo_l1b_goodtimes_20250415{first}", + f"imap_lo_l1b_goodtimes_20250416{second}", + ] + histrates = [ + f"imap_lo_l1b_histrates_20250415{first}", + f"imap_lo_l1b_histrates_20250416{second}", + ] + bgrates = [f"imap_lo_l1b_bgrates_20250415{first}"] + ancillary = "imap_lo_efficiency-factors_20250415_v001.csv" + + base_collection = ProcessingInputCollection( + ScienceInput(*goodtimes), + ScienceInput(*histrates), + ScienceInput(*bgrates), + AncillaryInput(ancillary), + ) + mock_super_pre_processing.return_value = base_collection + + instrument = Lo( + "l2", + # "ilo" rather than "l090": a map of no particular pivot angle + "ilo-ena-h-sf-nsp-ram-hae-6deg-3mo", + base_collection.serialize(), + "20250415", + "20250715", + "v001", + False, + ) + result = instrument.pre_processing() + + # The two pointings are at different pivot angles, and both are kept. + assert [ + [str(file_path.filename) for file_path in processing_input.imap_file_paths] + for processing_input in result.get_processing_inputs() + ] == [goodtimes, histrates, bgrates, [ancillary]] + # No goodtimes are read, there being no pivot angle to select them by + assert mock_load_cdf.call_count == 0 + + @mock.patch("imap_processing.cli.load_cdf") @mock.patch("imap_processing.cli.ProcessInstrument.pre_processing") def test_lo_pre_processing_drops_goodtimes_without_pivot( From e63f74b849a191bf2f1706c9b5af66c42fa46121 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 5 Aug 2026 16:06:57 -0400 Subject: [PATCH 12/14] using fillval for pixels with no exposure --- imap_processing/lo/l2/lo_l2.py | 27 ++++- imap_processing/tests/lo/test_lo_l2.py | 140 ++++++++++++++++--------- 2 files changed, 115 insertions(+), 52 deletions(-) diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index 4bac0477f..c5f29d2ba 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -53,6 +53,25 @@ "ena_intensity_sys_err_minus", ) +# The fill value of every floating point map variable for the L2 map. +# A pixel holding this is one the map has no measurement for. +FILLVAL_FLOAT = -1.0e31 + +# The map variables that are filled with FILLVAL_FLOAT where the map was never exposed. +FILLED_VARIABLES = ( + "ena_count_rate", + "ena_count_rate_stat_uncert", + "ena_intensity", + "ena_intensity_stat_uncert", + "ena_intensity_sys_err", + "ena_intensity_sys_err_plus", + "ena_intensity_sys_err_minus", + "bg_rate", + "bg_rate_stat_uncert", + "bg_intensity", + "bg_intensity_stat_uncert", +) + # The calibration ancillaries shipped with the package. ANCILLARY_DATA_DIR = Path(__file__).parent.parent / "ancillary_data" @@ -1600,7 +1619,8 @@ def _calculate_rates_and_intensities( """ Turn the accumulated counts and exposure into rates and intensities. - Every quantity is zero in the pixels that were never exposed. + Every derived quantity is filled in with a fill value in the pixels + that were never exposed. Parameters ---------- @@ -1784,7 +1804,10 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: # A masked pixel has no ENA measurement to report. if isn_mask is not None: for name in ISN_MASKED_VARIABLES: - variables[name] = np.where(isn_mask, np.nan, variables[name]) + variables[name] = np.where(isn_mask, FILLVAL_FLOAT, variables[name]) + + for name in FILLED_VARIABLES: + variables[name] = np.where(exposed, variables[name], FILLVAL_FLOAT) return variables diff --git a/imap_processing/tests/lo/test_lo_l2.py b/imap_processing/tests/lo/test_lo_l2.py index 5a10779d1..46ef837ae 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -16,6 +16,8 @@ ANCILLARY_DATA_DIR as PACKAGE_ANCILLARY_DIR, ) from imap_processing.lo.l2.lo_l2 import ( + FILLED_VARIABLES, + FILLVAL_FLOAT, ISN_MASKED_VARIABLES, LoSpinAnglePointingSet, _bootstrap_correct_intensity, @@ -112,6 +114,20 @@ # The ESA passband half-widths of the same file, a tenth of each center energy. ESA_ENERGY_DELTAS = {mode: energies / 10 for mode, energies in ESA_ENERGIES.items()} +# The fill value as it lands in the map, which writes its variables as float32. +FILL = np.float32(FILLVAL_FLOAT) + + +def is_fill(values): + """Whether each element of a map variable holds the fill value.""" + return np.asarray(values) == FILL + + +def measured(values): + """The elements of a map variable that are not filled.""" + values = np.asarray(values) + return values[~is_fill(values)] + @pytest.fixture(autouse=True) def use_test_geometric_factors(): @@ -497,15 +513,21 @@ def test_the_isn_mask_reads_the_pivot_angles_off_the_pointings( COMBINED_DESCRIPTOR, ) - masked = np.isnan(combined["ena_intensity"].values) + masked = is_fill(combined["ena_intensity"].values) intensity = unmasked["ena_intensity"].values + unexposed = is_fill(intensity) # The union of the three tunings: the widest band (pivot 75's 1 degree # loses to the 90 degrees of the others), the faintest brightness taken # as bright, and the shortest outlier tail. - peak = np.nanmax(intensity, axis=(-2, -1), keepdims=True) - median = np.nanpercentile(intensity, 50, axis=(-2, -1), keepdims=True) - expected = (intensity >= MASK_THRESHOLD_FRACTION * peak) | (intensity > median) + as_masked = np.where(unexposed, 0.0, intensity) + peak = np.max(as_masked, axis=(-2, -1), keepdims=True) + median = np.percentile(as_masked, 50, axis=(-2, -1), keepdims=True) + expected = ( + (as_masked >= MASK_THRESHOLD_FRACTION * peak) + | (as_masked > median) + | unexposed + ) assert expected.any(), "the tuning must mask something" np.testing.assert_array_equal(masked, expected) @@ -529,10 +551,10 @@ def test_a_combined_map_of_one_pivot_angle_uses_that_pivots_tuning( as_dependencies(one_pointing), anc_dependencies, MASK_DESCRIPTOR ) - assert np.isnan(combined["ena_intensity"].values).any() + assert is_fill(combined["ena_intensity"].values).any() np.testing.assert_array_equal( - np.isnan(combined["ena_intensity"].values), - np.isnan(by_descriptor["ena_intensity"].values), + is_fill(combined["ena_intensity"].values), + is_fill(by_descriptor["ena_intensity"].values), ) def test_an_unrecognisable_pivot_angle_cannot_tune_the_mask( @@ -710,7 +732,7 @@ def test_rate_and_intensity(self, full_map): counts[exposed] / exposure[exposed], rtol=1e-5, ) - assert np.all(dataset["ena_count_rate"].values[~exposed] == 0) + assert np.all(is_fill(dataset["ena_count_rate"].values[~exposed])) geometric_factor = GEO_FACTORS[0] energy = ESA_ENERGIES[0] @@ -748,7 +770,7 @@ def test_background_rate(self, full_map): np.testing.assert_allclose( bg_rate[exposed], pointing["background"][energy_index], rtol=1e-5 ) - assert np.all(bg_rate[~exposed] == 0) + assert np.all(is_fill(bg_rate[~exposed])) def test_systematic_error_bounds(self, full_map): """The systematic error is bracketed by the G-factor excursions.""" @@ -863,8 +885,8 @@ def test_rate_is_never_negative(self, sputter_maps): ] assert np.any(over_subtracted), "no pixel exercised the floor" - assert np.all(corrected["ena_count_rate"].values >= 0) - assert np.all(corrected["ena_intensity"].values >= 0) + assert np.all(measured(corrected["ena_count_rate"].values) >= 0) + assert np.all(measured(corrected["ena_intensity"].values) >= 0) def test_counts_stay_as_observed(self, sputter_maps): """The correction applies from the rate onward, not to the raw counts.""" @@ -903,16 +925,19 @@ def test_correction_removes_the_bled_intensity(self, bootstrap_maps): """A corrected step loses the scaled intensity of the steps above it.""" corrected, raw = bootstrap_maps - intensity = raw["ena_intensity"].values + # The correction runs before the unexposed pixels are filled, so it saw + # the zero they were divided down to. + intensity = np.where( + is_fill(raw["ena_intensity"].values), 0.0, raw["ena_intensity"].values + ) for target_esa in self.MAPPED_SOURCE_STEPS: expected = np.maximum( intensity[:, target_esa - 1] - self.bled(intensity, target_esa), 0.0 ) - np.testing.assert_allclose( - corrected["ena_intensity"].values[:, target_esa - 1], - expected, - rtol=1e-4, - ) + step = corrected["ena_intensity"].values[:, target_esa - 1] + keep = ~is_fill(step) + assert keep.any() + np.testing.assert_allclose(step[keep], expected[keep], rtol=1e-4) def test_uncorrected_steps_are_untouched(self, bootstrap_maps): """A step that nothing bleeds into keeps the intensity it already had.""" @@ -948,22 +973,22 @@ def test_uncertainty_gains_the_source_intensities(self, bootstrap_maps): """Subtracting a measured quantity can only add to the variance.""" corrected, raw = bootstrap_maps - variance = raw["ena_intensity_stat_uncert"].values ** 2 + uncert = raw["ena_intensity_stat_uncert"].values + variance = np.where(is_fill(uncert), 0.0, uncert) ** 2 for target_esa in self.MAPPED_SOURCE_STEPS: expected = variance[:, target_esa - 1] + self.bled( variance, target_esa, power=2 ) - np.testing.assert_allclose( - corrected["ena_intensity_stat_uncert"].values[:, target_esa - 1], - np.sqrt(expected), - rtol=1e-4, - ) + step = corrected["ena_intensity_stat_uncert"].values[:, target_esa - 1] + keep = ~is_fill(step) + assert keep.any() + np.testing.assert_allclose(step[keep], np.sqrt(expected[keep]), rtol=1e-4) def test_intensity_is_never_negative(self, bootstrap_maps): """The corrected map holds no negative intensity.""" corrected, _ = bootstrap_maps - assert np.all(corrected["ena_intensity"].values >= 0) + assert np.all(measured(corrected["ena_intensity"].values) >= 0) def test_over_subtraction_is_floored_at_zero(self): """Over-subtracting a pixel floors it rather than going below zero.""" @@ -997,8 +1022,8 @@ def test_systematic_error_brackets_the_correction(self, bootstrap_maps): symmetric = corrected["ena_intensity_sys_err"].values lit = corrected["ena_intensity"].values > 0 - assert np.all(plus >= 0) - assert np.all(minus >= 0) + assert np.all(measured(plus) >= 0) + assert np.all(measured(minus) >= 0) assert np.all(plus[lit] > 0) np.testing.assert_allclose( symmetric[lit], np.sqrt(plus[lit] * minus[lit]), rtol=1e-4 @@ -1121,8 +1146,8 @@ def test_intensities_are_corrected(self, cg_maps): ~np.isclose(corrected[variable].values[lit], raw[variable].values[lit]) ), f"{variable} was not corrected" - def test_map_holds_no_fill_values(self, cg_maps): - """A pixel the correction says nothing about comes out at zero.""" + def test_the_correction_leaks_no_nan(self, cg_maps): + """A pixel the correction says nothing about comes out at zero, not NaN.""" corrected, _ = cg_maps for variable in ( @@ -1135,13 +1160,14 @@ def test_map_holds_no_fill_values(self, cg_maps): ): values = corrected[variable].values assert np.isfinite(values).all(), f"{variable} holds NaN or inf" - assert np.all(values >= 0), f"{variable} went negative" + assert np.all(measured(values) >= 0), f"{variable} went negative" def test_dark_pixels_stay_dark(self, cg_maps): """A pixel that saw nothing has no spectrum to correct.""" corrected, raw = cg_maps - dark = raw["ena_intensity"].values <= 0 + intensity = raw["ena_intensity"].values + dark = (intensity <= 0) & ~is_fill(intensity) assert dark.any() assert np.all(corrected["ena_intensity"].values[dark] == 0) @@ -1468,20 +1494,21 @@ def test_the_bright_pixels_of_every_level_are_masked(self, masked_maps): masked, raw = masked_maps intensity = raw["ena_intensity"].values - peak = np.nanmax(intensity, axis=(-2, -1), keepdims=True) - expected = intensity >= MASK_THRESHOLD_FRACTION * peak + # The unexposed pixels are filled in both maps, mask or no mask. + unexposed = is_fill(intensity) + as_masked = np.where(unexposed, 0.0, intensity) + peak = np.max(as_masked, axis=(-2, -1), keepdims=True) + expected = (as_masked >= MASK_THRESHOLD_FRACTION * peak) | unexposed assert expected.any(), "the tuning must mask something" assert not expected.all(), "the tuning must leave something unmasked" - np.testing.assert_array_equal( - np.isnan(masked["ena_intensity"].values), expected - ) + np.testing.assert_array_equal(is_fill(masked["ena_intensity"].values), expected) def test_the_unmasked_pixels_keep_their_values(self, masked_maps): """Masking changes nothing about the pixels it does not blank out.""" masked, raw = masked_maps - keep = ~np.isnan(masked["ena_intensity"].values) + keep = ~is_fill(masked["ena_intensity"].values) np.testing.assert_allclose( masked["ena_intensity"].values[keep], raw["ena_intensity"].values[keep], @@ -1490,13 +1517,17 @@ def test_the_unmasked_pixels_keep_their_values(self, masked_maps): def test_every_species_variable_is_masked_together(self, masked_maps): """The counts, rates and uncertainties are blanked with the intensity.""" - masked, _ = masked_maps + masked, raw = masked_maps - expected = np.isnan(masked["ena_intensity"].values) + # What the mask blanked, as opposed to what the unexposed pixels of the + # unmasked map are already filled with. + expected = is_fill(masked["ena_intensity"].values) & ~is_fill( + raw["ena_intensity"].values + ) + assert expected.any(), "the tuning must mask something" for name in ISN_MASKED_VARIABLES: - np.testing.assert_array_equal( - np.isnan(masked[name].values), expected, err_msg=name - ) + blanked = is_fill(masked[name].values) & ~is_fill(raw[name].values) + np.testing.assert_array_equal(blanked, expected, err_msg=name) def test_the_exposure_and_background_are_left_alone(self, masked_maps): """The mask says nothing about what the map was pointed at.""" @@ -1517,8 +1548,15 @@ def test_an_unmasked_map_is_left_whole(self, full_map): """A descriptor without the mask code blanks nothing out.""" dataset, _ = full_map + # Only the pixels the map was never exposed in are filled, and the raw + # accumulators keep even those. + unexposed = dataset["exposure_factor"].values == 0 for name in ISN_MASKED_VARIABLES: - assert not np.isnan(dataset[name].values).any(), name + filled = name in FILLED_VARIABLES + expected = unexposed if filled else np.zeros_like(unexposed) + np.testing.assert_array_equal( + is_fill(dataset[name].values), expected, err_msg=name + ) def test_a_narrow_band_masks_nothing_off_the_ecliptic( self, one_pointing, anc_dependencies @@ -1535,8 +1573,9 @@ def test_a_narrow_band_masks_nothing_off_the_ecliptic( ) # The mocked sky pointing puts every pixel of the test map on the - # ecliptic, which the pivot 75 tuning is too narrow to reach. - assert not np.isnan(masked["ena_intensity"].values).any() + # ecliptic, which the pivot 75 tuning is too narrow to reach. The counts + # are blanked by the mask alone, never by the exposure. + assert not is_fill(masked["ena_count"].values).any() def test_the_outlier_tail_of_a_level_is_masked( self, one_pointing, anc_dependencies @@ -1556,15 +1595,16 @@ def test_the_outlier_tail_of_a_level_is_masked( ) # The pivot 105 tuning cannot mask a band, so only the top half of each - # level's own intensity distribution is masked. + # level's own intensity distribution is masked. The unexposed pixels are + # filled either way. intensity = raw["ena_intensity"].values - median = np.nanpercentile(intensity, 50, axis=(-2, -1), keepdims=True) - expected = intensity > median + unexposed = is_fill(intensity) + as_masked = np.where(unexposed, 0.0, intensity) + median = np.percentile(as_masked, 50, axis=(-2, -1), keepdims=True) + expected = (as_masked > median) | unexposed assert expected.any(), "the tuning must mask something" - np.testing.assert_array_equal( - np.isnan(masked["ena_intensity"].values), expected - ) + np.testing.assert_array_equal(is_fill(masked["ena_intensity"].values), expected) def test_an_untuned_pivot_angle_is_refused(self, one_pointing, anc_dependencies): """A map cannot be masked at a pivot the ancillary says nothing about.""" From 5e8b085ff0ac59b2b935562dbf6bf24f5031aefe Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 10 Aug 2026 10:42:28 -0400 Subject: [PATCH 13/14] merged dev --- imap_processing/lo/l2/lo_l2.py | 140 +++++++++++++++++++-------------- 1 file changed, 83 insertions(+), 57 deletions(-) diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index c5f29d2ba..d24475fde 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -1615,7 +1615,7 @@ def _calculate_rates_and_intensities( bootstrap_matrix: np.ndarray | None = None, flux_corrector: PowerLawFluxCorrector | None = None, isn_mask_parameters: pd.DataFrame | None = None, -) -> dict[str, np.ndarray]: +) -> dict[str, xr.DataArray]: """ Turn the accumulated counts and exposure into rates and intensities. @@ -1648,50 +1648,51 @@ def _calculate_rates_and_intensities( Returns ------- - dict[str, np.ndarray] + dict[str, xr.DataArray] The map variables, each of shape (epoch, esa level, pixel). """ - counts = sky_map.data_1d["ena_count"].values - exposure = sky_map.data_1d["exposure_factor"].values - bg_rate_exposure = sky_map.data_1d["bg_rate_exposure"].values + counts = sky_map.data_1d["ena_count"] + exposure = sky_map.data_1d["exposure_factor"] + bg_rate_exposure = sky_map.data_1d["bg_rate_exposure"] if sputter_matrix is None: rate_counts, rate_counts_var = counts, counts else: - rate_counts, rate_counts_var = _sputter_correct_counts( - counts, sky_map.data_1d["sputter_source_count"].values, sputter_matrix + corrected_counts, corrected_variance = _sputter_correct_counts( + counts.values, + sky_map.data_1d["sputter_source_count"].values, + sputter_matrix, ) + rate_counts = counts.copy(data=corrected_counts) + rate_counts_var = counts.copy(data=corrected_variance) - # Every ESA level quantity gets a pixel axis to broadcast over the map. - energy = calibration.energy[:, np.newaxis] - geometric_factor = calibration.geometric_factor[:, np.newaxis] - gf_low = calibration.geometric_factor_low[:, np.newaxis] - gf_high = calibration.geometric_factor_high[:, np.newaxis] + # Naming the energy dimension lets xarray broadcast during division + # without the need to introduce new dimensions + energy_dim = CoordNames.ENERGY_L2.value + energy = xr.DataArray(calibration.energy, dims=[energy_dim]) + geometric_factor = xr.DataArray(calibration.geometric_factor, dims=[energy_dim]) + gf_low = xr.DataArray(calibration.geometric_factor_low, dims=[energy_dim]) + gf_high = xr.DataArray(calibration.geometric_factor_high, dims=[energy_dim]) exposed = exposure > 0 - def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: + def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray: """ Divide only where the map was exposed, zero elsewhere. Parameters ---------- - numerator : np.ndarray + numerator : xr.DataArray The array being divided. - denominator : np.ndarray + denominator : xr.DataArray The array to divide it by. Returns ------- - np.ndarray + xr.DataArray The quotient, zero in the pixels that were never exposed. """ - return np.divide( - numerator, - denominator, - out=np.zeros_like(exposure), - where=exposed, - ) + return (numerator / denominator).where(exposed, 0) # Removing the sputtered counts can take a low-count pixel below zero, # which is not a rate the instrument can have observed. @@ -1708,10 +1709,15 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: isn_mask = None if isn_mask_parameters is not None: logger.info("Masking the ISN band out of the map") - isn_mask = _isn_mask( - _divide(_divide(counts, exposure), geometric_factor * energy), - sky_map.az_el_points.values[:, 1], - isn_mask_parameters, + uncorrected_intensity = _divide( + _divide(counts, exposure), geometric_factor * energy + ) + isn_mask = counts.copy( + data=_isn_mask( + uncorrected_intensity.values, + sky_map.az_el_points.values[:, 1], + isn_mask_parameters, + ) ) # The systematic error is the intensity excursion from the recalibrated @@ -1722,27 +1728,35 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: if not valid.all(): logger.warning( "The geometric factor of ESA levels " - f"{(np.flatnonzero(~valid[:, 0]) + 1).tolist()} is below its lower " + f"{(np.flatnonzero(~valid.values) + 1).tolist()} is below its lower " f"error bound; their systematic errors are left at zero." ) - intensity_upper = _divide(count_rate, np.where(valid, gf_low, 1.0) * energy) + intensity_upper = _divide(count_rate, gf_low.where(valid, 1.0) * energy) intensity_lower = _divide(count_rate, gf_high * energy) - intensity_sys_err_plus = np.where(valid, intensity_upper - intensity, 0.0) - intensity_sys_err_minus = np.where(valid, intensity - intensity_lower, 0.0) + intensity_sys_err_plus = (intensity_upper - intensity).where(valid, 0.0) + intensity_sys_err_minus = (intensity - intensity_lower).where(valid, 0.0) if bootstrap_matrix is not None: ( - intensity, - intensity_stat_uncert, - intensity_sys_err_plus, - intensity_sys_err_minus, + corrected, + corrected_stat_uncert, + corrected_sys_err_plus, + corrected_sys_err_minus, ) = _bootstrap_correct_intensity( - intensity, - intensity_stat_uncert**2, + intensity.values, + intensity_stat_uncert.values**2, calibration, bootstrap_matrix, sky_map.binning_grid_shape, - valid, + valid.values[:, np.newaxis], + ) + intensity = intensity.copy(data=corrected) + intensity_stat_uncert = intensity_stat_uncert.copy(data=corrected_stat_uncert) + intensity_sys_err_plus = intensity_sys_err_plus.copy( + data=corrected_sys_err_plus + ) + intensity_sys_err_minus = intensity_sys_err_minus.copy( + data=corrected_sys_err_minus ) bg_rate = _divide(bg_rate_exposure, exposure) @@ -1755,33 +1769,46 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: # of its own, so it is corrected on its own terms rather than with the # map's. if flux_corrector is not None: - cos_alpha = _divide(sky_map.data_1d["cos_alpha_exposure"].values, exposure) + cos_alpha = _divide(sky_map.data_1d["cos_alpha_exposure"], exposure) logger.info("Applying the Compton-Getting correction to the intensities") ( - intensity, + corrected, ( - intensity_stat_uncert, - intensity_sys_err_plus, - intensity_sys_err_minus, + corrected_stat_uncert, + corrected_sys_err_plus, + corrected_sys_err_minus, ), ) = _compton_getting_correct_intensity( - intensity, + intensity.values, ( - intensity_stat_uncert, - intensity_sys_err_plus, - intensity_sys_err_minus, + intensity_stat_uncert.values, + intensity_sys_err_plus.values, + intensity_sys_err_minus.values, ), - cos_alpha, + cos_alpha.values, calibration, flux_corrector, ) - bg_intensity, (bg_intensity_stat_uncert,) = _compton_getting_correct_intensity( - bg_intensity, - (bg_intensity_stat_uncert,), - cos_alpha, + intensity = intensity.copy(data=corrected) + intensity_stat_uncert = intensity_stat_uncert.copy(data=corrected_stat_uncert) + intensity_sys_err_plus = intensity_sys_err_plus.copy( + data=corrected_sys_err_plus + ) + intensity_sys_err_minus = intensity_sys_err_minus.copy( + data=corrected_sys_err_minus + ) + + bg_corrected, (bg_corrected_stat_uncert,) = _compton_getting_correct_intensity( + bg_intensity.values, + (bg_intensity_stat_uncert.values,), + cos_alpha.values, calibration, flux_corrector, ) + bg_intensity = bg_intensity.copy(data=bg_corrected) + bg_intensity_stat_uncert = bg_intensity_stat_uncert.copy( + data=bg_corrected_stat_uncert + ) variables = { "ena_count": counts, @@ -1804,17 +1831,17 @@ def _divide(numerator: np.ndarray, denominator: np.ndarray) -> np.ndarray: # A masked pixel has no ENA measurement to report. if isn_mask is not None: for name in ISN_MASKED_VARIABLES: - variables[name] = np.where(isn_mask, FILLVAL_FLOAT, variables[name]) + variables[name] = variables[name].where(~isn_mask, FILLVAL_FLOAT) for name in FILLED_VARIABLES: - variables[name] = np.where(exposed, variables[name], FILLVAL_FLOAT) + variables[name] = variables[name].where(exposed, FILLVAL_FLOAT) return variables def _build_map_dataset( sky_map: RectangularSkyMap, - variables: dict[str, np.ndarray], + variables: dict[str, xr.DataArray], calibration: EsaCalibration, ) -> xr.Dataset: """ @@ -1827,7 +1854,7 @@ def _build_map_dataset( ---------- sky_map : RectangularSkyMap The map being built. - variables : dict[str, np.ndarray] + variables : dict[str, xr.DataArray] The map variables, each of shape (epoch, esa level, pixel). calibration : EsaCalibration The energy response the map is binned in, read for the widths of the @@ -1839,9 +1866,8 @@ def _build_map_dataset( The map variables on the (epoch, energy, longitude, latitude) grid, with the energy coordinate and its widths. """ - dims = sky_map.data_1d["ena_count"].dims for name, values in variables.items(): - sky_map.data_1d[name] = xr.DataArray(values.astype(np.float32), dims=dims) + sky_map.data_1d[name] = values.astype(np.float32) # These are accumulators, not map variables. sky_map.data_1d = sky_map.data_1d.drop_vars( ["bg_rate_exposure", "sputter_source_count", "cos_alpha_exposure"], From 1d0336269f5d4cf6ba8271ad989fca0375ae745f Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 10 Aug 2026 11:04:07 -0400 Subject: [PATCH 14/14] explanatory note on einsums --- imap_processing/lo/l2/lo_l2.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/imap_processing/lo/l2/lo_l2.py b/imap_processing/lo/l2/lo_l2.py index d24475fde..4b2794197 100644 --- a/imap_processing/lo/l2/lo_l2.py +++ b/imap_processing/lo/l2/lo_l2.py @@ -1170,7 +1170,13 @@ def _sputter_correct_counts( square of its factor, and the variance only ever grows. """ logger.info("Applying the sputter correction to the accumulated counts") + # The einsum contracts the source level s away: target level t of the + # output is the sum over s of sputter_matrix[t, s] * source_counts[e, s, p], + # for every epoch e and pixel p. corrected = counts - np.einsum("ts,esp->etp", sputter_matrix, source_counts) + + # The same contraction against the squared matrix, which is how the + # variances of the scaled source terms add. variance = counts + np.einsum("ts,esp->etp", sputter_matrix**2, source_counts) return corrected, variance