diff --git a/imap_processing/ena_maps/utils/naming.py b/imap_processing/ena_maps/utils/naming.py index 4666e1afa..ada40d923 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: @@ -243,9 +249,7 @@ def _parse_principal_data(self) -> tuple[str, str]: extras is any extension on the principal data portion of the descriptor string. For example, "isnnbkgnd" would return ("isn", "nbkgnd"). """ - m = re.match( - r"^(drt|ena|int|isn|spx)(?:(?<=spx)\d+)?([^-_\s]*)$", self.principal_data - ) + m = PRINCIPAL_DATA_PATTERN.match(self.principal_data) if not m: raise ValueError( "Invalid principal_data format: " @@ -469,6 +473,96 @@ 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 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: + """ + 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") 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 + # 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/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/constants.py b/imap_processing/lo/constants.py index 87119c436..8d8c8b481 100644 --- a/imap_processing/lo/constants.py +++ b/imap_processing/lo/constants.py @@ -96,6 +96,36 @@ 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 + + # 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 4d2a630f6..4b2794197 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 ( @@ -14,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 @@ -38,6 +40,38 @@ # 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 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" @@ -63,7 +97,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 ---------- @@ -88,22 +125,21 @@ 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, 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 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 + # 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 ) logger.info("Step 1: Loading ancillary data") @@ -123,20 +159,60 @@ 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 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 map_descriptor.sputter_corrected + else (None, None) + ) + accumulators = ( + ACCUMULATED_VARIABLES + + (("sputter_source_count",) if sputter_source else ()) + + (("cos_alpha_exposure",) if map_descriptor.cg_corrected 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) + 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, + isn_mask_parameters, + ) dataset = _build_map_dataset(sky_map, variables, calibration) logger.info("IMAP-Lo L2 processing pipeline completed successfully") @@ -150,85 +226,6 @@ 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 # ============================================================================= @@ -266,40 +263,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: @@ -312,17 +341,199 @@ 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 _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_angles: list[int]) -> pd.DataFrame: + """ + 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_angles : list[int] + The nominal pivot angles [degrees] the map was built from. + + Returns + ------- + pd.DataFrame + The tuning of each ESA level, indexed by 1-based ESA step. + + Raises + ------ + ValueError + If the ancillary has no tuning for one of the pivot angles. + """ + parameters = load_isn_mask_parameters() + + 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 {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 tuning.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: """ @@ -373,6 +584,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]: @@ -517,7 +809,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 +826,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 +848,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 +869,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 +900,31 @@ 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) + ) + 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. + values["cos_alpha_exposure"] = ( + _ram_projection(spin_angles, pivot_angle) * pointing_exposure + ) + 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 +932,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, ) @@ -682,10 +995,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 # ============================================================================= @@ -807,13 +1145,488 @@ 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") + # 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 + + +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. + + 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 + 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 _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 + sky_map: RectangularSkyMap, + calibration: EsaCalibration, + 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, xr.DataArray]: """ 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 ---------- @@ -822,6 +1635,22 @@ 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. + 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. + 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. + 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 ------- @@ -832,6 +1661,17 @@ def _calculate_rates_and_intensities( 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: + 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) + # Naming the energy dimension lets xarray broadcast during division # without the need to introduce new dimensions energy_dim = CoordNames.ENERGY_L2.value @@ -860,13 +1700,32 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray: """ return (numerator / denominator).where(exposed, 0) - 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) + # 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") + 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 # 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 @@ -883,12 +1742,81 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray: 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: + ( + corrected, + corrected_stat_uncert, + corrected_sys_err_plus, + corrected_sys_err_minus, + ) = _bootstrap_correct_intensity( + intensity.values, + intensity_stat_uncert.values**2, + calibration, + bootstrap_matrix, + sky_map.binning_grid_shape, + 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) bg_rate_stat_uncert = np.sqrt(_divide(bg_rate, exposure)) bg_intensity = _divide(bg_rate, geometric_factor * energy) bg_intensity_stat_uncert = _divide(bg_rate_stat_uncert, geometric_factor * energy) - return { + # 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"], exposure) + logger.info("Applying the Compton-Getting correction to the intensities") + ( + corrected, + ( + corrected_stat_uncert, + corrected_sys_err_plus, + corrected_sys_err_minus, + ), + ) = _compton_getting_correct_intensity( + intensity.values, + ( + intensity_stat_uncert.values, + intensity_sys_err_plus.values, + intensity_sys_err_minus.values, + ), + cos_alpha.values, + calibration, + flux_corrector, + ) + 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, "exposure_factor": exposure, "ena_count_rate": count_rate, @@ -906,6 +1834,16 @@ def _divide(numerator: xr.DataArray, denominator: xr.DataArray) -> xr.DataArray: "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] = variables[name].where(~isn_mask, FILLVAL_FLOAT) + + for name in FILLED_VARIABLES: + variables[name] = variables[name].where(exposed, FILLVAL_FLOAT) + + return variables + def _build_map_dataset( sky_map: RectangularSkyMap, @@ -936,8 +1874,11 @@ def _build_map_dataset( """ for name, values in variables.items(): sky_map.data_1d[name] = values.astype(np.float32) - # `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", "cos_alpha_exposure"], + 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 6dcaf856e..29a6b9fdb 100644 --- a/imap_processing/tests/ena_maps/test_naming.py +++ b/imap_processing/tests/ena_maps/test_naming.py @@ -610,3 +610,114 @@ def test_get_duration_str_generic_unit_full(self): duration="5day", ) assert md._get_duration_str(full=True) == "5 day" + + @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, + # 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), + ("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 + + @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 + + @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_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_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_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 51321bd05..46ef837ae 100644 --- a/imap_processing/tests/lo/test_lo_l2.py +++ b/imap_processing/tests/lo/test_lo_l2.py @@ -9,29 +9,85 @@ 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 ( 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, _complete_pointings, + _compton_getting_correct_intensity, _dps_spin_angles, + _esa_calibration, + _extrapolate_top_intensity, + _spacecraft_frame_energy, _spin_phase_mask, 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 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" 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 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 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 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 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 + +# 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 @@ -58,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(): @@ -104,6 +174,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)), }, @@ -189,8 +263,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"] @@ -219,6 +293,295 @@ 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 + + +@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 + + +@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 + + +@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 + + +@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 = 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. + 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) + + 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 is_fill(combined["ena_intensity"].values).any() + np.testing.assert_array_equal( + is_fill(combined["ena_intensity"].values), + is_fill(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.""" @@ -369,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] @@ -407,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.""" @@ -429,6 +792,560 @@ 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(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.""" + 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 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 + + # 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 + ) + 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.""" + 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 + + 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 + ) + 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(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.""" + 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(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 + ) + # 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 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 + + 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_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 ( + "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(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 + + intensity = raw["ena_intensity"].values + dark = (intensity <= 0) & ~is_fill(intensity) + 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.""" @@ -529,14 +1446,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() @@ -544,7 +1461,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 @@ -569,6 +1486,160 @@ 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 + # 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(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 = ~is_fill(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, raw = masked_maps + + # 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: + 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.""" + 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 + + # 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: + 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 + ): + """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. 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 + ): + """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. The unexposed pixels are + # filled either way. + intensity = raw["ena_intensity"].values + 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(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.""" + with pytest.raises(ValueError, match=r"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.""" diff --git a/imap_processing/tests/test_cli.py b/imap_processing/tests/test_cli.py index e7a005bee..9bd1aa24c 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(